Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 <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 MSD that are obtained by the joint use of MSD variable and 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 alng the path and the distance from it in p1
51 :
52 : \verbatim
53 : t1: RMSD REFERENCE=frame_1.dat TYPE=OPTIMAL SQUARED
54 : t2: RMSD REFERENCE=frame_21.dat TYPE=OPTIMAL SQUARED
55 : t3: RMSD REFERENCE=frame_42.dat 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 : \endverbatim
59 :
60 : In this second example is shown how to define a PATH in the \ref CONTACTMAP space:
61 :
62 : \verbatim
63 : CONTACTMAP ...
64 : ATOMS1=1,2 REFERENCE1=0.1
65 : ATOMS2=3,4 REFERENCE2=0.5
66 : ATOMS3=4,5 REFERENCE3=0.25
67 : ATOMS4=5,6 REFERENCE4=0.0
68 : SWITCH={RATIONAL R_0=1.5}
69 : LABEL=c1
70 : CMDIST
71 : ... CONTACTMAP
72 :
73 : CONTACTMAP ...
74 : ATOMS1=1,2 REFERENCE1=0.3
75 : ATOMS2=3,4 REFERENCE2=0.9
76 : ATOMS3=4,5 REFERENCE3=0.45
77 : ATOMS4=5,6 REFERENCE4=0.1
78 : SWITCH={RATIONAL R_0=1.5}
79 : LABEL=c2
80 : CMDIST
81 : ... CONTACTMAP
82 :
83 : CONTACTMAP ...
84 : ATOMS1=1,2 REFERENCE1=1.0
85 : ATOMS2=3,4 REFERENCE2=1.0
86 : ATOMS3=4,5 REFERENCE3=1.0
87 : ATOMS4=5,6 REFERENCE4=1.0
88 : SWITCH={RATIONAL R_0=1.5}
89 : LABEL=c3
90 : CMDIST
91 : ... CONTACTMAP
92 :
93 : p1: FUNCPATHMSD ARG=c1,c2,c3 LAMBDA=500.0
94 : PRINT ARG=c1,c2,c3,p1.s,p1.z STRIDE=1 FILE=colvar FMT=%8.4f
95 : \endverbatim
96 :
97 : */
98 : //+ENDPLUMEDOC
99 :
100 4 : class FuncPathMSD : public Function {
101 : double lambda;
102 : int neigh_size;
103 : double neigh_stride;
104 : vector< pair<Value *,double> > neighpair;
105 : map<Value *,double > indexmap; // use double to allow isomaps
106 : vector <Value*> allArguments;
107 : // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
108 : // this below is useful when one wants to sort a vector of double and have back the order
109 : // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
110 : // create a custom sorter
111 : typedef vector<double>::const_iterator myiter;
112 : struct ordering {
113 : bool operator ()(pair<unsigned, myiter> const& a, pair<unsigned, myiter> const& b) {
114 : return *(a.second) < *(b.second);
115 : }
116 : };
117 : // sorting utility
118 : vector<int> increasingOrder( vector<double> &v) {
119 : // make a pair
120 : vector< pair<unsigned, myiter> > order(v.size());
121 : unsigned n = 0;
122 : for (myiter it = v.begin(); it != v.end(); ++it, ++n) {
123 : order[n] = make_pair(n, it); // note: heere i do not put the values but the addresses that point to the value
124 : }
125 : // now sort according the second value
126 : sort(order.begin(), order.end(), ordering());
127 : typedef vector< pair<unsigned, myiter> >::const_iterator pairiter;
128 : vector<int> vv(v.size()); n=0;
129 : for ( pairiter it = order.begin(); it != order.end(); ++it ) {
130 : vv[n]=(*it).first; n++;
131 : }
132 : return vv;
133 : }
134 : // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
135 : // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
136 :
137 : struct pairordering {
138 437 : bool operator ()(pair<Value *, double> const& a, pair<Value*, double> const& b) {
139 437 : return (a).second > (b).second;
140 : }
141 : };
142 :
143 : public:
144 : explicit FuncPathMSD(const ActionOptions&);
145 : // active methods:
146 : virtual void calculate();
147 : virtual void prepare();
148 : static void registerKeywords(Keywords& keys);
149 : };
150 :
151 2525 : PLUMED_REGISTER_ACTION(FuncPathMSD,"FUNCPATHMSD")
152 :
153 3 : void FuncPathMSD::registerKeywords(Keywords& keys) {
154 3 : Function::registerKeywords(keys);
155 3 : keys.use("ARG");
156 3 : keys.add("compulsory","LAMBDA","the lambda parameter is needed for smoothing, is in the units of plumed");
157 3 : keys.add("optional","NEIGH_SIZE","size of the neighbor list");
158 3 : keys.add("optional","NEIGH_STRIDE","how often the neighbor list needs to be calculated in time units");
159 3 : componentsAreNotOptional(keys);
160 3 : keys.addOutputComponent("s","default","the position on the path");
161 3 : keys.addOutputComponent("z","default","the distance from the path");
162 3 : }
163 2 : FuncPathMSD::FuncPathMSD(const ActionOptions&ao):
164 : Action(ao),
165 : Function(ao),
166 : neigh_size(-1),
167 2 : neigh_stride(-1.)
168 : {
169 :
170 2 : parse("LAMBDA",lambda);
171 2 : parse("NEIGH_SIZE",neigh_size);
172 2 : parse("NEIGH_STRIDE",neigh_stride);
173 2 : checkRead();
174 2 : log.printf(" lambda is %f\n",lambda);
175 : // list the action involved and check the type
176 2 : std::string myname=getPntrToArgument(0)->getPntrToAction()->getName();
177 2 : if(myname!="RMSD"&&myname!="CONTACTMAP"&&myname!="DISTANCE") error("One or more of your arguments is not of RMSD/CONTACTMAP/DISTANCE type!!!");
178 6 : for(unsigned i=1; i<getNumberOfArguments(); i++) {
179 : // for each value get the name and the label of the corresponding action
180 4 : if( getPntrToArgument(i)->getPntrToAction()->getName()!=myname ) error("mismatch between the types of arguments");
181 : }
182 2 : log.printf(" Consistency check completed! Your path cvs look good!\n");
183 : // do some neighbor printout
184 2 : if(neigh_stride>0. || neigh_size>0) {
185 1 : if(neigh_size>getNumberOfArguments()) {
186 0 : log.printf(" List size required ( %d ) is too large: resizing to the maximum number of arg required: %d \n",neigh_size,getNumberOfArguments());
187 0 : neigh_size=getNumberOfArguments();
188 : }
189 1 : log.printf(" Neighbor list enabled: \n");
190 1 : log.printf(" size : %d elements\n",neigh_size);
191 1 : log.printf(" stride : %f time \n",neigh_stride);
192 : } else {
193 1 : log.printf(" Neighbor list NOT enabled \n");
194 : }
195 :
196 2 : addComponentWithDerivatives("s"); componentIsNotPeriodic("s");
197 2 : addComponentWithDerivatives("z"); componentIsNotPeriodic("z");
198 :
199 : // now backup the arguments
200 2 : for(unsigned i=0; i<getNumberOfArguments(); i++)allArguments.push_back(getPntrToArgument(i));
201 2 : double i=1.;
202 8 : for(std::vector<Value*>:: const_iterator it=allArguments.begin(); it!=allArguments.end() ; ++it) {
203 6 : indexmap[(*it)]=i; i+=1.;
204 2 : }
205 :
206 2 : }
207 : // calculator
208 1092 : void FuncPathMSD::calculate() {
209 : // log.printf("NOW CALCULATE! \n");
210 1092 : double s_path=0.;
211 1092 : double partition=0.;
212 1092 : if(neighpair.empty()) { // at first step, resize it
213 0 : neighpair.resize(allArguments.size());
214 0 : for(unsigned i=0; i<allArguments.size(); i++)neighpair[i].first=allArguments[i];
215 : }
216 :
217 1092 : Value* val_s_path=getPntrToComponent("s");
218 1092 : Value* val_z_path=getPntrToComponent("z");
219 :
220 : typedef vector< pair< Value *,double> >::iterator pairiter;
221 3959 : for(pairiter it=neighpair.begin(); it!=neighpair.end(); ++it) {
222 2867 : (*it).second=exp(-lambda*((*it).first->get()));
223 2867 : s_path+=(indexmap[(*it).first])*(*it).second;
224 2867 : partition+=(*it).second;
225 : }
226 1092 : s_path/=partition;
227 1092 : val_s_path->set(s_path);
228 1092 : val_z_path->set(-(1./lambda)*std::log(partition));
229 1092 : int n=0;
230 3959 : for(pairiter it=neighpair.begin(); it!=neighpair.end(); ++it) {
231 2867 : double expval=(*it).second;
232 2867 : double tmp=lambda*expval*(s_path-(indexmap[(*it).first]))/partition;
233 2867 : setDerivative(val_s_path,n,tmp);
234 2867 : setDerivative(val_z_path,n,expval/partition);
235 2867 : n++;
236 : }
237 :
238 : // log.printf("CALCULATION DONE! \n");
239 1092 : }
240 : ///
241 : /// this function updates the needed argument list
242 : ///
243 1092 : void FuncPathMSD::prepare() {
244 :
245 : // neighbor list: rank and activate the chain for the next step
246 :
247 : // neighbor list: if neigh_size<0 never sort and keep the full vector
248 : // neighbor list: if neigh_size>0
249 : // if the size is full -> sort the vector and decide the dependencies for next step
250 : // if the size is not full -> check if next step will need the full dependency otherwise keep this dependencies
251 :
252 : // here just resize the neighpair. The real resizing of reinit will be done by the prepare stage that will modify the list of arguments
253 1092 : if (neigh_size>0) {
254 546 : if(neighpair.size()==allArguments.size()) { // I just did the complete round: need to sort, shorten and give it a go
255 : // sort the values
256 137 : sort(neighpair.begin(),neighpair.end(),pairordering());
257 : // resize the effective list
258 137 : neighpair.resize(neigh_size);
259 137 : log.printf(" NEIGH LIST NOW INCLUDE INDEXES: ");
260 137 : for(int i=0; i<neigh_size; ++i)log.printf(" %f ",indexmap[neighpair[i].first]); log.printf(" \n");
261 : } else {
262 409 : if( int(getStep())%int(neigh_stride/getTimeStep())==0 ) {
263 137 : log.printf(" Time %f : recalculating full neighlist \n",getStep()*getTimeStep());
264 137 : neighpair.resize(allArguments.size());
265 137 : for(unsigned i=0; i<allArguments.size(); i++)neighpair[i].first=allArguments[i];
266 : }
267 : }
268 : } else {
269 546 : if( int(getStep())==0) {
270 1 : neighpair.resize(allArguments.size());
271 1 : for(unsigned i=0; i<allArguments.size(); i++)neighpair[i].first=allArguments[i];
272 : }
273 : }
274 : typedef vector< pair<Value*,double> >::iterator pairiter;
275 1092 : vector<Value*> argstocall;
276 : //log.printf("PREPARING \n");
277 1092 : argstocall.clear();
278 1092 : if(!neighpair.empty()) {
279 3959 : for(pairiter it=neighpair.begin(); it!=neighpair.end(); ++it) {
280 2867 : argstocall.push_back( (*it).first );
281 : // log.printf("CALLING %p %f ",(*it).first ,indexmap[(*it).first] );
282 : }
283 : } else {
284 0 : for(unsigned i=0; i<allArguments.size(); i++) {
285 0 : argstocall.push_back(allArguments[i]);
286 : }
287 : }
288 : // now the list of argument changes
289 1092 : requestArguments(argstocall);
290 : //now resize the derivatives as well
291 : //for each value in this action
292 3276 : for(int i=0; i< getNumberOfComponents(); i++) {
293 : //resize the derivative to the number the
294 2184 : getPntrToComponent(i)->clearDerivatives();
295 2184 : getPntrToComponent(i)->resizeDerivatives(getNumberOfArguments());
296 1092 : }
297 : //log.printf("PREPARING DONE! \n");
298 1092 : }
299 :
300 : }
301 2523 : }
302 :
303 :
|
