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