Line data    Source code

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