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