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             : Rotating the full cell is useful if you want 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 using [DUMPATOMS](DUMPATOMS.md) directive to output the atomic positions.
      45             : 
      46             : Also notice that PLUMED propagate forces correctly so that you can add a bias on a CV computed
      47             : after rotation. You can read the documentation for [FIT_TO_TEMPLATE](FIT_TO_TEMPLATE.md) for more
      48             : detail.
      49             : 
      50             : Currently, only TYPE=TRIANGULAR is implemented, which allows one to reset
      51             : the cell to a lower triangular one. This command finds a proper rotation
      52             : rotates 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             : !!! caution ""
      57             : 
      58             :     The implementation of this action is available but should be considered in testing phase. Please report any
      59             :     strange behavior.
      60             : 
      61             : !!! caution ""
      62             : 
      63             :     This directive modifies the stored position at the precise moment
      64             :     it is executed. This means that only collective variables
      65             :     which are below it in the input script will see the corrected positions.
      66             :     Unless you
      67             :     know exactly what you are doing, leave the default stride (1), so that
      68             :     this action is performed at every MD step.
      69             : 
      70             : ## Examples
      71             : 
      72             : Reset cell to be triangular after a rototranslational fit
      73             : 
      74             : ```plumed
      75             : #SETTINGS INPUTFILES=regtest/basic/rt63/align.pdb
      76             : DUMPATOMS FILE=dump-original.xyz ATOMS=1-20
      77             : FIT_TO_TEMPLATE STRIDE=1 REFERENCE=regtest/basic/rt63/align.pdb TYPE=OPTIMAL
      78             : DUMPATOMS FILE=dump-fit.xyz ATOMS=1-20
      79             : RESET_CELL TYPE=TRIANGULAR STRIDE=1
      80             : DUMPATOMS FILE=dump-reset.xyz ATOMS=1-20
      81             : ```
      82             : 
      83             : */
      84             : //+ENDPLUMEDOC
      85             : 
      86             : 
      87             : class ResetCell:
      88             :   public ActionPilot,
      89             :   public ActionAtomistic {
      90             :   std::string type;
      91             :   Tensor rotation,newbox;
      92             :   Value* boxValue;
      93             :   PbcAction* pbc_action;
      94             : public:
      95             :   explicit ResetCell(const ActionOptions&ao);
      96             :   static void registerKeywords( Keywords& keys );
      97             :   void calculate() override;
      98             :   void apply() override;
      99             : };
     100             : 
     101             : PLUMED_REGISTER_ACTION(ResetCell,"RESET_CELL")
     102             : 
     103           6 : void ResetCell::registerKeywords( Keywords& keys ) {
     104           6 :   Action::registerKeywords( keys );
     105           6 :   ActionAtomistic::registerKeywords( keys );
     106           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!");
     107           6 :   keys.add("compulsory","TYPE","TRIANGULAR","the manner in which the cell is reset");
     108           6 : }
     109             : 
     110           2 : ResetCell::ResetCell(const ActionOptions&ao):
     111             :   Action(ao),
     112             :   ActionPilot(ao),
     113           2 :   ActionAtomistic(ao) {
     114           2 :   type.assign("TRIANGULAR");
     115           2 :   parse("TYPE",type);
     116             : 
     117           2 :   log<<"  type: "<<type<<"\n";
     118           2 :   if(type!="TRIANGULAR") {
     119           0 :     error("undefined type "+type);
     120             :   }
     121             : 
     122           2 :   pbc_action=plumed.getActionSet().selectWithLabel<PbcAction*>("Box");
     123           2 :   if( !pbc_action ) {
     124           0 :     error("cannot reset cell if box has not been set");
     125             :   }
     126           2 :   boxValue=pbc_action->copyOutput(0);
     127           2 : }
     128             : 
     129             : 
     130          17 : void ResetCell::calculate() {
     131          17 :   Tensor box=pbc_action->getPbc().getBox();
     132             : 
     133             : // moduli of lattice vectors
     134          17 :   double a=modulo(box.getRow(0));
     135          17 :   double b=modulo(box.getRow(1));
     136          17 :   double c=modulo(box.getRow(2));
     137             : // cos-angle between lattice vectors
     138          17 :   double ab=dotProduct(box.getRow(0),box.getRow(1))/(a*b);
     139          17 :   double ac=dotProduct(box.getRow(0),box.getRow(2))/(a*c);
     140          17 :   double bc=dotProduct(box.getRow(1),box.getRow(2))/(b*c);
     141             : 
     142             : // generate a new set of lattice vectors as a lower triangular matrix
     143          17 :   newbox[0][0]=a;
     144          17 :   newbox[1][0]=b*ab;
     145          17 :   newbox[1][1]=std::sqrt(b*b-newbox[1][0]*newbox[1][0]);
     146          17 :   newbox[2][0]=c*ac;
     147          17 :   newbox[2][1]=c*(bc-ac*ab)/std::sqrt(1-ab*ab);
     148          17 :   newbox[2][2]=std::sqrt(c*c-newbox[2][0]*newbox[2][0]-newbox[2][1]*newbox[2][1]);
     149             : 
     150          17 :   if(determinant(newbox)*determinant(box)<0) {
     151           1 :     newbox[2][2]=-newbox[2][2];
     152             :   }
     153             : 
     154             : // rotation matrix from old to new coordinates
     155          17 :   rotation=transpose(matmul(inverse(box),newbox));
     156             : 
     157             : // rotate all coordinates
     158          17 :   rotateAllGlobalPositions(rotation);
     159             : // rotate box
     160          17 :   pbc_action->getPbc().setBox(newbox);
     161          17 : }
     162             : 
     163          17 : void ResetCell::apply() {
     164             : // rotate back forces
     165          17 :   const Tensor rotationT=transpose(rotation);
     166          17 :   rotateForcesOnAllAtoms(rotationT);
     167             : 
     168             : // I have no mathematical derivation for this.
     169             : // The reasoning is the following.
     170             : // virial= h^T * dU/dh, where h is the box matrix and dU/dh its derivatives.
     171             : // The final virial should be rotationally invariant, that is symmetric.
     172             : // in the rotated frame, dU/dh elements [0][1], [0][2], and [1][2] should
     173             : // be changed so as to enforce rotational invariance. Thus we here have to
     174             : // make the virial matrix symmetric.
     175             : // Since h^T is upper triangular, it can be shown that any change in these elements
     176             : // will only affect the corresponding elements of the virial matrix.
     177             : // Thus, the only possibility is to set the corresponding elements
     178             : // of the virial matrix equal to their symmetric ones.
     179             : // GB
     180          17 :   Tensor virial;
     181          68 :   for(unsigned i=0; i<3; ++i)
     182         204 :     for(unsigned j=0; j<3; ++j) {
     183         153 :       virial[i][j]=boxValue->getForce( 3*i+j );
     184             :     }
     185          17 :   virial[0][1]=virial[1][0];
     186          17 :   virial[0][2]=virial[2][0];
     187          17 :   virial[1][2]=virial[2][1];
     188             : // rotate back virial
     189          17 :   virial=matmul(rotationT,matmul(virial,rotation));
     190          17 :   boxValue->clearInputForce();
     191          68 :   for(unsigned i=0; i<3; ++i)
     192         204 :     for(unsigned j=0; j<3; ++j) {
     193         153 :       boxValue->addForce( 3*i+j, virial(i,j) );
     194             :     }
     195             : 
     196             : 
     197          17 : }
     198             : 
     199             : }
     200             : }