Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-2019 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    plumed is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "core/ActionAtomistic.h"
      23             : #include "core/ActionPilot.h"
      24             : #include "core/ActionRegister.h"
      25             : #include "tools/Vector.h"
      26             : #include "tools/Matrix.h"
      27             : #include "tools/AtomNumber.h"
      28             : #include "tools/Tools.h"
      29             : #include "core/Atoms.h"
      30             : #include "tools/Pbc.h"
      31             : 
      32             : #include <vector>
      33             : #include <string>
      34             : 
      35             : using namespace std;
      36             : 
      37             : namespace PLMD {
      38             : namespace generic {
      39             : 
      40             : //+PLUMEDOC GENERIC RESET_CELL
      41             : /*
      42             : This action is used to rotate the full cell
      43             : 
      44             : This can be used to modify the periodic box. Notice that
      45             : this is done at fixed scaled coordinates,
      46             : so that also atomic coordinates for the entire system are affected.
      47             : To see what effect try
      48             : the \ref DUMPATOMS directive to output the atomic positions.
      49             : 
      50             : Also notice that PLUMED propagate forces correctly so that you can add a bias on a CV computed
      51             : after rotation. See also \ref FIT_TO_TEMPLATE
      52             : 
      53             : Currently, only TYPE=TRIANGULAR is implemented, which allows one to reset
      54             : the cell to a lower triangular one. Namely, a proper rotation is found that allows
      55             : rotating the box so that the first lattice vector is in the form (ax,0,0),
      56             : the second lattice vector is in the form (bx,by,0), and the third lattice vector is
      57             : arbitrary.
      58             : 
      59             : \attention
      60             : The implementation of this action is available but should be considered in testing phase. Please report any
      61             : strange behavior.
      62             : 
      63             : \attention
      64             : This directive modifies the stored position at the precise moment
      65             : it is executed. This means that only collective variables
      66             : which are below it in the input script will see the corrected positions.
      67             : Unless you
      68             : know exactly what you are doing, leave the default stride (1), so that
      69             : this action is performed at every MD step.
      70             : 
      71             : \par Examples
      72             : 
      73             : Reset cell to be triangular after a rototranslational fit
      74             : \plumedfile
      75             : DUMPATOMS FILE=dump-original.xyz ATOMS=1-20
      76             : FIT_TO_TEMPLATE STRIDE=1 REFERENCE=ref.pdb TYPE=OPTIMAL
      77             : DUMPATOMS FILE=dump-fit.xyz ATOMS=1-20
      78             : RESET_CELL TYPE=TRIANGULAR
      79             : DUMPATOMS FILE=dump-reset.xyz ATOMS=1-20
      80             : \endplumedfile
      81             : 
      82             : The reference file for the FIT_TO_TEMPLATE is just a normal pdb file with the format shown below:
      83             : 
      84             : \auxfile{ref.pdb}
      85             : ATOM      8  HT3 ALA     2      -1.480  -1.560   1.212  1.00  1.00      DIA  H
      86             : ATOM      9  CAY ALA     2      -0.096   2.144  -0.669  1.00  1.00      DIA  C
      87             : ATOM     10  HY1 ALA     2       0.871   2.385  -0.588  1.00  1.00      DIA  H
      88             : ATOM     12  HY3 ALA     2      -0.520   2.679  -1.400  1.00  1.00      DIA  H
      89             : ATOM     14  OY  ALA     2      -1.139   0.931  -0.973  1.00  1.00      DIA  O
      90             : END
      91             : \endauxfile
      92             : 
      93             : */
      94             : //+ENDPLUMEDOC
      95             : 
      96             : 
      97           8 : class ResetCell:
      98             :   public ActionPilot,
      99             :   public ActionAtomistic
     100             : {
     101             :   std::string type;
     102             :   Tensor rotation,newbox;
     103             : 
     104             : public:
     105             :   explicit ResetCell(const ActionOptions&ao);
     106             :   static void registerKeywords( Keywords& keys );
     107             :   void calculate() override;
     108             :   void apply() override;
     109             : };
     110             : 
     111        7836 : PLUMED_REGISTER_ACTION(ResetCell,"RESET_CELL")
     112             : 
     113           3 : void ResetCell::registerKeywords( Keywords& keys ) {
     114           3 :   Action::registerKeywords( keys );
     115           3 :   ActionAtomistic::registerKeywords( keys );
     116          15 :   keys.add("compulsory","STRIDE","1","the frequency with which molecules are reassembled.  Unless you are completely certain about what you are doing leave this set equal to 1!");
     117          15 :   keys.add("compulsory","TYPE","TRIANGULAR","the manner in which the cell is reset");
     118           3 : }
     119             : 
     120           2 : ResetCell::ResetCell(const ActionOptions&ao):
     121             :   Action(ao),
     122             :   ActionPilot(ao),
     123           4 :   ActionAtomistic(ao)
     124             : {
     125           2 :   type.assign("TRIANGULAR");
     126           4 :   parse("TYPE",type);
     127             : 
     128           2 :   log<<"  type: "<<type<<"\n";
     129           2 :   if(type!="TRIANGULAR") error("undefined type "+type);
     130             : 
     131           2 :   checkRead();
     132           2 : }
     133             : 
     134             : 
     135          17 : void ResetCell::calculate() {
     136             : 
     137             :   Pbc & pbc(modifyGlobalPbc());
     138             : 
     139          17 :   Tensor box=pbc.getBox();
     140             : 
     141             : // moduli of lattice vectors
     142          17 :   double a=modulo(box.getRow(0));
     143          17 :   double b=modulo(box.getRow(1));
     144          17 :   double c=modulo(box.getRow(2));
     145             : // cos-angle between lattice vectors
     146          17 :   double ab=dotProduct(box.getRow(0),box.getRow(1))/(a*b);
     147          17 :   double ac=dotProduct(box.getRow(0),box.getRow(2))/(a*c);
     148          17 :   double bc=dotProduct(box.getRow(1),box.getRow(2))/(b*c);
     149             : 
     150             : // generate a new set of lattice vectors as a lower triangular matrix
     151          17 :   newbox[0][0]=a;
     152          17 :   newbox[1][0]=b*ab;
     153          17 :   newbox[1][1]=std::sqrt(b*b-newbox[1][0]*newbox[1][0]);
     154          17 :   newbox[2][0]=c*ac;
     155          17 :   newbox[2][1]=c*(bc-ac*ab)/std::sqrt(1-ab*ab);
     156          17 :   newbox[2][2]=std::sqrt(c*c-newbox[2][0]*newbox[2][0]-newbox[2][1]*newbox[2][1]);
     157             : 
     158          17 :   if(determinant(newbox)*determinant(box)<0) newbox[2][2]=-newbox[2][2];
     159             : 
     160             : // rotation matrix from old to new coordinates
     161          17 :   rotation=transpose(matmul(inverse(box),newbox));
     162             : 
     163             : // rotate all coordinates
     164        3246 :   for(unsigned i=0; i<getTotAtoms(); i++) {
     165             :     Vector & ato (modifyGlobalPosition(AtomNumber::index(i)));
     166        1606 :     ato=matmul(rotation,ato);
     167             :   }
     168             : // rotate box
     169          17 :   pbc.setBox(newbox);
     170          17 : }
     171             : 
     172          17 : void ResetCell::apply() {
     173             : // rotate back forces
     174        3246 :   for(unsigned i=0; i<getTotAtoms(); i++) {
     175             :     Vector & f(modifyGlobalForce(AtomNumber::index(i)));
     176        1606 :     f=matmul(transpose(rotation),f);
     177             :   }
     178             : 
     179          17 :   Tensor& virial(modifyGlobalVirial());
     180             : // I have no mathematical derivation for this.
     181             : // The reasoning is the following.
     182             : // virial= h^T * dU/dh, where h is the box matrix and dU/dh its derivatives.
     183             : // The final virial should be rotationally invariant, that is symmetric.
     184             : // in the rotated frame, dU/dh elements [0][1], [0][2], and [1][2] should
     185             : // be changed so as to enforce rotational invariance. Thus we here have to
     186             : // make the virial matrix symmetric.
     187             : // Since h^T is upper triangular, it can be shown that any change in these elements
     188             : // will only affect the corresponding elements of the virial matrix.
     189             : // Thus, the only possibility is to set the corresponding elements
     190             : // of the virial matrix equal to their symmetric ones.
     191             : // GB
     192          17 :   virial[0][1]=virial[1][0];
     193          17 :   virial[0][2]=virial[2][0];
     194          17 :   virial[1][2]=virial[2][1];
     195             : // rotate back virial
     196          17 :   virial=matmul(transpose(rotation),matmul(virial,rotation));
     197             : 
     198             : 
     199             : 
     200          17 : }
     201             : 
     202             : }
     203        5874 : }