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