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