Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2018 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/ActionWithValue.h"
23 : #include "core/ActionPilot.h"
24 : #include "core/ActionWithArguments.h"
25 : #include "core/ActionRegister.h"
26 : #include "tools/Matrix.h"
27 : #include "PathProjectionCalculator.h"
28 :
29 : //+PLUMEDOC ANALYSIS AVERAGE_PATH_DISPLACEMENT
30 : /*
31 : Accumulate the distances between the reference frames in the paths and the configurations visited
32 :
33 : This action is used in the shortcut for [ADAPTIVE_PATH](ADAPTIVE_PATH.md) as it is helps us to adapt and refine
34 : the defintition of a path CV based on the sampling that is observed in the trajectory.
35 :
36 : By way of reminder path CVs were introduced by Branduardi _et al._ in the first paper cited in the bibliography below.
37 : In PLUMED this method is implemented in [PATH](PATH.md) and works by definition a position along the path as:
38 :
39 : $$
40 : s = \frac{ \sum_{i=1}^N i \exp( -\lambda R[X - X_i] ) }{ \sum_{i=1}^N \exp( -\lambda R[X - X_i] ) }
41 : $$
42 :
43 : while the distance from the path (z) is measured using:
44 :
45 : $$
46 : z = -\frac{1}{\lambda} \ln\left[ \sum_{i=1}^N \exp( -\lambda R[X - X_i] ) \right]
47 : $$
48 :
49 : In these expressions $N$ high-dimensional frames ($X_i$) are used to describe the path in the high-dimensional
50 : space. The two expressions above are then functions of the distances from each of the high-dimensional frames $R[X - X_i]$.
51 :
52 : When a simulator uses this method of path CVs the assumption is that when the system undergoes the transition from state A
53 : to state B it passes through a narrow (non-linear) tube that is centered on the $s$ coordinate that is defiend above in terms the
54 : various milestones that are used to define the path. In other words, when the system transitions from A to be B it does not pass along the $s$
55 : coordinate. If, however, we average the vector connecting each position the system passes through to the nearest
56 : point on the path that underpins the definition of the $s$ we should get zero as the $s$ coordinate is centered on the path connecting the two states.
57 :
58 : With that theory in mind this action can be used to collect the average displacement between the points that trajectories are passing through and the path.
59 : This action can thus be used to update the definitions of the milestones that are used in the definition of a [PATH](PATH.md) or [GPATH](GPATH.md) coordinate
60 : so that the PATH more clearly passes through the center of the path that the system has traversed along in passing from state A to state B. These path displacements
61 : are accumulated as follows. You first determine the vectors, $\mathbf{v}_1$ and $\mathbf{v}_3$ that connect the instaneous position to the closest and second closest milestones on the
62 : current path. You can then compute the following:
63 :
64 : $$
65 : \delta = \frac{ \sqrt{( \mathbf{v}_1\cdot\mathbf{v}_2 )^2 - |\mathbf{v}_2|^2(|\mathbf{v}_1|^2 - |\mathbf{v}_3|^2) } + \mathbf{v}_1\cdot\mathbf{v}_2 }{2|\mathbf{v}_2|^2} - \frac{1}{2}
66 : $$
67 :
68 : where $\mathbf{v}_2$ is the vector connecting the closest and second closest milestone on your path. Displacements for these two milestones are then computed as:
69 :
70 : $$
71 : d_1 = (1 + \delta) \left(\mathbf{v}_1 - \delta \mathbf{v}_2\right) \qquad d_2 = -\delta \left(\mathbf{v}_1 - \delta \mathbf{v}_2\right)
72 : $$
73 :
74 : These displacement vectors for the two closest nodes are computed on every accumulate step and are then averaged with the factors of $(1 + \delta)$ and $-\delta$ in these expressions
75 : serving as weights that are used when normalising these averages.
76 :
77 : ## Examples
78 :
79 : Suppose that you are interested in a transition between two states A and B and that you have representative structures for these two states in the files `start.pdb` and `end.pdb`. You can
80 : use [pathtools](pathtools.md) to create an initial (linear) path connecting these two states as follows:
81 :
82 : ```plumed
83 : plumed pathtools --start start.pdb --end end.pdb --nframes 20 --metric OPTIMAL --out path.pdb
84 : ```
85 :
86 : This path is likely not even close to the transition pathway that the system takes between these two states as the assumption that the path is linear is pretty severe. We can nevertheless use
87 : a path generated by such a command in a [MOVINGRESTRAINT](MOVINGRESTRAINT.md) command like the one shown below:
88 :
89 : ```plumed
90 : #SETTINGS INPUTFILES=regtest/trajectories/path_msd/all.pdb
91 : p1: PATH REFERENCE=regtest/trajectories/path_msd/all.pdb TYPE=OPTIMAL LAMBDA=69087
92 : mv: MOVINGRESTRAINT ...
93 : ARG=p1.s KAPPA0=1000
94 : STEP0=0 AT0=1
95 : STEP1=1000000 AT1=42
96 : STEP2=2000000 AT2=1
97 : STEP3=3000000 AT3=42
98 : STEP4=4000000 AT4=1
99 : STEP5=5000000 AT5=42
100 : STEP6=6000000 AT6=1
101 : STEP7=7000000 AT7=42
102 : STEP8=8000000 AT8=1
103 : ...
104 : PRINT ARG=p1.s,p1.z FILE=colvar
105 : ```
106 :
107 : This input hopefully (you need to check by looking at your trajectory) drives our system between the two states of interest 8 times - four of these transitions are from A to B and four are from B to A.
108 : You can then analyze the trajectories that are generated using [driver](driver.md) and reparameterize a new path that passes more closely through the sampled trajectories using:
109 :
110 : ```plumed
111 : #SETTINGS INPUTFILES=regtest/trajectories/path_msd/all.pdb
112 : rmsd: RMSD REFERENCE=regtest/trajectories/path_msd/all.pdb DISPLACEMENT TYPE=OPTIMAL
113 : # Accumulate the average displacement between the reference path and the trajectories that have sampled the transition
114 : disp: AVERAGE_PATH_DISPLACEMENT ARG=rmsd.disp STRIDE=1 METRIC={RMSD DISPLACEMENT TYPE=OPTIMAL ALIGN=1,1,1,1,1,1,1,1,1,1,1,1,1 DISPLACE=1,1,1,1,1,1,1,1,1,1,1,1,1} METRIC_COMPONENT=disp REFERENCE=rmsd_ref
115 : # Now displace the original path by the accumulated displacement and reparameterize so that all frames are equally spaced
116 : REPARAMETERIZE_PATH DISPLACE_FRAMES=disp FIXED=1,42 METRIC={RMSD DISPLACEMENT TYPE=OPTIMAL ALIGN=1,1,1,1,1,1,1,1,1,1,1,1,1 DISPLACE=1,1,1,1,1,1,1,1,1,1,1,1,1} METRIC_COMPONENT=disp REFERENCE=rmsd_ref
117 : # And output the final reparameterized path at the end of the simulation
118 : DUMPPDB DESCRIPTION=PATH STRIDE=0 FILE=outpatb.pdb ATOMS=rmsd_ref ATOM_INDICES=1-13
119 : ```
120 :
121 : You can then try the same calculation with the reparameterized path before eventually using a method such as [METAD](METAD.md) to obtain the free energy as a function of your $s$ path coordinate.
122 :
123 : Notice that the METRIC keyword appears here as the calculation involves calculating the vectors connecting the milestones on your path. This keyword operates in the same way that is described in
124 : the documentation for [GEOMETRIC_PATH](GEOMETRIC_PATH.md).
125 :
126 : */
127 : //+ENDPLUMEDOC
128 :
129 : namespace PLMD {
130 : namespace mapping {
131 :
132 : class PathDisplacements : public ActionWithValue, public ActionPilot, public ActionWithArguments {
133 : private:
134 : bool clearnextstep;
135 : unsigned clearstride;
136 : double fadefact;
137 : std::vector<double> wsum, displace_v;
138 : Matrix<double> displacements;
139 : PathProjectionCalculator path_projector;
140 : public:
141 : static void registerKeywords( Keywords& keys );
142 : explicit PathDisplacements(const ActionOptions&);
143 : unsigned getNumberOfDerivatives();
144 573 : void clearDerivatives( const bool& force=false ) {}
145 573 : void calculate() {}
146 573 : void apply() {}
147 : void update();
148 : };
149 :
150 : PLUMED_REGISTER_ACTION(PathDisplacements,"AVERAGE_PATH_DISPLACEMENT")
151 :
152 6 : void PathDisplacements::registerKeywords( Keywords& keys ) {
153 6 : Action::registerKeywords( keys );
154 6 : ActionWithValue::registerKeywords( keys );
155 6 : ActionPilot::registerKeywords( keys );
156 6 : ActionWithArguments::registerKeywords( keys );
157 6 : PathProjectionCalculator::registerKeywords( keys );
158 12 : keys.remove("NUMERICAL_DERIVATIVES");
159 6 : keys.add("compulsory","STRIDE","1","the frequency with which the average displacements should be collected and added to the average displacements");
160 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");
161 6 : keys.add("compulsory","CLEAR","0","the frequency with which to clear all the accumulated data. The default value "
162 : "of 0 implies that all the data will be used and that the grid will never be cleared");
163 12 : keys.setValueDescription("matrix","matrix containing the average displacement between the trajectory and each of the landmarks that makes up the path");
164 6 : keys.addDOI("10.1063/1.2432340");
165 6 : keys.addDOI("10.1063/1.3660208");
166 6 : }
167 :
168 2 : PathDisplacements::PathDisplacements(const ActionOptions& ao):
169 : Action(ao),
170 : ActionWithValue(ao),
171 : ActionPilot(ao),
172 : ActionWithArguments(ao),
173 2 : clearnextstep(false),
174 2 : path_projector(this) {
175 : // Read in clear instructions
176 2 : parse("CLEAR",clearstride);
177 2 : if( clearstride>0 ) {
178 2 : if( clearstride%getStride()!=0 ) {
179 0 : error("CLEAR parameter must be a multiple of STRIDE");
180 : }
181 2 : log.printf(" clearing average every %u steps \n",clearstride);
182 : }
183 : double halflife;
184 2 : parse("HALFLIFE",halflife);
185 2 : log.printf(" weight of contribution to frame halves every %f steps \n",halflife);
186 2 : if( halflife<0 ) {
187 2 : fadefact=1.0;
188 : } else {
189 0 : fadefact = exp( -0.693147180559945 / static_cast<double>(halflife) );
190 : }
191 : // Now create the weights vector and displacements matrix
192 2 : unsigned nrows = getPntrToArgument(0)->getShape()[0];
193 2 : unsigned ncols = getPntrToArgument(0)->getShape()[1];
194 2 : wsum.resize( nrows );
195 2 : displacements.resize( nrows, ncols );
196 64 : for(unsigned i=0; i<nrows; ++i) {
197 62 : wsum[i]=0;
198 1740 : for(unsigned j=0; j<ncols; ++j) {
199 1678 : displacements(i,j)=0;
200 : }
201 : }
202 : // Add bibliography
203 4 : log<<" Bibliography "<<plumed.cite("Diaz Leines and Ensing, Phys. Rev. Lett. 109, 020601 (2012)")<<"\n";
204 : // And create a value to hold the displacements
205 2 : std::vector<std::size_t> shape(2);
206 2 : shape[0]=nrows;
207 2 : shape[1]=ncols;
208 2 : addValue( shape );
209 2 : setNotPeriodic();
210 2 : getPntrToComponent(0)->reshapeMatrixStore( shape[1] );
211 2 : }
212 :
213 0 : unsigned PathDisplacements::getNumberOfDerivatives() {
214 0 : return 0;
215 : }
216 :
217 573 : void PathDisplacements::update() {
218 573 : unsigned nrows = getPntrToArgument(0)->getShape()[0];
219 573 : unsigned ncols = getPntrToArgument(0)->getShape()[1];
220 :
221 573 : if( clearnextstep ) {
222 : unsigned k=0;
223 1010 : for(unsigned i=0; i<nrows; ++i) {
224 38700 : for(unsigned j=0; j<ncols; ++j) {
225 37714 : displacements(i,j)=0;
226 37714 : getPntrToComponent(0)->set(k,0);
227 37714 : k++;
228 : }
229 : }
230 24 : clearnextstep=false;
231 : }
232 :
233 : unsigned k=0, iclose1=0, iclose2=0;
234 : double v1v1=0, v3v3=0;
235 22417 : for(unsigned i=0; i<nrows; ++i) {
236 : double dist = 0;
237 799020 : for(unsigned j=0; j<ncols; ++j) {
238 777176 : double tmp = getPntrToArgument(0)->get(k);
239 777176 : dist += tmp*tmp;
240 777176 : k++;
241 : }
242 21844 : if( i==0 ) {
243 : v1v1 = dist;
244 : iclose1 = 0;
245 21271 : } else if( dist<v1v1 ) {
246 : v3v3=v1v1;
247 : v1v1=dist;
248 : iclose2=iclose1;
249 : iclose1=i;
250 10991 : } else if( i==1 ) {
251 : v3v3=dist;
252 : iclose2=1;
253 10991 : } else if( dist<v3v3 ) {
254 : v3v3=dist;
255 : iclose2=i;
256 : }
257 : }
258 : // And find third closest point
259 573 : int isign = iclose1 - iclose2;
260 : if( isign>1 ) {
261 : isign=1;
262 : } else if( isign<-1 ) {
263 : isign=-1;
264 : }
265 573 : int iclose3 = iclose1 + isign;
266 573 : unsigned ifrom=iclose1, ito=iclose3;
267 573 : if( iclose3<0 || static_cast<unsigned>(iclose3)>=nrows ) {
268 0 : ifrom=iclose2;
269 0 : ito=iclose1;
270 : }
271 :
272 : // Calculate the dot product of v1 with v2
273 573 : path_projector.getDisplaceVector( ifrom, ito, displace_v );
274 : double v2v2=0, v1v2=0;
275 573 : unsigned kclose1 = iclose1*ncols;
276 19183 : for(unsigned i=0; i<displace_v.size(); ++i) {
277 18610 : v2v2 += displace_v[i]*displace_v[i];
278 18610 : v1v2 += displace_v[i]*getPntrToArgument(0)->get(kclose1+i);
279 : }
280 :
281 573 : double root = sqrt( v1v2*v1v2 - v2v2 * ( v1v1 - v3v3) );
282 573 : double dx = 0.5 * ( (root + v1v2) / v2v2 - 1.);
283 573 : double weight2 = -1.* dx;
284 573 : double weight1 = 1.0 + dx;
285 573 : if( weight1>1.0 ) {
286 : weight1=1.0;
287 : weight2=0.0;
288 573 : } else if( weight2>1.0 ) {
289 : weight1=0.0;
290 : weight2=1.0;
291 : }
292 :
293 : // Accumulate displacements for path
294 19183 : for(unsigned i=0; i<ncols; ++i) {
295 18610 : double displace = getPntrToArgument(0)->get(kclose1+i) - dx*displace_v[i];
296 18610 : displacements(iclose1,i) += weight1 * displace;
297 18610 : displacements(iclose2,i) += weight2 * displace;
298 : }
299 :
300 : // Update weight accumulators
301 573 : wsum[iclose1] *= fadefact;
302 573 : wsum[iclose2] *= fadefact;
303 573 : wsum[iclose1] += weight1;
304 573 : wsum[iclose2] += weight2;
305 :
306 : // Update numbers in values
307 573 : if( wsum[iclose1] > epsilon ) {
308 19183 : for(unsigned i=0; i<ncols; ++i) {
309 18610 : getPntrToComponent(0)->set( kclose1+i, displacements(iclose1,i) / wsum[iclose1] );
310 : }
311 : }
312 573 : if( wsum[iclose2] > epsilon ) {
313 573 : unsigned kclose2 = iclose2*ncols;
314 19183 : for(unsigned i=0; i<ncols; ++i) {
315 18610 : getPntrToComponent(0)->set( kclose2+i, displacements(iclose2,i) / wsum[iclose2] );
316 : }
317 : }
318 :
319 : // Clear if required
320 573 : if( (getStep()>0 && clearstride>0 && getStep()%clearstride==0) ) {
321 25 : clearnextstep=true;
322 : }
323 573 : }
324 :
325 : }
326 : }
|
