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             : {
      96             :   std::string type;
      97             :   Tensor rotation,newbox;
      98             : 
      99             : public:
     100             :   explicit ResetCell(const ActionOptions&ao);
     101             :   static void registerKeywords( Keywords& keys );
     102             :   void calculate() override;
     103             :   void apply() override;
     104             : };
     105             : 
     106       10423 : PLUMED_REGISTER_ACTION(ResetCell,"RESET_CELL")
     107             : 
     108           3 : void ResetCell::registerKeywords( Keywords& keys ) {
     109           3 :   Action::registerKeywords( keys );
     110           3 :   ActionAtomistic::registerKeywords( keys );
     111           6 :   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!");
     112           6 :   keys.add("compulsory","TYPE","TRIANGULAR","the manner in which the cell is reset");
     113           3 : }
     114             : 
     115           2 : ResetCell::ResetCell(const ActionOptions&ao):
     116             :   Action(ao),
     117             :   ActionPilot(ao),
     118           2 :   ActionAtomistic(ao)
     119             : {
     120           2 :   type.assign("TRIANGULAR");
     121           2 :   parse("TYPE",type);
     122             : 
     123           2 :   log<<"  type: "<<type<<"\n";
     124           2 :   if(type!="TRIANGULAR") error("undefined type "+type);
     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) newbox[2][2]=-newbox[2][2];
     154             : 
     155             : // rotation matrix from old to new coordinates
     156          17 :   rotation=transpose(matmul(inverse(box),newbox));
     157             : 
     158             : // rotate all coordinates
     159        1623 :   for(unsigned i=0; i<getTotAtoms(); i++) {
     160             :     Vector & ato (modifyGlobalPosition(AtomNumber::index(i)));
     161        1606 :     ato=matmul(rotation,ato);
     162             :   }
     163             : // rotate box
     164          17 :   pbc.setBox(newbox);
     165          17 : }
     166             : 
     167          17 : void ResetCell::apply() {
     168             : // rotate back forces
     169        1623 :   for(unsigned i=0; i<getTotAtoms(); i++) {
     170             :     Vector & f(modifyGlobalForce(AtomNumber::index(i)));
     171        1606 :     f=matmul(transpose(rotation),f);
     172             :   }
     173             : 
     174             :   Tensor& virial(modifyGlobalVirial());
     175             : // I have no mathematical derivation for this.
     176             : // The reasoning is the following.
     177             : // virial= h^T * dU/dh, where h is the box matrix and dU/dh its derivatives.
     178             : // The final virial should be rotationally invariant, that is symmetric.
     179             : // in the rotated frame, dU/dh elements [0][1], [0][2], and [1][2] should
     180             : // be changed so as to enforce rotational invariance. Thus we here have to
     181             : // make the virial matrix symmetric.
     182             : // Since h^T is upper triangular, it can be shown that any change in these elements
     183             : // will only affect the corresponding elements of the virial matrix.
     184             : // Thus, the only possibility is to set the corresponding elements
     185             : // of the virial matrix equal to their symmetric ones.
     186             : // GB
     187          17 :   virial[0][1]=virial[1][0];
     188          17 :   virial[0][2]=virial[2][0];
     189          17 :   virial[1][2]=virial[2][1];
     190             : // rotate back virial
     191          17 :   virial=matmul(transpose(rotation),matmul(virial,rotation));
     192             : 
     193             : 
     194             : 
     195          17 : }
     196             : 
     197             : }
     198             : }