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 );
99 6 : Path::registerInputFileKeywords( keys );
100 6 : keys.add("optional","PROPERTY","read in path coordinates by finding option with this label in remark of pdb frames");
101 6 : 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");
102 6 : 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");
103 6 : keys.add("compulsory","UPDATE","the frequency with which the path should be updated");
104 6 : keys.add("compulsory","TOLERANCE","1E-6","the tolerance to use for the path updating algorithm that makes all frames equidistant");
105 6 : keys.add("optional","WFILE","file on which to write out the path");
106 6 : keys.add("compulsory","FMT","%f","the format to use for output files");
107 6 : keys.add("compulsory","WSTRIDE","0,","frequency with which to write out the path");
108 12 : keys.setValueDescription("scalar","the position along and from the adaptive path");
109 6 : keys.needsAction("GEOMETRIC_PATH");
110 6 : keys.needsAction("AVERAGE_PATH_DISPLACEMENT");
111 6 : keys.needsAction("REPARAMETERIZE_PATH");
112 6 : keys.needsAction("DUMPPDB");
113 6 : keys.needsAction("PDB2CONSTANT");
114 6 : keys.needsAction("DISPLACEMENT");
115 6 : keys.needsAction("CONSTANT");
116 6 : }
117 :
118 2 : AdaptivePath::AdaptivePath(const ActionOptions& ao):
119 : Action(ao),
120 2 : ActionShortcut(ao) {
121 : // Read in the arguments
122 : std::vector<std::string> argnames;
123 4 : parseVector("ARG",argnames);
124 : std::string reference_data, metric, mtype;
125 4 : parse("TYPE", mtype);
126 : std::string reference;
127 4 : parse("REFERENCE",reference);
128 2 : FILE* fp=std::fopen(reference.c_str(),"r");
129 2 : PDB mypdb;
130 2 : if(!fp) {
131 0 : error("could not open reference file " + reference );
132 : }
133 2 : bool do_read=mypdb.readFromFilepointer(fp,false,0.1);
134 2 : if( !do_read ) {
135 0 : error("missing file " + reference );
136 : }
137 : // Create list of reference configurations that PLUMED will use
138 2 : Path::readInputFrames( reference, mtype, argnames, true, this, reference_data );
139 : // Now get coordinates on spath
140 : std::vector<std::string> pnames;
141 2 : parseVector("PROPERTY",pnames);
142 2 : Path::readPropertyInformation( pnames, getShortcutLabel(), reference, this );
143 : // Create action that computes the geometric path variablesa
144 2 : std::string propstr = getShortcutLabel() + "_ind";
145 2 : if( pnames.size()>0 ) {
146 2 : propstr = pnames[0] + "_ref";
147 : }
148 2 : if( argnames.size()>0 ) {
149 2 : readInputLine( getShortcutLabel() + ": GEOMETRIC_PATH ARG=" + getShortcutLabel() + "_data " + " PROPERTY=" + propstr + " REFERENCE=" + reference_data + " METRIC={DIFFERENCE}");
150 : } else {
151 : std::string num, align_str, displace_str;
152 1 : Tools::convert( mypdb.getOccupancy()[0], align_str );
153 1 : Tools::convert( mypdb.getBeta()[0], displace_str );
154 13 : for(unsigned j=1; j<mypdb.getAtomNumbers().size(); ++j ) {
155 12 : Tools::convert( mypdb.getOccupancy()[j], num );
156 12 : align_str += "," + num;
157 12 : Tools::convert( mypdb.getBeta()[0], num );
158 24 : displace_str += "," + num;
159 : }
160 2 : metric = "RMSD_VECTOR DISPLACEMENT TYPE=" + mtype + " ALIGN=" + align_str + " DISPLACE=" + displace_str;
161 2 : readInputLine( getShortcutLabel() + ": GEOMETRIC_PATH ARG=" + getShortcutLabel() + "_data.disp " + " PROPERTY=" + propstr + " REFERENCE=" + reference_data + " METRIC={" + metric + "} METRIC_COMPONENT=disp");
162 : }
163 : // Create the object to accumulate the average path displacements
164 : std::string update, halflife;
165 2 : parse("HALFLIFE",halflife);
166 2 : parse("UPDATE",update);
167 4 : std::string refframes = " REFERENCE=" + getShortcutLabel() + "_pos";
168 2 : if( argnames.size()>0 ) {
169 2 : readInputLine( getShortcutLabel() + "_disp: AVERAGE_PATH_DISPLACEMENT ARG=" + getShortcutLabel() + "_data HALFLIFE=" + halflife + " CLEAR=" + update + " METRIC={DIFFERENCE} REFERENCE=" + reference_data );
170 : } else {
171 2 : readInputLine( getShortcutLabel() + "_disp: AVERAGE_PATH_DISPLACEMENT ARG=" + getShortcutLabel() + "_data.disp HALFLIFE=" + halflife + " CLEAR=" + update + " METRIC={" + metric + "} METRIC_COMPONENT=disp REFERENCE=" + reference_data );
172 : }
173 :
174 : // Create the object to update the path
175 : std::string fixedn;
176 2 : parse("FIXED",fixedn);
177 2 : std::string component="METRIC_COMPONENT=disp";
178 2 : if( argnames.size()>0 ) {
179 : metric="DIFFERENCE";
180 : component="";
181 : }
182 4 : readInputLine("REPARAMETERIZE_PATH DISPLACE_FRAMES=" + getShortcutLabel() + "_disp FIXED=" + fixedn + " STRIDE=" + update + " METRIC={" + metric + "} " + component + " REFERENCE=" + reference_data );
183 :
184 : // Information for write out
185 : std::string wfilename;
186 4 : parse("WFILE",wfilename);
187 2 : if( wfilename.length()>0 ) {
188 : // This just gets the atom numbers for output
189 : std::string atomstr;
190 2 : if( argnames.size()==0 ) {
191 1 : FILE* fp=std::fopen(reference.c_str(),"r");
192 : double fake_unit=0.1;
193 1 : PDB mypdb;
194 1 : bool do_read=mypdb.readFromFilepointer(fp,false,fake_unit);
195 : std::string num;
196 1 : Tools::convert( mypdb.getAtomNumbers()[0].serial(), atomstr );
197 13 : for(unsigned j=1; j<mypdb.getAtomNumbers().size(); ++j ) {
198 12 : Tools::convert( mypdb.getAtomNumbers()[j].serial(), num );
199 24 : atomstr += "," + num;
200 : }
201 1 : Tools::convert( mypdb.getResidueNumber( mypdb.getAtomNumbers()[0] ), num );
202 1 : atomstr += " RESIDUE_INDICES=" + num;
203 13 : for(unsigned i=1; i<mypdb.getAtomNumbers().size(); ++i ) {
204 12 : Tools::convert( mypdb.getResidueNumber( mypdb.getAtomNumbers()[i] ), num );
205 24 : atomstr += "," + num;
206 : }
207 : std::string anum, dnum;
208 1 : Tools::convert( mypdb.getOccupancy()[0], anum );
209 1 : Tools::convert( mypdb.getBeta()[0], dnum );
210 1 : atomstr += " OCCUPANCY=" + anum;
211 13 : for(unsigned i=1; i<mypdb.getOccupancy().size(); ++i) {
212 12 : Tools::convert( mypdb.getOccupancy()[i], anum );
213 24 : atomstr += "," + anum;
214 : }
215 1 : atomstr += " BETA=" + dnum;
216 13 : for(unsigned i=1; i<mypdb.getBeta().size(); ++i) {
217 12 : Tools::convert( mypdb.getBeta()[i], dnum );
218 24 : atomstr += "," + dnum;
219 : }
220 1 : }
221 :
222 2 : if( wfilename.find(".pdb")==std::string::npos ) {
223 0 : error("output must be to a pdb file");
224 : }
225 : std::string ofmt, pframes, wstride;
226 2 : parse("WSTRIDE",wstride);
227 4 : parse("FMT",ofmt);
228 2 : if( argnames.size()>0 ) {
229 1 : std::string argstr = argnames[0];
230 2 : for(unsigned i=1; i<argnames.size(); ++i) {
231 2 : argstr += "," + argnames[i];
232 : }
233 2 : readInputLine("DUMPPDB DESCRIPTION=PATH STRIDE=" + wstride + " FMT=" + ofmt + " FILE=" + wfilename + " ARG=" + reference_data );
234 : } else {
235 2 : readInputLine("DUMPPDB DESCRIPTION=PATH STRIDE=" + wstride + " FMT=" + ofmt + " FILE=" + wfilename + " ATOMS=" + reference_data + " ATOM_INDICES=" + atomstr );
236 : }
237 : }
238 4 : log<<" Bibliography "<<plumed.cite("Diaz Leines and Ensing, Phys. Rev. Lett. 109, 020601 (2012)")<<"\n";
239 4 : }
240 :
241 : }
242 : }
|
