Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-2020 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/ActionShortcut.h"
      23             : #include "core/ActionRegister.h"
      24             : #include "core/PlumedMain.h"
      25             : #include "tools/PDB.h"
      26             : #include "Path.h"
      27             : 
      28             : //+PLUMEDOC COLVAR ADAPTIVE_PATH
      29             : /*
      30             : Compute path collective variables that adapt to the lowest free energy path connecting states A and B.
      31             : 
      32             : The Path Collective Variables developed by Branduardi and co-workers \cite brand07 allow one
      33             : to compute the progress along a high-dimensional path and the distance from the high-dimensional
      34             : path.  The progress along the path (s) is computed using:
      35             : 
      36             : \f[
      37             : s = i_2 + \textrm{sign}(i_2-i_1) \frac{ \sqrt{( \mathbf{v}_1\cdot\mathbf{v}_2 )^2 - |\mathbf{v}_3|^2(|\mathbf{v}_1|^2 - |\mathbf{v}_2|^2) } }{2|\mathbf{v}_3|^2} - \frac{\mathbf{v}_1\cdot\mathbf{v}_3 - |\mathbf{v}_3|^2}{2|\mathbf{v}_3|^2}
      38             : \f]
      39             : 
      40             : In this expression \f$\mathbf{v}_1\f$ and \f$\mathbf{v}_3\f$ are the vectors connecting the current position to the closest and second closest node of the path,
      41             : respectfully and \f$i_1\f$ and \f$i_2\f$ are the projections of the closest and second closest frames of the path. \f$\mathbf{v}_2\f$, meanwhile, is the
      42             : vector connecting the closest frame to the second closest frame.  The distance from the path, \f$z\f$ is calculated using:
      43             : 
      44             : \f[
      45             : z = \sqrt{ \left[ |\mathbf{v}_1|^2 - |\mathbf{v}_2| \left( \frac{ \sqrt{( \mathbf{v}_1\cdot\mathbf{v}_2 )^2 - |\mathbf{v}_3|^2(|\mathbf{v}_1|^2 - |\mathbf{v}_2|^2) } }{2|\mathbf{v}_3|^2} - \frac{\mathbf{v}_1\cdot\mathbf{v}_3 - |\mathbf{v}_3|^2}{2|\mathbf{v}_3|^2} \right) \right]^2 }
      46             : \f]
      47             : 
      48             : Notice that these are the definitions of \f$s\f$ and \f$z\f$ that are used by \ref PATH when the GPATH option is employed.  The reason for this is that
      49             : the adaptive path method implemented in this action was inspired by the work of Diaz and Ensing in which these formula were used \cite BerndAdaptivePath.
      50             : To learn more about how the path is adapted we strongly recommend reading this paper.
      51             : 
      52             : \par Examples
      53             : 
      54             : The input below provides an example that shows how the adaptive path works. The path is updated every 50 steps of
      55             : MD based on the data accumulated during the preceding 50 time steps.
      56             : 
      57             : \plumedfile
      58             : d1: DISTANCE ATOMS=1,2 COMPONENTS
      59             : pp: ADAPTIVE_PATH TYPE=EUCLIDEAN FIXED=2,5 UPDATE=50 WFILE=out-path.pdb WSTRIDE=50 REFERENCE=mypath.pdb
      60             : PRINT ARG=d1.x,d1.y,pp.* FILE=colvar
      61             : \endplumedfile
      62             : 
      63             : In the case above the distance between frames is calculated based on the \f$x\f$ and \f$y\f$ components of the vector connecting
      64             : atoms 1 and 2.  As such an extract from the input reference path (mypath.pdb) would look as follows:
      65             : 
      66             : \auxfile{mypath.pdb}
      67             : REMARK ARG=d1.x,d1.y d1.x=1.12 d1.y=-.60
      68             : END
      69             : REMARK ARG=d1.x,d1.y d1.x=.99 d1.y=-.45
      70             : END
      71             : REMARK ARG=d1.x,d1.y d1.x=.86 d1.y=-.30
      72             : END
      73             : REMARK ARG=d1.x,d1.y d1.x=.73 d1.y=-.15
      74             : END
      75             : REMARK ARG=d1.x,d1.y d1.x=.60 d1.y=0
      76             : END
      77             : REMARK ARG=d1.x,d1.y d1.x=.47 d1.y=.15
      78             : END
      79             : \endauxfile
      80             : 
      81             : Notice that one can also use RMSD frames in place of arguments like those above.
      82             : 
      83             : */
      84             : //+ENDPLUMEDOC
      85             : 
      86             : namespace PLMD {
      87             : namespace mapping {
      88             : 
      89             : class AdaptivePath : public ActionShortcut {
      90             : public:
      91             :   static void registerKeywords( Keywords& keys );
      92             :   explicit AdaptivePath(const ActionOptions&);
      93             : };
      94             : 
      95             : PLUMED_REGISTER_ACTION(AdaptivePath,"ADAPTIVE_PATH")
      96             : 
      97           6 : void AdaptivePath::registerKeywords( Keywords& keys ) {
      98           6 :   ActionShortcut::registerKeywords( keys ); Path::registerInputFileKeywords( keys );
      99          12 :   keys.add("optional","PROPERTY","read in path coordinates by finding option with this label in remark of pdb frames");
     100          12 :   keys.add("compulsory","FIXED","the positions in the list of input frames of the two path nodes whose positions remain fixed during the path optimization");
     101          12 :   keys.add("compulsory","HALFLIFE","-1","the number of MD steps after which a previously measured path distance weighs only 50 percent in the average. This option may increase convergence by allowing to forget the memory of a bad initial guess path. The default is to set this to infinity");
     102          12 :   keys.add("compulsory","UPDATE","the frequency with which the path should be updated");
     103          12 :   keys.add("compulsory","TOLERANCE","1E-6","the tolerance to use for the path updating algorithm that makes all frames equidistant");
     104          12 :   keys.add("optional","WFILE","file on which to write out the path");
     105          12 :   keys.add("compulsory","FMT","%f","the format to use for output files");
     106          12 :   keys.add("compulsory","WSTRIDE","frequency with which to write out the path");
     107           6 :   keys.setValueDescription("the position along and from the adaptive path");
     108          12 :   keys.needsAction("GEOMETRIC_PATH"); keys.needsAction("AVERAGE_PATH_DISPLACEMENT");
     109          12 :   keys.needsAction("REPARAMETERIZE_PATH"); keys.needsAction("DUMPPDB");
     110          18 :   keys.needsAction("PDB2CONSTANT"); keys.needsAction("DISPLACEMENT"); keys.needsAction("CONSTANT");
     111           6 : }
     112             : 
     113           2 : AdaptivePath::AdaptivePath(const ActionOptions& ao):
     114             :   Action(ao),
     115           2 :   ActionShortcut(ao)
     116             : {
     117             :   // Read in the arguments
     118           4 :   std::vector<std::string> argnames; parseVector("ARG",argnames);
     119           4 :   std::string reference_data, metric, mtype; parse("TYPE", mtype);
     120           4 :   std::string reference; parse("REFERENCE",reference);
     121           2 :   FILE* fp=std::fopen(reference.c_str(),"r"); PDB mypdb; if(!fp) error("could not open reference file " + reference );
     122           2 :   bool do_read=mypdb.readFromFilepointer(fp,false,0.1); if( !do_read ) error("missing file " + reference );
     123             :   // Create list of reference configurations that PLUMED will use
     124           2 :   Path::readInputFrames( reference, mtype, argnames, true, this, reference_data );
     125             :   // Now get coordinates on spath
     126           4 :   std::vector<std::string> pnames; parseVector("PROPERTY",pnames); Path::readPropertyInformation( pnames, getShortcutLabel(), reference, this );
     127             :   // Create action that computes the geometric path variablesa
     128           3 :   std::string propstr = getShortcutLabel() + "_ind"; if( pnames.size()>0 ) propstr = pnames[0] + "_ref";
     129           3 :   if( argnames.size()>0 ) readInputLine( getShortcutLabel() + ": GEOMETRIC_PATH ARG=" + getShortcutLabel() + "_data " + " PROPERTY=" + propstr + " REFERENCE=" + reference_data + " METRIC={DIFFERENCE}");
     130             :   else {
     131           1 :     std::string num, align_str, displace_str; Tools::convert( mypdb.getOccupancy()[0], align_str ); Tools::convert( mypdb.getBeta()[0], displace_str );
     132          25 :     for(unsigned j=1; j<mypdb.getAtomNumbers().size(); ++j ) { Tools::convert( mypdb.getOccupancy()[j], num ); align_str += "," + num; Tools::convert( mypdb.getBeta()[0], num ); displace_str += "," + num; }
     133           2 :     metric = "RMSD_VECTOR DISPLACEMENT TYPE=" + mtype + " ALIGN=" + align_str + " DISPLACE=" + displace_str;
     134           2 :     readInputLine( getShortcutLabel() + ": GEOMETRIC_PATH ARG=" + getShortcutLabel() + "_data.disp " + " PROPERTY=" +  propstr + " REFERENCE=" + reference_data + " METRIC={" + metric + "} METRIC_COMPONENT=disp");
     135             :   }
     136             :   // Create the object to accumulate the average path displacements
     137           8 :   std::string update, halflife; parse("HALFLIFE",halflife); parse("UPDATE",update); std::string refframes = " REFERENCE=" + getShortcutLabel() + "_pos";
     138           3 :   if( argnames.size()>0 ) readInputLine( getShortcutLabel() + "_disp: AVERAGE_PATH_DISPLACEMENT ARG=" + getShortcutLabel() + "_data HALFLIFE=" + halflife + " CLEAR=" + update + " METRIC={DIFFERENCE} REFERENCE=" + reference_data );
     139           2 :   else readInputLine( getShortcutLabel() + "_disp: AVERAGE_PATH_DISPLACEMENT ARG=" + getShortcutLabel() + "_data.disp HALFLIFE=" + halflife + " CLEAR=" + update + " METRIC={" + metric + "} METRIC_COMPONENT=disp REFERENCE=" + reference_data );
     140             : 
     141             :   // Create the object to update the path
     142           4 :   std::string fixedn; parse("FIXED",fixedn); std::string component="METRIC_COMPONENT=disp"; if( argnames.size()>0 ) { metric="DIFFERENCE"; component=""; }
     143           4 :   if( fixedn.length()>0 ) readInputLine("REPARAMETERIZE_PATH DISPLACE_FRAMES=" + getShortcutLabel() + "_disp FIXED=" + fixedn + " STRIDE=" + update + " METRIC={" + metric + "} " + component + " REFERENCE=" + reference_data );
     144           0 :   else readInputLine("REPARAMETERIZE_PATH DISPLACE_FRAMES=" + getShortcutLabel() + "_disp STRIDE=" + update + " METRIC={" + metric + "} " + component + " REFERENCE=" + reference_data );
     145             : 
     146             :   // Information for write out
     147           4 :   std::string wfilename; parse("WFILE",wfilename);
     148           2 :   if( wfilename.length()>0 ) {
     149             :     // This just gets the atom numbers for output
     150             :     std::string atomstr;
     151           2 :     if( argnames.size()==0 ) {
     152           1 :       FILE* fp=std::fopen(reference.c_str(),"r"); double fake_unit=0.1; PDB mypdb; bool do_read=mypdb.readFromFilepointer(fp,false,fake_unit);
     153           1 :       std::string num; Tools::convert( mypdb.getAtomNumbers()[0].serial(), atomstr );
     154          13 :       for(unsigned j=1; j<mypdb.getAtomNumbers().size(); ++j ) { Tools::convert( mypdb.getAtomNumbers()[j].serial(), num ); atomstr += "," + num; }
     155           1 :     }
     156             : 
     157           2 :     if( wfilename.find(".pdb")==std::string::npos ) error("output must be to a pdb file");
     158           6 :     std::string ofmt, pframes, wstride; parse("WSTRIDE",wstride); parse("FMT",ofmt);
     159           2 :     if( argnames.size()>0 ) {
     160           3 :       std::string argstr = argnames[0]; for(unsigned i=1; i<argnames.size(); ++i) argstr += "," + argnames[i];
     161           2 :       readInputLine("DUMPPDB DESCRIPTION=PATH STRIDE=" + wstride + " FMT=" + ofmt + " FILE=" + wfilename + " ARG=" + reference_data );
     162           2 :     } else readInputLine("DUMPPDB DESCRIPTION=PATH STRIDE=" + wstride + " FMT=" + ofmt + " FILE=" + wfilename + " ATOMS=" + reference_data + " ATOM_INDICES=" + atomstr );
     163             :   }
     164           4 :   log<<"  Bibliography "<<plumed.cite("Diaz Leines and Ensing, Phys. Rev. Lett. 109, 020601 (2012)")<<"\n";
     165           4 : }
     166             : 
     167             : }
     168             : }