Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2023 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 :
23 : #include "Function.h"
24 : #include "core/ActionRegister.h"
25 :
26 : namespace PLMD {
27 : namespace function {
28 :
29 : //+PLUMEDOC FUNCTION FUNCPATHMSD
30 : /*
31 : This function calculates path collective variables.
32 :
33 : This is the Path Collective Variables implementation described in the paper from the bibliography.
34 : This variable computes the progress along a given set of frames that is provided
35 : in input ("s" component) and the distance from them ("z" component).
36 : It is a function of mean squared displacement that are obtained by the joint use of mean squared displacement variables with the SQUARED flag
37 : (see below).
38 :
39 : ## Examples
40 :
41 : Here below is a case where you have defined three frames and you want to
42 : calculate the progress along the path and the distance from it in p1
43 :
44 : ```plumed
45 : #SETTINGS INPUTFILES=regtest/trajectories/path_msd/frame_1.dat,regtest/trajectories/path_msd/frame_21.dat,regtest/trajectories/path_msd/frame_42.dat
46 : t1: RMSD REFERENCE=regtest/trajectories/path_msd/frame_1.dat TYPE=OPTIMAL SQUARED
47 : t2: RMSD REFERENCE=regtest/trajectories/path_msd/frame_21.dat TYPE=OPTIMAL SQUARED
48 : t3: RMSD REFERENCE=regtest/trajectories/path_msd/frame_42.dat TYPE=OPTIMAL SQUARED
49 : p1: FUNCPATHMSD ARG=t1,t2,t3 LAMBDA=500.0
50 : PRINT ARG=t1,t2,t3,p1.s,p1.z STRIDE=1 FILE=colvar FMT=%8.4f
51 : ```
52 :
53 : You can see that the reference coordinates here are described in three separate pdb files that you view above.
54 :
55 : This second example shows how to define a PATH in [CONTACTMAP](CONTACTMAP.md) space:
56 :
57 : ```plumed
58 : CONTACTMAP ...
59 : ATOMS1=1,2 REFERENCE1=0.1
60 : ATOMS2=3,4 REFERENCE2=0.5
61 : ATOMS3=4,5 REFERENCE3=0.25
62 : ATOMS4=5,6 REFERENCE4=0.0
63 : SWITCH={RATIONAL R_0=1.5}
64 : LABEL=c1
65 : CMDIST
66 : ... CONTACTMAP
67 :
68 : CONTACTMAP ...
69 : ATOMS1=1,2 REFERENCE1=0.3
70 : ATOMS2=3,4 REFERENCE2=0.9
71 : ATOMS3=4,5 REFERENCE3=0.45
72 : ATOMS4=5,6 REFERENCE4=0.1
73 : SWITCH={RATIONAL R_0=1.5}
74 : LABEL=c2
75 : CMDIST
76 : ... CONTACTMAP
77 :
78 : CONTACTMAP ...
79 : ATOMS1=1,2 REFERENCE1=1.0
80 : ATOMS2=3,4 REFERENCE2=1.0
81 : ATOMS3=4,5 REFERENCE3=1.0
82 : ATOMS4=5,6 REFERENCE4=1.0
83 : SWITCH={RATIONAL R_0=1.5}
84 : LABEL=c3
85 : CMDIST
86 : ... CONTACTMAP
87 :
88 : p1: FUNCPATHMSD ARG=c1,c2,c3 LAMBDA=500.0
89 : PRINT ARG=c1,c2,c3,p1.s,p1.z STRIDE=1 FILE=colvar FMT=%8.4f
90 : ```
91 :
92 : This third example shows how to define a PATH in [PIV](PIV.md) space:
93 :
94 : ```plumed
95 : #SETTINGS INPUTFILES=regtest/piv/rt-piv-distance/Liq.pdb,regtest/piv/rt-piv-distance/Ice.pdb
96 : PIV ...
97 : LABEL=c1
98 : PRECISION=1000
99 : NLIST
100 : REF_FILE=regtest/piv/rt-piv-distance/Liq.pdb
101 : PIVATOMS=2
102 : ATOMTYPES=A,B
103 : ONLYDIRECT
104 : SFACTOR=1.0,0.2
105 : SORT=1,1
106 : SWITCH1={RATIONAL R_0=0.6 MM=12 NN=4}
107 : SWITCH2={RATIONAL R_0=0.5 MM=10 NN=5}
108 : NL_CUTOFF=1.2,1.2
109 : NL_STRIDE=10,10
110 : NL_SKIN=0.1,0.1
111 : ... PIV
112 : PIV ...
113 : LABEL=c2
114 : PRECISION=1000
115 : NLIST
116 : REF_FILE=regtest/piv/rt-piv-distance/Ice.pdb
117 : PIVATOMS=2
118 : ATOMTYPES=A,B
119 : ONLYDIRECT
120 : SFACTOR=1.0,0.2
121 : SORT=1,1
122 : SWITCH1={RATIONAL R_0=0.6 MM=12 NN=4}
123 : SWITCH2={RATIONAL R_0=0.4 MM=10 NN=5}
124 : NL_CUTOFF=1.2,1.2
125 : NL_STRIDE=10,10
126 : NL_SKIN=0.1,0.1
127 : ... PIV
128 :
129 : p1: FUNCPATHMSD ARG=c1,c2 LAMBDA=0.180338
130 : METAD ARG=p1.s,p1.z SIGMA=0.01,0.2 HEIGHT=0.8 PACE=500 LABEL=res
131 : PRINT ARG=c1,c2,p1.s,p1.z,res.bias STRIDE=500 FILE=colvar FMT=%15.6f
132 : ```
133 :
134 : */
135 : //+ENDPLUMEDOC
136 :
137 : class FuncPathMSD : public Function {
138 : double lambda;
139 : int neigh_size;
140 : double neigh_stride;
141 : std::vector< std::pair<Value *,double> > neighpair;
142 : std::map<Value *,double > indexmap; // use double to allow isomaps
143 : /// the index (as in indexmap) of each element of neighpair, so the map is not searched at every step
144 : std::vector<double> neighindex;
145 : std::vector <Value*> allArguments;
146 : /// rebuild neighindex after neighpair has changed
147 : void updateNeighIndex();
148 : // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
149 : // this below is useful when one wants to sort a vector of double and have back the order
150 : // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
151 : // create a custom sorter
152 : typedef std::vector<double>::const_iterator myiter;
153 : struct ordering {
154 : bool operator ()(std::pair<unsigned, myiter> const& a, std::pair<unsigned, myiter> const& b) {
155 : return *(a.second) < *(b.second);
156 : }
157 : };
158 : // sorting utility
159 : std::vector<int> increasingOrder( std::vector<double> &v) {
160 : // make a pair
161 : std::vector< std::pair<unsigned, myiter> > order(v.size());
162 : unsigned n = 0;
163 : for (myiter it = v.begin(); it != v.end(); ++it, ++n) {
164 : order[n] = make_pair(n, it); // note: heere i do not put the values but the addresses that point to the value
165 : }
166 : // now sort according the second value
167 : std::sort(order.begin(), order.end(), ordering());
168 : std::vector<int> vv(v.size());
169 : n=0;
170 : for (const auto & it : order) {
171 : vv[n]=it.first;
172 : n++;
173 : }
174 : return vv;
175 : }
176 : // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
177 : // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
178 :
179 : struct pairordering {
180 : bool operator ()(std::pair<Value *, double> const& a, std::pair<Value*, double> const& b) {
181 452 : return (a).second > (b).second;
182 : }
183 : };
184 :
185 : public:
186 : explicit FuncPathMSD(const ActionOptions&);
187 : // active methods:
188 : void calculate() override;
189 : void prepare() override;
190 : static void registerKeywords(Keywords& keys);
191 : };
192 :
193 : PLUMED_REGISTER_ACTION(FuncPathMSD,"FUNCPATHMSD")
194 :
195 8 : void FuncPathMSD::registerKeywords(Keywords& keys) {
196 8 : Function::registerKeywords(keys);
197 8 : keys.add("compulsory","LAMBDA","the lambda parameter is needed for smoothing, is in the units of plumed");
198 8 : keys.add("optional","NEIGH_SIZE","size of the neighbor list");
199 8 : keys.add("optional","NEIGH_STRIDE","how often the neighbor list needs to be calculated in time units");
200 16 : keys.addOutputComponent("s","default","scalar","the position on the path");
201 16 : keys.addOutputComponent("z","default","scalar","the distance from the path");
202 8 : keys.addDOI("10.1063/1.2432340");
203 8 : }
204 4 : FuncPathMSD::FuncPathMSD(const ActionOptions&ao):
205 : Action(ao),
206 : Function(ao),
207 4 : neigh_size(-1),
208 4 : neigh_stride(-1.) {
209 :
210 4 : parse("LAMBDA",lambda);
211 4 : parse("NEIGH_SIZE",neigh_size);
212 4 : parse("NEIGH_STRIDE",neigh_stride);
213 4 : checkRead();
214 4 : log.printf(" lambda is %f\n",lambda);
215 4 : if( !(lambda>0) ) {
216 0 : error("LAMBDA should be positive");
217 : }
218 4 : if( neigh_size>0 && !(neigh_stride>0) ) {
219 0 : error("with NEIGH_SIZE you should also give a positive NEIGH_STRIDE");
220 : }
221 : // list the action involved and check the type
222 4 : std::string myname=getPntrToArgument(0)->getPntrToAction()->getName();
223 : // DISTANCE is a shortcut: the scalar distance is computed by DISTANCE_SCALAR
224 4 : if(myname!="RMSD_SCALAR"&&myname!="CONTACTMAP"&&myname!="DISTANCE"&&myname!="DISTANCE_SCALAR"&&myname!="PIV") {
225 0 : error("One or more of your arguments is not of RMSD/CONTACTMAP/DISTANCE/PIV type!!!");
226 : }
227 14 : for(unsigned i=1; i<getNumberOfArguments(); i++) {
228 : // for each value get the name and the label of the corresponding action
229 10 : if( getPntrToArgument(i)->getPntrToAction()->getName()!=myname ) {
230 0 : error("mismatch between the types of arguments");
231 : }
232 : }
233 4 : log.printf(" Consistency check completed! Your path cvs look good!\n");
234 : // do some neighbor printout
235 4 : if(neigh_size>0) {
236 2 : if(neigh_size>static_cast<int>(getNumberOfArguments())) {
237 0 : log.printf(" List size required ( %d ) is too large: resizing to the maximum number of arg required: %d \n",neigh_size,getNumberOfArguments());
238 0 : neigh_size=getNumberOfArguments();
239 : }
240 2 : log.printf(" Neighbor list enabled: \n");
241 2 : log.printf(" size : %d elements\n",neigh_size);
242 2 : log.printf(" stride : %f time \n",neigh_stride);
243 : } else {
244 2 : log.printf(" Neighbor list NOT enabled \n");
245 : }
246 :
247 8 : addComponentWithDerivatives("s");
248 8 : componentIsNotPeriodic("s");
249 8 : addComponentWithDerivatives("z");
250 4 : componentIsNotPeriodic("z");
251 :
252 : // now backup the arguments
253 18 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
254 14 : allArguments.push_back(getPntrToArgument(i));
255 : }
256 : double i=1.;
257 18 : for(const auto & it : allArguments) {
258 14 : indexmap[it]=i;
259 14 : i+=1.;
260 : }
261 :
262 4 : }
263 1102 : void FuncPathMSD::updateNeighIndex() {
264 1102 : neighindex.resize(neighpair.size());
265 4005 : for(unsigned i=0; i<neighpair.size(); ++i) {
266 2903 : neighindex[i]=indexmap[neighpair[i].first];
267 : }
268 1102 : }
269 :
270 : // calculator
271 1102 : void FuncPathMSD::calculate() {
272 : double s_path=0.;
273 : double partition=0.;
274 1102 : if(neighpair.empty()) { // at first step, resize it
275 0 : neighpair.resize(allArguments.size());
276 0 : for(unsigned i=0; i<allArguments.size(); i++) {
277 0 : neighpair[i].first=allArguments[i];
278 : }
279 : }
280 1102 : if( neighindex.size()!=neighpair.size() ) {
281 0 : updateNeighIndex();
282 : }
283 :
284 1102 : Value* val_s_path=getPntrToComponent("s");
285 2204 : Value* val_z_path=getPntrToComponent("z");
286 :
287 : // The weights are computed relative to the closest frame, exp(-lambda*(d_i-d_min)),
288 : // so that they do not all underflow to zero far from the path.
289 1102 : double dmin=neighpair[0].first->get();
290 4005 : for(const auto & it : neighpair) {
291 3624 : dmin=std::min(dmin,it.first->get());
292 : }
293 4005 : for(unsigned i=0; i<neighpair.size(); ++i) {
294 : auto & it=neighpair[i];
295 2903 : it.second=std::exp(-lambda*(it.first->get()-dmin));
296 2903 : s_path+=neighindex[i]*it.second;
297 2903 : partition+=it.second;
298 : }
299 : // the closest frame has weight one, so partition>=1; the check makes this explicit
300 : // (and keeps static analysis from flagging the logarithm below)
301 1102 : if( !(partition>=1.0) ) {
302 0 : error("the weights of the frames of the path are not finite: check LAMBDA and the arguments");
303 : return;
304 : }
305 1102 : s_path/=partition;
306 : val_s_path->set(s_path);
307 1102 : val_z_path->set(dmin-(1./lambda)*std::log(partition));
308 4005 : for(unsigned n=0; n<neighpair.size(); ++n) {
309 2903 : double expval=neighpair[n].second;
310 2903 : double tmp=lambda*expval*(s_path-neighindex[n])/partition;
311 : setDerivative(val_s_path,n,tmp);
312 2903 : setDerivative(val_z_path,n,expval/partition);
313 : }
314 : }
315 : ///
316 : /// this function updates the needed argument list
317 : ///
318 1102 : void FuncPathMSD::prepare() {
319 :
320 : // neighbor list: rank and activate the chain for the next step
321 :
322 : // neighbor list: if neigh_size<0 never sort and keep the full vector
323 : // neighbor list: if neigh_size>0
324 : // if the size is full -> sort the vector and decide the dependencies for next step
325 : // if the size is not full -> check if next step will need the full dependency otherwise keep this dependencies
326 :
327 : // here just resize the neighpair. The real resizing of reinit will be done by the prepare stage that will modify the list of arguments
328 1102 : if (neigh_size>0) {
329 551 : if(neighpair.size()==allArguments.size()) { // I just did the complete round: need to sort, shorten and give it a go
330 : // sort the values
331 139 : std::sort(neighpair.begin(),neighpair.end(),pairordering());
332 : // resize the effective list
333 139 : neighpair.resize(neigh_size);
334 139 : log.printf(" NEIGH LIST NOW INCLUDE INDEXES: ");
335 417 : for(int i=0; i<neigh_size; ++i) {
336 278 : log.printf(" %f ",indexmap[neighpair[i].first]);
337 : }
338 139 : log.printf(" \n");
339 : } else {
340 412 : const long long nstride=std::max(1LL,static_cast<long long>(neigh_stride/getTimeStep()));
341 412 : if( getStep()%nstride==0 ) {
342 140 : log.printf(" Time %f : recalculating full neighlist \n",getStep()*getTimeStep());
343 140 : neighpair.resize(allArguments.size());
344 563 : for(unsigned i=0; i<allArguments.size(); i++) {
345 423 : neighpair[i].first=allArguments[i];
346 : }
347 : }
348 : }
349 : } else {
350 551 : if( int(getStep())==0) {
351 2 : neighpair.resize(allArguments.size());
352 9 : for(unsigned i=0; i<allArguments.size(); i++) {
353 7 : neighpair[i].first=allArguments[i];
354 : }
355 : }
356 : }
357 1102 : std::vector<Value*> argstocall;
358 1102 : if(!neighpair.empty()) {
359 4005 : for(const auto & it : neighpair) {
360 2903 : argstocall.push_back( it.first );
361 : }
362 : } else {
363 0 : argstocall=allArguments;
364 : }
365 1102 : updateNeighIndex();
366 : // the list of arguments (and the dependencies) only needs to be updated when it changes
367 1102 : if( argstocall!=getArguments() ) {
368 277 : requestArguments(argstocall);
369 831 : for(unsigned i=0; i< getNumberOfComponents(); i++) {
370 554 : getPntrToComponent(i)->clearDerivatives();
371 554 : getPntrToComponent(i)->resizeDerivatives(getNumberOfArguments());
372 : }
373 : }
374 1102 : }
375 :
376 : }
377 : }
378 :
379 :
|
