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 defined 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 instantaneous 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 1665 : void clearDerivatives( const bool& force=false ) {}
145 1665 : void calculate() {}
146 1665 : void apply() {}
147 : void update();
148 : };
149 :
150 : PLUMED_REGISTER_ACTION(PathDisplacements,"AVERAGE_PATH_DISPLACEMENT")
151 :
152 12 : void PathDisplacements::registerKeywords( Keywords& keys ) {
153 12 : Action::registerKeywords( keys );
154 12 : ActionWithValue::registerKeywords( keys );
155 12 : ActionPilot::registerKeywords( keys );
156 12 : ActionWithArguments::registerKeywords( keys );
157 12 : PathProjectionCalculator::registerKeywords( keys );
158 24 : keys.remove("NUMERICAL_DERIVATIVES");
159 12 : keys.add("compulsory","STRIDE","1","the frequency with which the average displacements should be collected and added to the average displacements");
160 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");
161 12 : 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 24 : keys.setValueDescription("matrix","matrix containing the average displacement between the trajectory and each of the landmarks that makes up the path");
164 12 : keys.addDOI("10.1063/1.2432340");
165 12 : keys.addDOI("10.1063/1.3660208");
166 12 : }
167 :
168 4 : PathDisplacements::PathDisplacements(const ActionOptions& ao):
169 : Action(ao),
170 : ActionWithValue(ao),
171 : ActionPilot(ao),
172 : ActionWithArguments(ao),
173 4 : clearnextstep(false),
174 4 : path_projector(this) {
175 : // Read in clear instructions
176 4 : parse("CLEAR",clearstride);
177 4 : if( clearstride>0 ) {
178 4 : if( clearstride%getStride()!=0 ) {
179 0 : error("CLEAR parameter must be a multiple of STRIDE");
180 : }
181 4 : log.printf(" clearing average every %u steps \n",clearstride);
182 : }
183 : double halflife;
184 4 : parse("HALFLIFE",halflife);
185 4 : log.printf(" weight of contribution to frame halves every %f steps \n",halflife);
186 4 : if( halflife<0 ) {
187 3 : fadefact=1.0;
188 : } else {
189 1 : fadefact = exp( -0.693147180559945 / static_cast<double>(halflife) );
190 : }
191 : // Now create the weights vector and displacements matrix
192 4 : unsigned nrows = getPntrToArgument(0)->getShape()[0];
193 4 : unsigned ncols = getPntrToArgument(0)->getShape()[1];
194 4 : wsum.resize( nrows );
195 4 : displacements.resize( nrows, ncols );
196 84 : for(unsigned i=0; i<nrows; ++i) {
197 80 : wsum[i]=0;
198 2460 : for(unsigned j=0; j<ncols; ++j) {
199 2380 : 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 4 : std::vector<std::size_t> shape(2);
206 4 : shape[0]=nrows;
207 4 : shape[1]=ncols;
208 4 : addValue( shape );
209 4 : setNotPeriodic();
210 4 : getPntrToComponent(0)->reshapeMatrixStore( shape[1] );
211 4 : }
212 :
213 0 : unsigned PathDisplacements::getNumberOfDerivatives() {
214 0 : return 0;
215 : }
216 :
217 1665 : void PathDisplacements::update() {
218 1665 : unsigned nrows = getPntrToArgument(0)->getShape()[0];
219 1665 : unsigned ncols = getPntrToArgument(0)->getShape()[1];
220 :
221 1665 : if( clearnextstep ) {
222 : unsigned k=0;
223 1820 : for(unsigned i=0; i<nrows; ++i) {
224 67860 : for(unsigned j=0; j<ncols; ++j) {
225 66145 : displacements(i,j)=0;
226 66145 : getPntrToComponent(0)->set(k,0);
227 66145 : k++;
228 : }
229 : }
230 105 : clearnextstep=false;
231 : }
232 :
233 : unsigned k=0, iclose1=0, iclose2=0;
234 : double v1v1=0, v3v3=0;
235 33337 : for(unsigned i=0; i<nrows; ++i) {
236 : double dist = 0;
237 1192140 : for(unsigned j=0; j<ncols; ++j) {
238 1160468 : double tmp = getPntrToArgument(0)->get(k);
239 1160468 : dist += tmp*tmp;
240 1160468 : k++;
241 : }
242 31672 : if( i==0 ) {
243 : v1v1 = dist;
244 : iclose1 = 0;
245 30007 : } else if( dist<v1v1 ) {
246 : v3v3=v1v1;
247 : v1v1=dist;
248 : iclose2=iclose1;
249 : iclose1=i;
250 15285 : } else if( i==1 ) {
251 : v3v3=dist;
252 : iclose2=1;
253 15179 : } else if( dist<v3v3 ) {
254 : v3v3=dist;
255 : iclose2=i;
256 : }
257 : }
258 : // And find third closest point
259 1665 : int isign = iclose1 - iclose2;
260 : if( isign>1 ) {
261 : isign=1;
262 1665 : } else if( isign<-1 ) {
263 : isign=-1;
264 : }
265 1665 : int iclose3 = iclose1 + isign;
266 1665 : unsigned ifrom=iclose1, ito=iclose3;
267 1665 : if( iclose3<0 || static_cast<unsigned>(iclose3)>=nrows ) {
268 230 : ifrom=iclose2;
269 230 : ito=iclose1;
270 : }
271 :
272 : // Calculate the dot product of v1 with v2
273 1665 : path_projector.getDisplaceVector( ifrom, ito, displace_v );
274 : double v2v2=0, v1v2=0;
275 1665 : unsigned kclose1 = iclose1*ncols;
276 62863 : for(unsigned i=0; i<displace_v.size(); ++i) {
277 61198 : v2v2 += displace_v[i]*displace_v[i];
278 61198 : v1v2 += displace_v[i]*getPntrToArgument(0)->get(kclose1+i);
279 : }
280 :
281 1665 : double root = sqrt( v1v2*v1v2 - v2v2 * ( v1v1 - v3v3) );
282 1665 : double dx = 0.5 * ( (root + v1v2) / v2v2 - 1.);
283 1665 : double weight2 = -1.* dx;
284 1665 : double weight1 = 1.0 + dx;
285 1665 : if( weight1>1.0 ) {
286 : weight1=1.0;
287 : weight2=0.0;
288 1523 : } else if( weight2>1.0 ) {
289 : weight1=0.0;
290 : weight2=1.0;
291 : }
292 :
293 : // Accumulate displacements for path
294 62863 : for(unsigned i=0; i<ncols; ++i) {
295 61198 : double displace = getPntrToArgument(0)->get(kclose1+i) - dx*displace_v[i];
296 61198 : displacements(iclose1,i) += weight1 * displace;
297 61198 : displacements(iclose2,i) += weight2 * displace;
298 : }
299 :
300 : // Update weight accumulators
301 1665 : wsum[iclose1] *= fadefact;
302 1665 : wsum[iclose2] *= fadefact;
303 1665 : wsum[iclose1] += weight1;
304 1665 : wsum[iclose2] += weight2;
305 :
306 : // Update numbers in values
307 1665 : if( wsum[iclose1] > epsilon ) {
308 62863 : for(unsigned i=0; i<ncols; ++i) {
309 61198 : getPntrToComponent(0)->set( kclose1+i, displacements(iclose1,i) / wsum[iclose1] );
310 : }
311 : }
312 1665 : if( wsum[iclose2] > epsilon ) {
313 1665 : unsigned kclose2 = iclose2*ncols;
314 62863 : for(unsigned i=0; i<ncols; ++i) {
315 61198 : getPntrToComponent(0)->set( kclose2+i, displacements(iclose2,i) / wsum[iclose2] );
316 : }
317 : }
318 :
319 : // Clear if required
320 1665 : if( (getStep()>0 && clearstride>0 && getStep()%clearstride==0) ) {
321 106 : clearnextstep=true;
322 : }
323 1665 : }
324 :
325 : }
326 : }
|
