Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2020-2021 of Michele Invernizzi.
3 :
4 : This file is part of the OPES plumed module.
5 :
6 : The OPES plumed module is free software: you can redistribute it and/or modify
7 : it under the terms of the GNU Lesser General Public License as published by
8 : the Free Software Foundation, either version 3 of the License, or
9 : (at your option) any later version.
10 :
11 : The OPES plumed module is distributed in the hope that it will be useful,
12 : but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : GNU Lesser General Public License for more details.
15 :
16 : You should have received a copy of the GNU Lesser General Public License
17 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
18 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
19 : #include "bias/Bias.h"
20 : #include "core/PlumedMain.h"
21 : #include "core/ActionRegister.h"
22 : #include "core/ActionSet.h"
23 : #include "tools/Communicator.h"
24 : #include "tools/File.h"
25 : #include "tools/OpenMP.h"
26 :
27 : #include "ExpansionCVs.h"
28 :
29 : namespace PLMD {
30 : namespace opes {
31 :
32 : //+PLUMEDOC BIAS OPES_EXPANDED
33 : /*
34 : On-the-fly probability enhanced sampling with expanded ensembles for the target distribution.
35 :
36 : This method is similar to the [OPES method](module_opes.md) with expanded ensembles target distribution (replica-exchange-like) that is discussed in the paper cited below.
37 :
38 : An expanded ensemble is obtained by summing a set of ensembles at slightly different termodynamic conditions, or with slightly different Hamiltonians.
39 : Such ensembles can be sampled via methods like replica exchange, or this [OPES_EXPANDED](OPES_EXPANDED.md) bias action.
40 : A typical example is a multicanonical simulation, in which a whole range of temperatures is sampled instead of a single one.
41 :
42 : In oreder to define an expanded target ensemble we use the EXPANSION_CV (ECVs) that are discussed on the [module page](module_opes.md), $\Delta u_\lambda(\mathbf{x})$.
43 : The bias at step $n$ is
44 :
45 : $$
46 : V_n(\mathbf{x})=-\frac{1}{\beta}\log \left(\frac{1}{N_{\{\lambda\}}}\sum_\lambda e^{-\Delta u_\lambda(\mathbf{x})+\beta\Delta F_n(\lambda)}\right)\, .
47 : $$
48 :
49 : See the paper cited below for more details on the method.
50 :
51 : Notice that the estimates in the DELTAFS file are expressed in energy units, and should be multiplied by $\beta$ to be dimensionless as in the paper cited below.
52 : The DELTAFS file also contains an estimate of $c(t)=\frac{1}{\beta} \log \langle e^{\beta V}\rangle$.
53 : Similarly to [OPES_METAD](OPES_METAD.md), it is printed only for reference, since $c(t)$ reaches a fixed value as the bias converges, and should NOT be used for reweighting.
54 : Its value is also needed for restarting a simulation.
55 :
56 : You can store the instantaneous $\Delta F_n(\lambda)$ estimates also in a more readable format using STATE_WFILE and STATE_WSTRIDE.
57 : Restart can be done either from a DELTAFS file or from a STATE_RFILE, it is equivalent.
58 : To continue a simulation with a static bias, for example once the $\Delta F(\lambda)$ estimates have converged, restart and add UPDATE_UNTIL=0.
59 : The $\Delta F(\lambda)$ estimates are then no longer updated, and the bias remains the one read at restart.
60 :
61 : Contrary to [OPES_METAD](OPES_METAD.md), OPES_EXPANDED does not use kernel density estimation.
62 :
63 : ## Examples
64 :
65 : ```plumed
66 : # simulate multiple temperatures, as in parallel tempering
67 : ene: ENERGY
68 : ecv: ECV_MULTITHERMAL ARG=ene TEMP_MAX=1000
69 : opes: OPES_EXPANDED ARG=ecv.* PACE=500
70 : PRINT FILE=COLVAR STRIDE=500 ARG=ene,opes.bias
71 : ```
72 :
73 : You can easily combine multiple ECVs.
74 : The OPES_EXPANDED bias will create a multidimensional target grid to sample all the combinations.
75 :
76 : ```plumed
77 : # simulate multiple temperatures while biasing a CV
78 : ene: ENERGY
79 : dst: DISTANCE ATOMS=1,2
80 :
81 : ecv1: ECV_MULTITHERMAL ARG=ene TEMP_SET_ALL=200,300,500,1000
82 : ecv2: ECV_UMBRELLAS_LINE ARG=dst CV_MIN=1.2 CV_MAX=4.3 SIGMA=0.5
83 : opes: OPES_EXPANDED ARG=ecv1.*,ecv2.* PACE=500 OBSERVATION_STEPS=1
84 :
85 : PRINT FILE=COLVAR STRIDE=500 ARG=ene,dst,opes.bias
86 : ```
87 :
88 : If an ECV is based on more than one CV you must provide all the output component, in the proper order.
89 : You can use the regular expressions discussed [here](specifying_arguments.md) for that, like in the following example.
90 :
91 : ```plumed
92 : # simulate multiple temperatures and pressures while biasing a two-CVs linear path
93 : ene: ENERGY
94 : vol: VOLUME
95 : ecv_mtp: ECV_MULTITHERMAL_MULTIBARIC ...
96 : ARG=ene,vol
97 : TEMP=300
98 : TEMP_MIN=200
99 : TEMP_MAX=800
100 : PRESSURE=0.06022140857*1000 #1 kbar
101 : PRESSURE_MIN=0
102 : PRESSURE_MAX=0.06022140857*2000 #2 kbar
103 : ...
104 :
105 : cv1: DISTANCE ATOMS=1,2
106 : cv2: DISTANCE ATOMS=3,4
107 : ecv_umb: ECV_UMBRELLAS_LINE ARG=cv1,cv2 TEMP=300 CV_MIN=0.1,0.1 CV_MAX=1.5,1.5 SIGMA=0.2 BARRIER=70
108 :
109 : opes: OPES_EXPANDED ARG=(ecv_.*) PACE=500 WALKERS_MPI PRINT_STRIDE=1000
110 :
111 : PRINT FILE=COLVAR STRIDE=500 ARG=ene,vol,cv1,cv2,opes.bias
112 : ```
113 :
114 : */
115 : //+ENDPLUMEDOC
116 :
117 : class OPESexpanded : public bias::Bias {
118 :
119 : private:
120 : bool isFirstStep_;
121 : unsigned NumOMP_;
122 : unsigned NumParallel_;
123 : unsigned rank_;
124 : unsigned NumWalkers_;
125 : unsigned walker_rank_;
126 : unsigned long long counter_;
127 : std::size_t ncv_;
128 :
129 : std::vector<const double *> ECVs_;
130 : std::vector<const double *> derECVs_;
131 : std::vector<opes::ExpansionCVs*> pntrToECVsClass_;
132 : std::vector< std::vector<unsigned> > index_k_;
133 : // A note on indexes usage:
134 : // j -> underlying CVs
135 : // i -> DeltaFs
136 : // k -> single ECVs, which might not be trivially numbered
137 : // l -> groups of ECVs, pntrToECVsClass
138 : // h -> subgroups of ECVs, arguments in ECVsClass
139 : // w -> walkers
140 :
141 : double kbt_;
142 : unsigned stride_;
143 : unsigned deltaF_size_; //different from deltaF_.size() if NumParallel_>1
144 : std::vector<double> deltaF_;
145 : std::vector<double> diff_;
146 : double rct_;
147 :
148 : std::vector<double> all_deltaF_;
149 : std::vector<int> all_size_;
150 : std::vector<int> disp_;
151 :
152 : unsigned obs_steps_;
153 : std::vector<double> obs_cvs_;
154 :
155 : bool calc_work_;
156 : double work_;
157 :
158 : unsigned print_stride_;
159 : OFile deltaFsOfile_;
160 : std::vector<std::string> deltaF_name_;
161 :
162 : OFile stateOfile_;
163 : int wStateStride_;
164 : bool storeOldStates_;
165 :
166 : void init_pntrToECVsClass();
167 : void init_linkECVs();
168 : void init_fromObs();
169 :
170 : void printDeltaF();
171 : void dumpStateToFile();
172 : void updateDeltaF(double);
173 : double getExpansion(const unsigned) const;
174 :
175 : public:
176 : explicit OPESexpanded(const ActionOptions&);
177 : void calculate() override;
178 : void update() override;
179 : static void registerKeywords(Keywords& keys);
180 : };
181 :
182 : PLUMED_REGISTER_ACTION(OPESexpanded,"OPES_EXPANDED")
183 :
184 40 : void OPESexpanded::registerKeywords(Keywords& keys) {
185 40 : Bias::registerKeywords(keys);
186 80 : keys.addInputKeyword("compulsory","ARG","scalar","the label of the ECVs that define the expansion. You can use an * to make sure all the output components of the ECVs are used, as in the examples above");
187 40 : keys.add("compulsory","PACE","how often the bias is updated");
188 40 : keys.add("compulsory","OBSERVATION_STEPS","100","number of unbiased initial PACE steps to collect statistics for initialization");
189 : //DeltaFs and state files
190 40 : keys.add("compulsory","FILE","DELTAFS","a file with the estimate of the relative Delta F for each component of the target and of the global c(t)");
191 40 : keys.add("compulsory","PRINT_STRIDE","100","stride for printing to DELTAFS file, in units of PACE");
192 40 : keys.add("optional","FMT","specify format for DELTAFS file");
193 40 : keys.add("optional","STATE_RFILE","read from this file the Delta F estimates and all the info needed to RESTART the simulation");
194 40 : keys.add("optional","STATE_WFILE","write to this file the Delta F estimates and all the info needed to RESTART the simulation");
195 40 : keys.add("optional","STATE_WSTRIDE","number of MD steps between writing the STATE_WFILE. Default is only on CPT events (but not all MD codes set them)");
196 40 : keys.addFlag("STORE_STATES",false,"append to STATE_WFILE instead of overwriting it each time");
197 : //miscellaneous
198 40 : keys.addFlag("CALC_WORK",false,"calculate the total accumulated work done by the bias since last restart");
199 40 : keys.addFlag("WALKERS_MPI",false,"switch on MPI version of multiple walkers");
200 40 : keys.addFlag("SERIAL",false,"perform calculations in serial");
201 40 : keys.use("RESTART");
202 40 : keys.use("UPDATE_FROM");
203 40 : keys.use("UPDATE_UNTIL");
204 :
205 : //output components
206 80 : keys.addOutputComponent("work","CALC_WORK","scalar","total accumulated work done by the bias");
207 40 : keys.addDOI("10.1103/PhysRevX.10.041034");
208 40 : }
209 :
210 36 : OPESexpanded::OPESexpanded(const ActionOptions&ao)
211 : : PLUMED_BIAS_INIT(ao)
212 36 : , isFirstStep_(true)
213 36 : , counter_(0)
214 36 : , ncv_(getNumberOfArguments())
215 36 : , deltaF_size_(0)
216 36 : , rct_(0)
217 36 : , work_(0) {
218 : //set pace
219 36 : long long pace=0;
220 36 : parse("PACE",pace);
221 36 : plumed_massert(pace>0 && static_cast<unsigned long long>(pace)<=std::numeric_limits<unsigned>::max(),
222 : "PACE must be greater than zero and no larger than "+std::to_string(std::numeric_limits<unsigned>::max()));
223 36 : stride_=static_cast<unsigned>(pace);
224 36 : parse("OBSERVATION_STEPS",obs_steps_);
225 36 : plumed_massert(obs_steps_!=0,"minimum is OBSERVATION_STEPS=1");
226 36 : obs_cvs_.resize(obs_steps_*ncv_);
227 :
228 : //deltaFs file
229 : std::string deltaFsFileName;
230 36 : parse("FILE",deltaFsFileName);
231 72 : parse("PRINT_STRIDE",print_stride_);
232 : std::string fmt;
233 72 : parse("FMT",fmt);
234 : //output checkpoint of current state
235 : std::string restartFileName;
236 72 : parse("STATE_RFILE",restartFileName);
237 : std::string stateFileName;
238 36 : parse("STATE_WFILE",stateFileName);
239 36 : wStateStride_=0;
240 36 : parse("STATE_WSTRIDE",wStateStride_);
241 36 : storeOldStates_=false;
242 36 : parseFlag("STORE_STATES",storeOldStates_);
243 36 : if(wStateStride_!=0 || storeOldStates_) {
244 7 : plumed_massert(stateFileName.length()>0,"filename for storing simulation status not specified, use STATE_WFILE");
245 : }
246 36 : if(wStateStride_>0) {
247 7 : plumed_massert(wStateStride_>=(int)stride_,"STATE_WSTRIDE is in units of MD steps, thus should be a multiple of PACE");
248 : }
249 36 : if(stateFileName.length()>0 && wStateStride_==0) {
250 1 : wStateStride_=-1; //will print only on CPT events (checkpoints set by some MD engines, like gromacs)
251 : }
252 :
253 : //work flag
254 36 : parseFlag("CALC_WORK",calc_work_);
255 :
256 : //multiple walkers //external MW for cp2k not supported, but anyway cp2k cannot put bias on energy!
257 36 : bool walkers_mpi=false;
258 36 : parseFlag("WALKERS_MPI",walkers_mpi);
259 36 : if(walkers_mpi) {
260 : //If this Action is not compiled with MPI the user is informed and we exit gracefully
261 6 : plumed_massert(Communicator::plumedHasMPI(),"Invalid walkers configuration: WALKERS_MPI flag requires MPI compilation");
262 6 : plumed_massert(Communicator::initialized(),"Invalid walkers configuration: WALKERS_MPI needs the communicator correctly initialized.");
263 :
264 6 : NumWalkers_=getNumberOfReplicas();
265 6 : walker_rank_=getReplicaIndex();
266 : } else {
267 30 : NumWalkers_=1;
268 30 : walker_rank_=0;
269 : }
270 :
271 : //parallelization stuff
272 36 : NumOMP_=OpenMP::getNumThreads();
273 36 : NumParallel_=comm.Get_size();
274 36 : rank_=comm.Get_rank();
275 36 : bool serial=false;
276 36 : parseFlag("SERIAL",serial);
277 36 : if(serial) {
278 5 : NumOMP_=1;
279 5 : NumParallel_=1;
280 5 : rank_=0;
281 : }
282 :
283 36 : checkRead();
284 :
285 : //check ECVs and link them
286 36 : init_pntrToECVsClass();
287 : //set kbt_
288 36 : kbt_=pntrToECVsClass_[0]->getKbT();
289 81 : for(unsigned l=0; l<pntrToECVsClass_.size(); l++) {
290 45 : plumed_massert(std::abs(kbt_-pntrToECVsClass_[l]->getKbT())<1e-4,"must set same TEMP for each ECV");
291 : }
292 :
293 : //restart if needed
294 36 : if(getRestart()) {
295 : bool stateRestart=true;
296 12 : if(restartFileName.length()==0) {
297 : stateRestart=false;
298 : restartFileName=deltaFsFileName;
299 : }
300 12 : IFile ifile;
301 12 : ifile.link(*this);
302 12 : if(ifile.FileExist(restartFileName)) {
303 12 : log.printf(" RESTART - make sure all ECVs used are the same as before\n");
304 12 : log.printf(" restarting from: %s\n",restartFileName.c_str());
305 12 : ifile.open(restartFileName);
306 12 : if(stateRestart) { //get all info
307 3 : log.printf(" it should be a STATE file (not a DELTAFS file)\n");
308 : double time; //not used
309 3 : ifile.scanField("time",time);
310 3 : ifile.scanField("counter",counter_);
311 6 : ifile.scanField("rct",rct_);
312 : std::string tmp_lambda;
313 69 : while(ifile.scanField(getPntrToArgument(0)->getName(),tmp_lambda)) {
314 66 : std::string subs="DeltaF_"+tmp_lambda;
315 132 : for(unsigned jj=1; jj<ncv_; jj++) {
316 : tmp_lambda.clear();
317 66 : ifile.scanField(getPntrToArgument(jj)->getName(),tmp_lambda);
318 132 : subs+="_"+tmp_lambda;
319 : }
320 66 : deltaF_name_.push_back(subs);
321 : double tmp_deltaF;
322 66 : ifile.scanField("DeltaF",tmp_deltaF);
323 66 : deltaF_.push_back(tmp_deltaF);
324 66 : ifile.scanField();
325 : tmp_lambda.clear();
326 : }
327 3 : log.printf(" successfully read %lu DeltaF values\n",deltaF_name_.size());
328 3 : if(NumParallel_>1) {
329 2 : all_deltaF_=deltaF_;
330 : }
331 : } else { //get just deltaFs names
332 9 : ifile.scanFieldList(deltaF_name_);
333 9 : plumed_massert(deltaF_name_.size()>=4,"RESTART - fewer than expected FIELDS found in '"+deltaFsFileName+"' file");
334 9 : plumed_massert(deltaF_name_[deltaF_name_.size()-1]=="print_stride","RESTART - could not find expected FIELDS in '"+deltaFsFileName+"' file");
335 9 : plumed_massert(deltaF_name_[0]=="time","RESTART - could not find expected FIELDS in '"+deltaFsFileName+"' file");
336 9 : plumed_massert(deltaF_name_[1]=="rct","RESTART - could not find expected FIELDS in '"+deltaFsFileName+"' file");
337 : deltaF_name_.pop_back();
338 : deltaF_name_.erase(deltaF_name_.begin(),deltaF_name_.begin()+2);
339 : std::size_t pos=5; //each name starts with "DeltaF"
340 25 : for(unsigned j=0; j<ncv_; j++) {
341 16 : pos=deltaF_name_[0].find("_",pos+1); //checking only first one, hopefully is enough
342 : }
343 9 : plumed_massert(pos<deltaF_name_[0].length(),"RESTART - fewer '_' than expected in DeltaF fields: did you remove any CV?");
344 9 : pos=deltaF_name_[0].find("_",pos+1);
345 9 : plumed_massert(pos>deltaF_name_[0].length(),"RESTART - more '_' than expected in DeltaF fields: did you add new CV?");
346 : }
347 : //get lambdas, init ECVs and Link them
348 12 : deltaF_size_=deltaF_name_.size();
349 517 : auto getLambdaName=[](const std::string& name,const unsigned start,const unsigned dim) {
350 : std::size_t pos_start=5; //each name starts with "DeltaF"
351 1052 : for(unsigned j=0; j<=start; j++) {
352 535 : pos_start=name.find("_",pos_start+1);
353 : }
354 : std::size_t pos_end=pos_start;
355 1503 : for(unsigned j=0; j<dim; j++) {
356 986 : pos_end=name.find("_",pos_end+1);
357 : }
358 517 : pos_start++; //do not include heading "_"
359 517 : return name.substr(pos_start,pos_end-pos_start);
360 : };
361 12 : unsigned index_j=ncv_;
362 : unsigned sizeSkip=1;
363 26 : for(int l=pntrToECVsClass_.size()-1; l>=0; l--) {
364 14 : const unsigned dim_l=pntrToECVsClass_[l]->getNumberOfArguments();
365 14 : index_j-=dim_l;
366 14 : std::vector<std::string> lambdas_l(1);
367 14 : lambdas_l[0]=getLambdaName(deltaF_name_[0],index_j,dim_l);
368 515 : for(unsigned i=sizeSkip; i<deltaF_size_; i+=sizeSkip) {
369 503 : std::string tmp_lambda=getLambdaName(deltaF_name_[i],index_j,dim_l);
370 503 : if(tmp_lambda==lambdas_l[0]) {
371 : break;
372 : }
373 501 : lambdas_l.push_back(tmp_lambda);
374 : }
375 14 : pntrToECVsClass_[l]->initECVs_restart(lambdas_l);
376 14 : sizeSkip*=lambdas_l.size();
377 14 : }
378 12 : plumed_massert(sizeSkip==deltaF_size_,"RESTART - this should not happen");
379 12 : init_linkECVs(); //link ECVs and initializes index_k_
380 12 : log.printf(" ->%4u DeltaFs in total\n",deltaF_size_);
381 12 : obs_steps_=0; //avoid initializing again
382 12 : if(stateRestart) {
383 3 : if(NumParallel_>1) {
384 2 : const unsigned start=(deltaF_size_/NumParallel_)*rank_+std::min(rank_,deltaF_size_%NumParallel_);
385 : unsigned iter=0;
386 34 : for(unsigned i=start; i<start+deltaF_.size(); i++) {
387 32 : deltaF_[iter++]=all_deltaF_[i];
388 : }
389 : }
390 : } else { //read each step
391 9 : counter_=1;
392 : unsigned count_lines=0;
393 9 : ifile.allowIgnoredFields(); //this allows for multiple restart, but without checking for consistency between them!
394 : double time;
395 54 : while(ifile.scanField("time",time)) { //only number of lines and last line is important
396 : unsigned restart_stride;
397 18 : ifile.scanField("print_stride",restart_stride);
398 18 : ifile.scanField("rct",rct_);
399 18 : if(NumParallel_==1) {
400 996 : for(unsigned i=0; i<deltaF_size_; i++) {
401 978 : ifile.scanField(deltaF_name_[i],deltaF_[i]);
402 : }
403 : } else {
404 0 : const unsigned start=(deltaF_size_/NumParallel_)*rank_+std::min(rank_,deltaF_size_%NumParallel_);
405 : unsigned iter=0;
406 0 : for(unsigned i=start; i<start+deltaF_.size(); i++) {
407 0 : ifile.scanField(deltaF_name_[i],deltaF_[iter++]);
408 : }
409 : }
410 18 : ifile.scanField();
411 18 : if(count_lines>0) {
412 9 : counter_+=restart_stride;
413 : }
414 18 : count_lines++;
415 : }
416 9 : counter_*=NumWalkers_;
417 9 : log.printf(" successfully read %u lines, up to t=%g\n",count_lines,time);
418 : }
419 12 : ifile.reset(false);
420 12 : ifile.close();
421 : } else { //same behaviour as METAD
422 0 : plumed_merror("RESTART requested, but file '"+restartFileName+"' was not found!\n Set RESTART=NO or provide a restart file");
423 : }
424 12 : if(NumWalkers_>1) { //make sure that all walkers are doing the same thing
425 2 : std::vector<unsigned long long> all_counter(NumWalkers_);
426 2 : if(comm.Get_rank()==0) {
427 2 : multi_sim_comm.Allgather(counter_,all_counter);
428 : }
429 2 : comm.Bcast(all_counter,0);
430 : bool same_number_of_steps=true;
431 4 : for(unsigned w=1; w<NumWalkers_; w++)
432 2 : if(all_counter[0]!=all_counter[w]) {
433 : same_number_of_steps=false;
434 : }
435 2 : plumed_massert(same_number_of_steps,"RESTART - not all walkers are reading the same file!");
436 : }
437 36 : } else if(restartFileName.length()>0) {
438 0 : log.printf(" +++ WARNING +++ the provided STATE_RFILE will be ignored, since RESTART was not requested\n");
439 : }
440 :
441 : //sync all walkers to avoid opening files before reding is over (see also METAD)
442 36 : comm.Barrier();
443 36 : if(comm.Get_rank()==0 && walkers_mpi) {
444 6 : multi_sim_comm.Barrier();
445 : }
446 :
447 : //setup DeltaFs file
448 36 : deltaFsOfile_.link(*this);
449 36 : if(NumWalkers_>1) {
450 6 : if(walker_rank_>0) {
451 : deltaFsFileName="/dev/null"; //only first walker writes on file
452 : }
453 12 : deltaFsOfile_.enforceSuffix("");
454 : }
455 36 : deltaFsOfile_.open(deltaFsFileName);
456 36 : if(fmt.length()>0) {
457 72 : deltaFsOfile_.fmtField(" "+fmt);
458 : }
459 : deltaFsOfile_.setHeavyFlush(); //do I need it?
460 36 : deltaFsOfile_.addConstantField("print_stride");
461 36 : deltaFsOfile_.printField("print_stride",print_stride_);
462 :
463 : //open file for storing state
464 36 : if(wStateStride_!=0) {
465 8 : stateOfile_.link(*this);
466 8 : if(NumWalkers_>1) {
467 0 : if(walker_rank_>0) {
468 : stateFileName="/dev/null"; //only first walker writes on file
469 : }
470 0 : stateOfile_.enforceSuffix("");
471 : }
472 8 : stateOfile_.open(stateFileName);
473 8 : if(fmt.length()>0) {
474 16 : stateOfile_.fmtField(" "+fmt);
475 : }
476 : }
477 :
478 : //add output components
479 36 : if(calc_work_) {
480 16 : addComponent("work");
481 16 : componentIsNotPeriodic("work");
482 : }
483 :
484 : //printing some info
485 36 : log.printf(" updating the bias with PACE = %u\n",stride_);
486 36 : log.printf(" initial unbiased OBSERVATION_STEPS = %u (in units of PACE)\n",obs_steps_);
487 36 : if(wStateStride_>0) {
488 7 : log.printf(" state checkpoints are written on file %s every %d MD steps\n",stateFileName.c_str(),wStateStride_);
489 : }
490 36 : if(wStateStride_==-1) {
491 1 : log.printf(" state checkpoints are written on file '%s' only on CPT events (or never if MD code does define them!)\n",stateFileName.c_str());
492 : }
493 36 : if(walkers_mpi) {
494 6 : log.printf(" -- WALKERS_MPI: if multiple replicas are present, they will share the same bias via MPI\n");
495 : }
496 36 : if(NumWalkers_>1) {
497 6 : log.printf(" using multiple walkers\n");
498 6 : log.printf(" number of walkers: %u\n",NumWalkers_);
499 6 : log.printf(" walker rank: %u\n",walker_rank_);
500 : }
501 36 : int mw_warning=0;
502 36 : if(!walkers_mpi && comm.Get_rank()==0 && multi_sim_comm.Get_size()>(int)NumWalkers_) {
503 0 : mw_warning=1;
504 : }
505 36 : comm.Bcast(mw_warning,0);
506 36 : if(mw_warning) { //log.printf messes up with comm, so never use it without Bcast!
507 0 : log.printf(" +++ WARNING +++ multiple replicas will NOT communicate unless the flag WALKERS_MPI is used\n");
508 : }
509 36 : if(NumParallel_>1) {
510 4 : log.printf(" using multiple MPI processes per simulation: %u\n",NumParallel_);
511 : }
512 36 : if(NumOMP_>1) {
513 31 : log.printf(" using multiple OpenMP threads per simulation: %u\n",NumOMP_);
514 : }
515 36 : if(serial) {
516 5 : log.printf(" -- SERIAL: no loop parallelization, despite %d MPI processes and %u OpenMP threads available\n",comm.Get_size(),OpenMP::getNumThreads());
517 : }
518 36 : log.printf(" Bibliography: ");
519 36 : log<<plumed.cite("M. Invernizzi, P.M. Piaggi, and M. Parrinello, Phys. Rev. X 10, 041034 (2020)");
520 36 : log.printf("\n");
521 36 : }
522 :
523 1794 : void OPESexpanded::calculate() {
524 1794 : if(deltaF_size_==0) { //no bias before initialization
525 409 : return;
526 : }
527 :
528 : //get diffMax, to avoid over/underflow
529 1385 : double diffMax=-std::numeric_limits<double>::max();
530 1385 : #pragma omp parallel num_threads(NumOMP_)
531 : {
532 : #pragma omp for reduction(max:diffMax)
533 : for(unsigned i=0; i<deltaF_.size(); i++) {
534 : diff_[i]=(-getExpansion(i)+deltaF_[i]/kbt_);
535 : if(diff_[i]>diffMax) {
536 : diffMax=diff_[i];
537 : }
538 : }
539 : }
540 1385 : if(NumParallel_>1) {
541 162 : comm.Max(diffMax);
542 : }
543 :
544 : //calculate the bias and the forces
545 1385 : double sum=0;
546 1385 : std::vector<double> der_sum_cv(ncv_,0);
547 1385 : if(NumOMP_==1) {
548 2730 : for(unsigned i=0; i<deltaF_.size(); i++) {
549 2520 : double add_i=std::exp(diff_[i]-diffMax);
550 2520 : sum+=add_i;
551 : //set derivatives
552 6960 : for(unsigned j=0; j<ncv_; j++) {
553 4440 : der_sum_cv[j]-=derECVs_[j][index_k_[i][j]]*add_i;
554 : }
555 : }
556 : } else {
557 1175 : #pragma omp parallel num_threads(NumOMP_)
558 : {
559 : std::vector<double> omp_der_sum_cv(ncv_,0);
560 : #pragma omp for reduction(+:sum) nowait
561 : for(unsigned i=0; i<deltaF_.size(); i++) {
562 : double add_i=std::exp(diff_[i]-diffMax);
563 : sum+=add_i;
564 : //set derivatives
565 : for(unsigned j=0; j<ncv_; j++) {
566 : omp_der_sum_cv[j]-=derECVs_[j][index_k_[i][j]]*add_i;
567 : }
568 : }
569 : #pragma omp critical
570 : for(unsigned j=0; j<ncv_; j++) {
571 : der_sum_cv[j]+=omp_der_sum_cv[j];
572 : }
573 : }
574 : }
575 1385 : if(NumParallel_>1) {
576 : //each MPI process has part of the full deltaF_ vector, so must Sum
577 162 : comm.Sum(sum);
578 162 : comm.Sum(der_sum_cv);
579 : }
580 :
581 : //set bias and forces
582 1385 : const double bias=-kbt_*(diffMax+std::log(sum/deltaF_size_));
583 1385 : setBias(bias);
584 3823 : for(unsigned j=0; j<ncv_; j++) {
585 2438 : setOutputForce(j,kbt_*der_sum_cv[j]/sum);
586 : }
587 : }
588 :
589 1794 : void OPESexpanded::update() {
590 1794 : if(isFirstStep_) { //skip very first step, as in METAD
591 36 : isFirstStep_=false;
592 36 : if(obs_steps_!=1) { //if obs_steps_==1 go on with initialization
593 : return;
594 : }
595 : }
596 1762 : if(getStep()%stride_==0) {
597 887 : if(obs_steps_>0) {
598 583 : for(unsigned j=0; j<ncv_; j++) {
599 384 : obs_cvs_[counter_*ncv_+j]=getArgument(j);
600 : }
601 199 : counter_++;
602 199 : if(counter_==obs_steps_) {
603 24 : log.printf("\nAction %s\n",getName().c_str());
604 24 : init_fromObs();
605 24 : log.printf("Finished initialization\n\n");
606 24 : counter_=NumWalkers_; //all preliminary observations count 1
607 24 : obs_steps_=0; //no more observation
608 : }
609 199 : return;
610 : }
611 :
612 : //update averages
613 688 : const double current_bias=getOutputQuantity(0); //the first value is always the bias
614 688 : if(NumWalkers_==1) {
615 580 : updateDeltaF(current_bias);
616 : } else {
617 108 : std::vector<double> cvs(ncv_);
618 324 : for(unsigned j=0; j<ncv_; j++) {
619 216 : cvs[j]=getArgument(j);
620 : }
621 108 : std::vector<double> all_bias(NumWalkers_);
622 108 : std::vector<double> all_cvs(NumWalkers_*ncv_);
623 108 : if(comm.Get_rank()==0) {
624 108 : multi_sim_comm.Allgather(current_bias,all_bias);
625 108 : multi_sim_comm.Allgather(cvs,all_cvs);
626 : }
627 108 : comm.Bcast(all_bias,0);
628 108 : comm.Bcast(all_cvs,0);
629 324 : for(unsigned w=0; w<NumWalkers_; w++) {
630 : //calculate ECVs
631 216 : unsigned index_wj=w*ncv_;
632 548 : for(unsigned k=0; k<pntrToECVsClass_.size(); k++) {
633 332 : pntrToECVsClass_[k]->calculateECVs(&all_cvs[index_wj]);
634 332 : index_wj+=pntrToECVsClass_[k]->getNumberOfArguments();
635 : }
636 216 : updateDeltaF(all_bias[w]);
637 : }
638 108 : if(calc_work_) { //the ECVs must be those of this walker
639 : unsigned index_j=0;
640 174 : for(unsigned k=0; k<pntrToECVsClass_.size(); k++) {
641 116 : pntrToECVsClass_[k]->calculateECVs(&cvs[index_j]);
642 116 : index_j+=pntrToECVsClass_[k]->getNumberOfArguments();
643 : }
644 : }
645 : }
646 :
647 : //write DeltaFs to file
648 688 : if((counter_/NumWalkers_-1)%print_stride_==0) {
649 52 : printDeltaF();
650 : }
651 :
652 : //calculate work if requested
653 688 : if(calc_work_) {
654 : //some copy and paste from calculate()
655 : //get diffMax, to avoid over/underflow
656 138 : double diffMax=-std::numeric_limits<double>::max();
657 138 : #pragma omp parallel num_threads(NumOMP_)
658 : {
659 : #pragma omp for reduction(max:diffMax)
660 : for(unsigned i=0; i<deltaF_.size(); i++) {
661 : diff_[i]=(-getExpansion(i)+deltaF_[i]/kbt_);
662 : if(diff_[i]>diffMax) {
663 : diffMax=diff_[i];
664 : }
665 : }
666 : }
667 138 : if(NumParallel_>1) {
668 50 : comm.Max(diffMax);
669 : }
670 : //calculate the bias
671 138 : double sum=0;
672 138 : #pragma omp parallel num_threads(NumOMP_)
673 : {
674 : #pragma omp for reduction(+:sum) nowait
675 : for(unsigned i=0; i<deltaF_.size(); i++) {
676 : sum+=std::exp(diff_[i]-diffMax);
677 : }
678 : }
679 138 : if(NumParallel_>1) {
680 50 : comm.Sum(sum);
681 : }
682 138 : const double new_bias=-kbt_*(diffMax+std::log(sum/deltaF_size_));
683 : //accumulate work
684 138 : work_+=new_bias-current_bias;
685 276 : getPntrToComponent("work")->set(work_);
686 : }
687 : }
688 :
689 : //dump state if requested, but not before initialization
690 1563 : if(obs_steps_==0 && ((wStateStride_>0 && getStep()%wStateStride_==0) || (wStateStride_==-1 && getCPT()))) {
691 31 : dumpStateToFile();
692 : }
693 : }
694 :
695 36 : void OPESexpanded::init_pntrToECVsClass() {
696 36 : std::vector<opes::ExpansionCVs*> all_pntrToECVsClass=plumed.getActionSet().select<opes::ExpansionCVs*>();
697 36 : plumed_massert(all_pntrToECVsClass.size()>0,"no Expansion CVs found");
698 81 : for(unsigned j=0; j<ncv_; j++) {
699 90 : std::string error_notECV("all the ARGs of "+getName()+" must be Expansion Collective Variables (ECV)");
700 45 : const unsigned dot_pos=getPntrToArgument(j)->getName().find(".");
701 45 : plumed_massert(dot_pos<getPntrToArgument(j)->getName().size(),error_notECV+", thus contain a dot in the name");
702 45 : unsigned foundECV_l=all_pntrToECVsClass.size();
703 55 : for(unsigned l=0; l<all_pntrToECVsClass.size(); l++) {
704 55 : if(getPntrToArgument(j)->getName().substr(0,dot_pos)==all_pntrToECVsClass[l]->getLabel()) {
705 : foundECV_l=l;
706 45 : pntrToECVsClass_.push_back(all_pntrToECVsClass[l]);
707 45 : std::string missing_arg="some ECV component is missing from ARG";
708 45 : plumed_massert(j+all_pntrToECVsClass[l]->getNumberOfArguments()<=getNumberOfArguments(),missing_arg);
709 110 : for(unsigned h=0; h<all_pntrToECVsClass[l]->getNumberOfArguments(); h++) {
710 65 : std::string argName=getPntrToArgument(j+h)->getName();
711 65 : std::string expectedECVname=all_pntrToECVsClass[l]->getComponentsVector()[h];
712 65 : plumed_massert(argName==expectedECVname,missing_arg+", or is in wrong order: given ARG="+argName+" expected ARG="+expectedECVname);
713 : }
714 45 : j+=all_pntrToECVsClass[l]->getNumberOfArguments()-1;
715 : break;
716 : }
717 : }
718 45 : plumed_massert(foundECV_l<all_pntrToECVsClass.size(),error_notECV);
719 : }
720 81 : for(unsigned l=0; l<pntrToECVsClass_.size(); l++)
721 54 : for(unsigned ll=l+1; ll<pntrToECVsClass_.size(); ll++) {
722 9 : plumed_massert(pntrToECVsClass_[l]->getLabel()!=pntrToECVsClass_[ll]->getLabel(),"cannot use same ECV twice");
723 : }
724 36 : }
725 :
726 36 : void OPESexpanded::init_linkECVs() {
727 : //TODO It should be possible to make all of this more straightforward (and probably also faster):
728 : // - get rid of index_k_, making it trivial for each ECV
729 : // - store the ECVs_ and derECVs_ vectors here as a contiguous vector, and use pointers in the ECV classes
730 : // Some caveats:
731 : // - ECVmultiThermalBaric has a nontrivial index_k_ to avoid duplicates. use duplicates instead
732 : // - can the ECVs be MPI parallel or it's too complicated?
733 36 : plumed_massert(deltaF_size_>0,"must set deltaF_size_ before calling init_linkECVs()");
734 36 : if(NumParallel_==1) {
735 32 : deltaF_.resize(deltaF_size_);
736 : } else {
737 4 : const unsigned extra=(rank_<(deltaF_size_%NumParallel_)?1:0);
738 4 : deltaF_.resize(deltaF_size_/NumParallel_+extra);
739 : //these are used when printing deltaF_ to file
740 4 : all_deltaF_.resize(deltaF_size_);
741 4 : all_size_.resize(NumParallel_,deltaF_size_/NumParallel_);
742 4 : disp_.resize(NumParallel_);
743 8 : for(unsigned r=0; r<NumParallel_-1; r++) {
744 4 : if(r<deltaF_size_%NumParallel_) {
745 0 : all_size_[r]++;
746 : }
747 4 : disp_[r+1]=disp_[r]+all_size_[r];
748 : }
749 : }
750 36 : diff_.resize(deltaF_.size());
751 36 : ECVs_.resize(ncv_);
752 36 : derECVs_.resize(ncv_);
753 36 : index_k_.resize(deltaF_.size(),std::vector<unsigned>(ncv_));
754 : unsigned index_j=0;
755 36 : unsigned sizeSkip=deltaF_size_;
756 81 : for(unsigned l=0; l<pntrToECVsClass_.size(); l++) {
757 45 : std::vector< std::vector<unsigned> > l_index_k(pntrToECVsClass_[l]->getIndex_k());
758 45 : plumed_massert(deltaF_size_%l_index_k.size()==0,"buggy ECV: mismatch between getTotNumECVs() and getIndex_k().size()");
759 45 : plumed_massert(l_index_k[0].size()==pntrToECVsClass_[l]->getNumberOfArguments(),"buggy ECV: mismatch between number of ARG and underlying CVs");
760 45 : sizeSkip/=l_index_k.size();
761 110 : for(unsigned h=0; h<pntrToECVsClass_[l]->getNumberOfArguments(); h++) {
762 65 : ECVs_[index_j+h]=pntrToECVsClass_[l]->getPntrToECVs(h);
763 65 : derECVs_[index_j+h]=pntrToECVsClass_[l]->getPntrToDerECVs(h);
764 65 : if(NumParallel_==1) {
765 44811 : for(unsigned i=0; i<deltaF_size_; i++) {
766 44754 : index_k_[i][index_j+h]=l_index_k[(i/sizeSkip)%l_index_k.size()][h];
767 : }
768 : } else {
769 8 : const unsigned start=(deltaF_size_/NumParallel_)*rank_+std::min(rank_,deltaF_size_%NumParallel_);
770 : unsigned iter=0;
771 76 : for(unsigned i=start; i<start+deltaF_.size(); i++) {
772 68 : index_k_[iter++][index_j+h]=l_index_k[(i/sizeSkip)%l_index_k.size()][h];
773 : }
774 : }
775 : }
776 45 : index_j+=pntrToECVsClass_[l]->getNumberOfArguments();
777 45 : }
778 36 : plumed_massert(sizeSkip==1,"this should not happen!");
779 36 : }
780 :
781 24 : void OPESexpanded::init_fromObs() { //This could probably be faster and/or require less memory...
782 : //in case of multiple walkers gather all the statistics
783 24 : if(NumWalkers_>1) {
784 4 : std::vector<double> all_obs_cv(ncv_*obs_steps_*NumWalkers_);
785 4 : if(comm.Get_rank()==0) {
786 4 : multi_sim_comm.Allgather(obs_cvs_,all_obs_cv);
787 : }
788 4 : comm.Bcast(all_obs_cv,0);
789 4 : obs_cvs_=all_obs_cv; //could this lead to memory issues?
790 4 : obs_steps_*=NumWalkers_;
791 : }
792 :
793 : //initialize ECVs from observations
794 : unsigned index_j=0;
795 24 : deltaF_size_=1;
796 55 : for(unsigned l=0; l<pntrToECVsClass_.size(); l++) {
797 31 : pntrToECVsClass_[l]->initECVs_observ(obs_cvs_,ncv_,index_j);
798 31 : deltaF_size_*=pntrToECVsClass_[l]->getTotNumECVs(); //ECVs from different exansions will be combined
799 31 : index_j+=pntrToECVsClass_[l]->getNumberOfArguments();
800 : }
801 24 : plumed_massert(index_j==getNumberOfArguments(),"mismatch between number of linked CVs and number of ARG");
802 : //link ECVs and initialize index_k_, mapping each deltaF to a single ECVs set
803 24 : init_linkECVs();
804 :
805 : //initialize deltaF_ from obs
806 : //for the first point, t=0, the ECVs are calculated by initECVs_observ, setting also any initial guess
807 12221 : for(unsigned i=0; i<deltaF_.size(); i++)
808 55995 : for(unsigned j=0; j<ncv_; j++) {
809 43798 : deltaF_[i]+=kbt_*ECVs_[j][index_k_[i][j]];
810 : }
811 239 : for(unsigned t=1; t<obs_steps_; t++) { //starts from t=1
812 : index_j=0;
813 533 : for(unsigned l=0; l<pntrToECVsClass_.size(); l++) {
814 318 : pntrToECVsClass_[l]->calculateECVs(&obs_cvs_[t*ncv_+index_j]);
815 318 : index_j+=pntrToECVsClass_[l]->getNumberOfArguments();
816 : }
817 101915 : for(unsigned i=0; i<deltaF_.size(); i++) {
818 101700 : const double diff_i=(-getExpansion(i)+deltaF_[i]/kbt_-std::log(t));
819 101700 : if(diff_i>0) { //save exp from overflow
820 15849 : deltaF_[i]-=kbt_*(diff_i+std::log1p(std::exp(-diff_i))+std::log1p(-1./(1.+t)));
821 : } else {
822 85851 : deltaF_[i]-=kbt_*(std::log1p(std::exp(diff_i))+std::log1p(-1./(1.+t)));
823 : }
824 : }
825 : }
826 24 : obs_cvs_.clear();
827 :
828 : //set deltaF_name_
829 24 : deltaF_name_.resize(deltaF_size_,"DeltaF");
830 24 : unsigned sizeSkip=deltaF_size_;
831 55 : for(unsigned l=0; l<pntrToECVsClass_.size(); l++) {
832 31 : std::vector<std::string> lambdas_l=pntrToECVsClass_[l]->getLambdas();
833 31 : plumed_massert(lambdas_l.size()==pntrToECVsClass_[l]->getTotNumECVs(),"buggy ECV: mismatch between getTotNumECVs() and getLambdas().size()");
834 31 : sizeSkip/=lambdas_l.size();
835 22146 : for(unsigned i=0; i<deltaF_size_; i++) {
836 44230 : deltaF_name_[i]+="_"+lambdas_l[(i/sizeSkip)%lambdas_l.size()];
837 : }
838 31 : }
839 :
840 : //print initialization to file
841 24 : log.printf(" ->%4u DeltaFs in total\n",deltaF_size_);
842 24 : printDeltaF();
843 24 : }
844 :
845 76 : void OPESexpanded::printDeltaF() {
846 76 : deltaFsOfile_.printField("time",getTime());
847 76 : deltaFsOfile_.printField("rct",rct_);
848 76 : if(NumParallel_==1) {
849 23652 : for(unsigned i=0; i<deltaF_.size(); i++) {
850 23582 : deltaFsOfile_.printField(deltaF_name_[i],deltaF_[i]);
851 : }
852 : } else {
853 6 : comm.Allgatherv(deltaF_,all_deltaF_,&all_size_[0],&disp_[0]); //can we avoid using this big vector?
854 78 : for(unsigned i=0; i<deltaF_size_; i++) {
855 72 : deltaFsOfile_.printField(deltaF_name_[i],all_deltaF_[i]);
856 : }
857 : }
858 76 : deltaFsOfile_.printField();
859 76 : }
860 :
861 31 : void OPESexpanded::dumpStateToFile() {
862 : //only the first walker writes the state, the others would write on /dev/null,
863 : //and their fields would pile up since they never rewind
864 31 : if(walker_rank_>0) {
865 : return;
866 : }
867 : //rewrite header or rewind file
868 31 : if(storeOldStates_) {
869 3 : stateOfile_.clearFields();
870 : } else {
871 28 : stateOfile_.rewind();
872 : }
873 : //define fields
874 31 : stateOfile_.addConstantField("time");
875 31 : stateOfile_.addConstantField("counter");
876 31 : stateOfile_.addConstantField("rct");
877 : //print
878 31 : stateOfile_.printField("time",getTime());
879 31 : stateOfile_.printField("counter",counter_);
880 31 : stateOfile_.printField("rct",rct_);
881 31 : if(NumParallel_>1) {
882 20 : comm.Allgatherv(deltaF_,all_deltaF_,&all_size_[0],&disp_[0]); //can we avoid using this big vector?
883 : }
884 299 : for(unsigned i=0; i<deltaF_size_; i++) {
885 : std::size_t pos_start=7; //skip "DeltaF_"
886 804 : for(unsigned j=0; j<ncv_; j++) {
887 : plumed_dbg_massert(pos_start>6,"not enought _ in deltaF_name_"+std::to_string(i-1)+" string?");
888 536 : const std::size_t pos_end=deltaF_name_[i].find("_",pos_start);
889 536 : stateOfile_.printField(getPntrToArgument(j)->getName()," "+deltaF_name_[i].substr(pos_start,pos_end-pos_start));
890 536 : pos_start=pos_end+1;
891 : }
892 268 : if(NumParallel_==1) {
893 456 : stateOfile_.printField("DeltaF",deltaF_[i]);
894 : } else {
895 80 : stateOfile_.printField("DeltaF",all_deltaF_[i]);
896 : }
897 268 : stateOfile_.printField();
898 : }
899 : //make sure file is written even if small
900 31 : if(!storeOldStates_) {
901 28 : stateOfile_.flush();
902 : }
903 : }
904 :
905 796 : void OPESexpanded::updateDeltaF(double bias) {
906 : plumed_dbg_massert(counter_>0,"deltaF_ must be initialized");
907 796 : counter_++;
908 796 : const double arg=(bias-rct_)/kbt_-std::log(counter_-1.);
909 : double increment;
910 796 : if(arg>0) { //save exp from overflow
911 32 : increment=kbt_*(arg+std::log1p(std::exp(-arg)));
912 : } else {
913 764 : increment=kbt_*(std::log1p(std::exp(arg)));
914 : }
915 796 : #pragma omp parallel num_threads(NumOMP_)
916 : {
917 : #pragma omp for
918 : for(unsigned i=0; i<deltaF_.size(); i++) {
919 : const double diff_i=(-getExpansion(i)+(bias-rct_+deltaF_[i])/kbt_-std::log(counter_-1.));
920 : if(diff_i>0) { //save exp from overflow
921 : deltaF_[i]+=increment-kbt_*(diff_i+std::log1p(std::exp(-diff_i)));
922 : } else {
923 : deltaF_[i]+=increment-kbt_*std::log1p(std::exp(diff_i));
924 : }
925 : }
926 : }
927 796 : rct_+=increment+kbt_*std::log1p(-1./counter_);
928 796 : }
929 :
930 637117 : double OPESexpanded::getExpansion(unsigned i) const {
931 : double expansion=0;
932 2957138 : for(unsigned j=0; j<ncv_; j++) {
933 2320021 : expansion+=ECVs_[j][index_k_[i][j]]; //the index_k could be trivially guessed for most ECVs, but unfourtunately not all
934 : }
935 637117 : return expansion;
936 : }
937 :
938 : }
939 : }
|
