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 "tools/FileTools.h"
27 : #include "Path.h"
28 :
29 : //+PLUMEDOC COLVAR ADAPTIVE_PATH
30 : /*
31 : Compute path collective variables that adapt to the lowest free energy path connecting states A and B.
32 :
33 : This shortcut can be used to implement the method discussed in the second paper cited below. This method allows you to
34 : run a simulation with [PATH](PATH.md) collective variables that forces the system to transition between state A and state B.
35 : If you run with ADAPTIVE_PATH rather than [PATH](PATH.md) the CV adapts as the simulation progresses and should find the
36 : lowest free energy path that connects the two states.
37 :
38 : The Path Collective Variables developed by Branduardi that are described in the first paper cited below allow one
39 : to compute the progress along a high-dimensional path and the distance from the high-dimensional
40 : path. Instead of computing progress along the path ($s$) using Branduardi's method we use the expressions that are used in
41 : [GEOMETRIC_PATH](GEOMETRIC_PATH.md) and that are introduced in the second paper that is cited below. The progress along the path is thus computed using:
42 :
43 : $$
44 : 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}
45 : $$
46 :
47 : In this expression $\mathbf{v}_1$ and $\mathbf{v}_3$ are the vectors connecting the current position to the closest and second closest node of the path,
48 : respectfully and $i_1$ and $i_2$ are the projections of the closest and second closest frames of the path. $\mathbf{v}_2$, meanwhile, is the
49 : vector connecting the closest frame to the second closest frame. The distance from the path, $z$ is calculated using:
50 :
51 : $$
52 : 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 }
53 : $$
54 :
55 : Notice that these are the definitions of $s$ and $z$ from the [GEOMETRIC_PATH](GEOMETRIC_PATH.md). The reason for this is that
56 : the adaptive path method implemented in this action was inspired by the work of Diaz and Ensing that was introduced in the second paper cited below.
57 : To learn more about how the path is adapted we strongly recommend reading this paper.
58 :
59 : The example input below shows how to use the adaptive path.
60 :
61 : ```plumed
62 : #SETTINGS INPUTFILES=regtest/mapping/rt-adapt/mypath.pdb
63 : d1: DISTANCE ATOMS=1,2 COMPONENTS
64 : pp: ADAPTIVE_PATH ...
65 : ARG=d1.x,d1.y UPDATE=50 FIXED=5,15
66 : WFILE=out-path.pdb WSTRIDE=50
67 : REFERENCE=regtest/mapping/rt-adapt/mypath.pdb
68 : PROPERTY=path
69 : ...
70 : PRINT ARG=d1.x,d1.y,pp.* FILE=colvar
71 : ```
72 :
73 : The curved path here is defined using a series of points in a two dimensional space. Furthermore, the coordinates that these points should be projected at
74 : on the low dimensional path (the $i$ values in the equation for $s$ above) are specified in the `path` property of the input pdb file.
75 :
76 : If you expand the shortcuts in the input above you can see how
77 : we use a [GEOMETRIC_PATH](GEOMETRIC_PATH.md) shortcut to calculate our position along the path and distance from the path. In addition, we also collect
78 : information on the average displacement from the path using an [AVERAGE_PATH_DISPLACEMENT](AVERAGE_PATH_DISPLACEMENT.md) so that we can update the path every step.
79 : During these update steps the reference points on the path are shifted by the values that are stored in the [AVERAGE_PATH_DISPLACEMENT](AVERAGE_PATH_DISPLACEMENT.md)
80 : action. We then use the [REPARAMETERIZE_PATH](REPARAMETERIZE_PATH.md) to ensure that the new set of reference points on the path are all equally spaced. After the adaption
81 : steps the new reference points are used when calculating are position on and distance from the reference path.
82 :
83 : The input below illustrates how to the positions of atoms to define the adaptive path:
84 :
85 : ```plumed
86 : ##SETTINGS INPUTFILES=regtest/trajectories/path_msd/all.pdb
87 : p1b: ADAPTIVE_PATH REFERENCE=regtest/trajectories/path_msd/all.pdb FIXED=1,42 UPDATE=1000 WFILE=out-path.pdb WSTRIDE=20000
88 : PRINT ARG=p1b.* FILE=colvar_b STRIDE=1
89 : ```
90 :
91 : As we have not used the PROPERTY keyword in the input the $i$ values in the equation for $s$ above are set equal to 1, 2, 3...
92 : Furthermore, when an input like the one above is used the vectors $\mathbf{v}_1$, $\mathbf{v}_2$ and $\mathbf{v}_3$ from the expression above are computed using an [RMSD](RMSD.md)
93 : action with the DISPLACEMENT flag enabled. The instantaneous structure is thus aligned with the reference structures so as to remove motions due to translation of the center
94 : of mass and rotation of the reference frame. These vector of atomic displacements that do not include the translation and rotation are also used in the
95 : [AVERAGE_PATH_DISPLACEMENT](AVERAGE_PATH_DISPLACEMENT.md) action that is used to update the path.
96 :
97 : As we have not used the PROPERTY keyword in the input the $i$ values in the equation for $s$ above are set equal to 1, 2, 3...
98 :
99 :
100 : */
101 : //+ENDPLUMEDOC
102 :
103 : namespace PLMD {
104 : namespace mapping {
105 :
106 : class AdaptivePath : public ActionShortcut {
107 : public:
108 : static void registerKeywords( Keywords& keys );
109 : explicit AdaptivePath(const ActionOptions&);
110 : };
111 :
112 : PLUMED_REGISTER_ACTION(AdaptivePath,"ADAPTIVE_PATH")
113 :
114 10 : void AdaptivePath::registerKeywords( Keywords& keys ) {
115 10 : ActionShortcut::registerKeywords( keys );
116 10 : Path::registerInputFileKeywords( keys );
117 10 : keys.add("optional","PROPERTY","read in path coordinates by finding option with this label in remark of pdb frames");
118 10 : 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");
119 10 : 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");
120 10 : keys.add("compulsory","UPDATE","the frequency with which the path should be updated");
121 10 : keys.add("compulsory","TOLERANCE","1E-4","the tolerance to use for the path updating algorithm that makes all frames equidistant");
122 10 : keys.add("optional","WFILE","file on which to write out the path");
123 10 : keys.add("compulsory","FMT","%f","the format to use for output files");
124 10 : keys.add("compulsory","WSTRIDE","0,","frequency with which to write out the path");
125 20 : keys.setValueDescription("scalar","the position along and from the adaptive path");
126 10 : keys.addDOI("10.1063/1.2432340");
127 10 : keys.addDOI("10.1103/PhysRevLett.109.020601");
128 10 : keys.needsAction("GEOMETRIC_PATH");
129 10 : keys.needsAction("AVERAGE_PATH_DISPLACEMENT");
130 10 : keys.needsAction("REPARAMETERIZE_PATH");
131 10 : keys.needsAction("DUMPPDB");
132 10 : keys.needsAction("PDB2CONSTANT");
133 10 : keys.needsAction("DISPLACEMENT");
134 10 : keys.needsAction("CONSTANT");
135 10 : }
136 :
137 4 : AdaptivePath::AdaptivePath(const ActionOptions& ao):
138 : Action(ao),
139 4 : ActionShortcut(ao) {
140 : // Read in the arguments
141 4 : std::vector<std::string> argnames;
142 8 : parseVector("ARG",argnames);
143 : std::string reference_data, metric, mtype;
144 8 : parse("TYPE", mtype);
145 : std::string reference;
146 8 : parse("REFERENCE",reference);
147 : {
148 : //extra scope for fp and mypdb
149 4 : unique_FILE fp{std::fopen(reference.c_str(),"r")};
150 4 : PDB mypdb;
151 4 : if(!fp.get()) {
152 0 : error("could not open reference file " + reference );
153 : }
154 4 : bool do_read=mypdb.readFromFilepointer(fp.get(),false,0.1);
155 4 : if( !do_read ) {
156 0 : error("missing file " + reference );
157 : }
158 : // Create list of reference configurations that PLUMED will use
159 4 : Path::readInputFrames( reference, mtype, argnames, true, this, reference_data );
160 : // Now get coordinates on spath
161 4 : std::vector<std::string> pnames;
162 4 : parseVector("PROPERTY",pnames);
163 4 : Path::readPropertyInformation( pnames, getShortcutLabel(), reference, this );
164 : // Create action that computes the geometric path variablesa
165 4 : std::string propstr = getShortcutLabel() + "_ind";
166 4 : if( pnames.size()>0 ) {
167 2 : propstr = pnames[0] + "_ref";
168 : }
169 4 : if( argnames.size()>0 ) {
170 5 : readInputLine( getShortcutLabel() + ": GEOMETRIC_PATH ARG=" + getShortcutLabel() + "_data " + " PROPERTY=" + propstr + " REFERENCE=" + reference_data + " METRIC={DIFFERENCE}");
171 : } else {
172 : std::string num, align_str, displace_str;
173 3 : Tools::convert( mypdb.getOccupancy()[0], align_str );
174 3 : Tools::convert( mypdb.getBeta()[0], displace_str );
175 39 : for(unsigned j=1; j<mypdb.getAtomNumbers().size(); ++j ) {
176 36 : Tools::convert( mypdb.getOccupancy()[j], num );
177 36 : align_str += "," + num;
178 36 : Tools::convert( mypdb.getBeta()[j], num );
179 36 : displace_str += "," + num;
180 : }
181 6 : metric = "RMSD DISPLACEMENT TYPE=" + mtype + " ALIGN=" + align_str + " DISPLACE=" + displace_str;
182 18 : readInputLine( getShortcutLabel() + ": GEOMETRIC_PATH ARG=" + getShortcutLabel() + "_data.disp " + " PROPERTY=" + propstr + " REFERENCE=" + reference_data + " METRIC={" + metric + "} METRIC_COMPONENT=disp");
183 : }
184 4 : }
185 : // Create the object to accumulate the average path displacements
186 : std::string update, halflife;
187 4 : parse("HALFLIFE",halflife);
188 4 : parse("UPDATE",update);
189 8 : std::string refframes = " REFERENCE=" + getShortcutLabel() + "_pos";
190 4 : if( argnames.size()>0 ) {
191 5 : readInputLine( getShortcutLabel() + "_disp: AVERAGE_PATH_DISPLACEMENT ARG=" + getShortcutLabel() + "_data HALFLIFE=" + halflife + " CLEAR=" + update + " METRIC={DIFFERENCE} REFERENCE=" + reference_data );
192 : } else {
193 18 : readInputLine( getShortcutLabel() + "_disp: AVERAGE_PATH_DISPLACEMENT ARG=" + getShortcutLabel() + "_data.disp HALFLIFE=" + halflife + " CLEAR=" + update + " METRIC={" + metric + "} METRIC_COMPONENT=disp REFERENCE=" + reference_data );
194 : }
195 :
196 : // Create the object to update the path
197 : std::string fixedn, tolerance;
198 4 : parse("FIXED",fixedn);
199 4 : parse("TOLERANCE",tolerance);
200 4 : std::string component="METRIC_COMPONENT=disp";
201 4 : if( argnames.size()>0 ) {
202 : metric="DIFFERENCE";
203 : component="";
204 : }
205 32 : readInputLine("REPARAMETERIZE_PATH DISPLACE_FRAMES=" + getShortcutLabel() + "_disp FIXED=" + fixedn + " STRIDE=" + update + " TOL=" + tolerance + " METRIC={" + metric + "} " + component + " REFERENCE=" + reference_data );
206 :
207 : // Information for write out
208 : std::string wfilename;
209 8 : parse("WFILE",wfilename);
210 4 : if( wfilename.length()>0 ) {
211 : // This just gets the atom numbers for output
212 : std::string atomstr;
213 3 : if( argnames.size()==0 ) {
214 2 : unique_FILE fp{std::fopen(reference.c_str(),"r")};
215 : double fake_unit=0.1;
216 2 : PDB mypdb;
217 2 : bool do_read=mypdb.readFromFilepointer(fp.get(),false,fake_unit);
218 2 : if( !do_read ) {
219 0 : error("missing file " + reference );
220 : }
221 : std::string num;
222 2 : Tools::convert( mypdb.getAtomNumbers()[0].serial(), atomstr );
223 26 : for(unsigned j=1; j<mypdb.getAtomNumbers().size(); ++j ) {
224 24 : Tools::convert( mypdb.getAtomNumbers()[j].serial(), num );
225 24 : atomstr += "," + num;
226 : }
227 2 : Tools::convert( mypdb.getResidueNumber( mypdb.getAtomNumbers()[0] ), num );
228 2 : atomstr += " RESIDUE_INDICES=" + num;
229 26 : for(unsigned i=1; i<mypdb.getAtomNumbers().size(); ++i ) {
230 24 : Tools::convert( mypdb.getResidueNumber( mypdb.getAtomNumbers()[i] ), num );
231 24 : atomstr += "," + num;
232 : }
233 : std::string anum, dnum;
234 2 : Tools::convert( mypdb.getOccupancy()[0], anum );
235 2 : Tools::convert( mypdb.getBeta()[0], dnum );
236 2 : atomstr += " OCCUPANCY=" + anum;
237 26 : for(unsigned i=1; i<mypdb.getOccupancy().size(); ++i) {
238 24 : Tools::convert( mypdb.getOccupancy()[i], anum );
239 24 : atomstr += "," + anum;
240 : }
241 2 : atomstr += " BETA=" + dnum;
242 26 : for(unsigned i=1; i<mypdb.getBeta().size(); ++i) {
243 24 : Tools::convert( mypdb.getBeta()[i], dnum );
244 24 : atomstr += "," + dnum;
245 : }
246 2 : }
247 :
248 3 : if( wfilename.find(".pdb")==std::string::npos ) {
249 0 : error("output must be to a pdb file");
250 : }
251 : std::string ofmt, pframes, wstride;
252 3 : parse("WSTRIDE",wstride);
253 6 : parse("FMT",ofmt);
254 3 : if( argnames.size()>0 ) {
255 1 : std::string argstr = argnames[0];
256 2 : for(unsigned i=1; i<argnames.size(); ++i) {
257 2 : argstr += "," + argnames[i];
258 : }
259 5 : readInputLine("DUMPPDB DESCRIPTION=PATH STRIDE=" + wstride + " FMT=" + ofmt + " FILE=" + wfilename + " ARG=" + reference_data );
260 : } else {
261 12 : readInputLine("DUMPPDB DESCRIPTION=PATH STRIDE=" + wstride + " FMT=" + ofmt + " FILE=" + wfilename + " ATOMS=" + reference_data + " ATOM_INDICES=" + atomstr );
262 : }
263 : }
264 8 : log<<" Bibliography "<<plumed.cite("Diaz Leines and Ensing, Phys. Rev. Lett. 109, 020601 (2012)")<<"\n";
265 181 : }
266 :
267 : }
268 : }
|
