Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2021 The VES code team
3 : (see the PEOPLE-VES file at the root of this folder for a list of names)
4 :
5 : See http://www.ves-code.org for more information.
6 :
7 : This file is part of VES code module.
8 :
9 : The VES code module 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 : The VES code module 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 the VES code module. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 :
23 : #include "bias/Bias.h"
24 : #include "core/PlumedMain.h"
25 : #include "core/ActionRegister.h"
26 : #include "core/Atoms.h"
27 : #include "tools/Communicator.h"
28 : #include "tools/Grid.h"
29 : #include "tools/File.h"
30 : //#include <algorithm> //std::fill
31 :
32 : namespace PLMD {
33 : namespace ves {
34 :
35 : //+PLUMEDOC VES_BIAS VES_DELTA_F
36 : /*
37 : Implementation of VES\f$\Delta F\f$ method \cite Invernizzi2019vesdeltaf (step two only).
38 :
39 : \warning
40 : Notice that this is a stand-alone bias Action, it does not need any of the other VES module components
41 :
42 : First you should create some estimate of the local free energy basins of your system,
43 : using e.g. multiple \ref METAD short runs, and combining them with the \ref sum_hills utility.
44 : Once you have them, you can use this bias Action to perform the VES optimization part of the method.
45 :
46 : These \f$N+1\f$ local basins are used to model the global free energy.
47 : In particular, given the conditional probabilities \f$P(\mathbf{s}|i)\propto e^{-\beta F_i(\mathbf{s})}\f$
48 : and the probabilities of being in a given basin \f$P_i\f$, we can write:
49 : \f[
50 : e^{-\beta F(\mathbf{s})}\propto P(\mathbf{s})=\sum_{i=0}^N P(\mathbf{s}|i)P_i \, .
51 : \f]
52 : We use this free energy model and the chosen bias factor \f$\gamma\f$ to build the bias potential:
53 : \f$V(\mathbf{s})=-(1-1/\gamma)F(\mathbf{s})\f$.
54 : Or, more explicitly:
55 : \f[
56 : V(\mathbf{s})=(1-1/\gamma)\frac{1}{\beta}\log\left[e^{-\beta F_0(\mathbf{s})}
57 : +\sum_{i=1}^{N} e^{-\beta F_i(\mathbf{s})} e^{-\beta \alpha_i}\right] \, ,
58 : \f]
59 : where the parameters \f$\boldsymbol{\alpha}\f$ are the \f$N\f$ free energy differences (see below) from the \f$F_0\f$ basin.
60 :
61 : By default the \f$F_i(\mathbf{s})\f$ are shifted so that \f$\min[F_i(\mathbf{s})]=0\f$ for all \f$i=\{0,...,N\}\f$.
62 : In this case the optimization parameters \f$\alpha_i\f$ are the difference in height between the minima of the basins.
63 : Using the keyword `NORMALIZE`, you can also decide to normalize the local free energies so that
64 : \f$\int d\mathbf{s}\, e^{-\beta F_i(\mathbf{s})}=1\f$.
65 : In this case the parameters will represent not the difference in height (which depends on the chosen CVs),
66 : but the actual free energy difference, \f$\alpha_i=\Delta F_i\f$.
67 :
68 : However, as discussed in Ref. \cite Invernizzi2019vesdeltaf, a better estimate of \f$\Delta F_i\f$ should be obtained through the reweighting procedure.
69 :
70 : \par Examples
71 :
72 : The following performs the optimization of the free energy difference between two metastable basins:
73 :
74 : \plumedfile
75 : cv: DISTANCE ATOMS=1,2
76 : ves: VES_DELTA_F ...
77 : ARG=cv
78 : TEMP=300
79 : FILE_F0=fesA.data
80 : FILE_F1=fesB.data
81 : BIASFACTOR=10.0
82 : M_STEP=0.1
83 : AV_STRIDE=500
84 : PRINT_STRIDE=100
85 : ...
86 : PRINT FMT=%g STRIDE=500 FILE=Colvar.data ARG=cv,ves.bias,ves.rct
87 : \endplumedfile
88 :
89 : The local FES files can be obtained as described in Sec. 4.2 of Ref. \cite Invernizzi2019vesdeltaf, i.e. for example:
90 : - run 4 independent metad runs, all starting from basin A, and kill them as soon as they make the first transition (see e.g. \ref COMMITTOR)
91 : - \verbatim cat HILLS* > all_HILLS \endverbatim
92 : - \verbatim plumed sum_hills --hills all_HILLS --outfile all_fesA.dat --mintozero --min 0 --max 1 --bin 100 \endverbatim
93 : - \verbatim awk -v n_rep=4 '{if($1!="#!" && $1!="") {for(i=1+(NF-1)/2; i<=NF; i++) $i/=n_rep;} print $0}' all_fesA.dat > fesA.data \endverbatim
94 :
95 : The header of both FES files must be identical, and should be similar to the following:
96 :
97 : \auxfile{fesA.data}
98 : #! FIELDS cv file.free der_cv
99 : #! SET min_cv 0
100 : #! SET max_cv 1
101 : #! SET nbins_cv 100
102 : #! SET periodic_cv false
103 : 0 0 0
104 : \endauxfile
105 : \auxfile{fesB.data}
106 : #! FIELDS cv file.free der_cv
107 : #! SET min_cv 0
108 : #! SET max_cv 1
109 : #! SET nbins_cv 100
110 : #! SET periodic_cv false
111 : 0 0 0
112 : \endauxfile
113 :
114 : */
115 : //+ENDPLUMEDOC
116 :
117 : class VesDeltaF : public bias::Bias {
118 :
119 : private:
120 : double beta_;
121 : unsigned NumParallel_;
122 : unsigned rank_;
123 : unsigned NumWalkers_;
124 : bool isFirstStep_;
125 : bool afterCalculate_;
126 :
127 : //prob
128 : double tot_prob_;
129 : std::vector<double> prob_;
130 : std::vector< std::vector<double> > der_prob_;
131 :
132 : //local basins
133 : std::vector< std::unique_ptr<Grid> > grid_p_; //pointers because of GridBase::create
134 : std::vector<double> norm_;
135 :
136 : //optimizer-related stuff
137 : long unsigned mean_counter_;
138 : unsigned mean_weight_tau_;
139 : unsigned alpha_size_;
140 : unsigned sym_alpha_size_;
141 : std::vector<double> mean_alpha_;
142 : std::vector<double> inst_alpha_;
143 : std::vector<double> past_increment2_;
144 : double minimization_step_;
145 : bool damping_off_;
146 : //'tg' -> 'target distribution'
147 : double inv_gamma_;
148 : unsigned tg_counter_;
149 : unsigned tg_stride_;
150 : std::vector<double> tg_dV_dAlpha_;
151 : std::vector<double> tg_d2V_dAlpha2_;
152 : //'av' -> 'ensemble average'
153 : unsigned av_counter_;
154 : unsigned av_stride_;
155 : std::vector<double> av_dV_dAlpha_;
156 : std::vector<double> av_dV_dAlpha_prod_;
157 : std::vector<double> av_d2V_dAlpha2_;
158 : //printing
159 : unsigned print_stride_;
160 : OFile alphaOfile_;
161 : //other
162 : std::vector<double> exp_alpha_;
163 : std::vector<double> prev_exp_alpha_;
164 : double work_;
165 :
166 : //functions
167 : void update_alpha();
168 : void update_tg_and_rct();
169 : inline unsigned get_index(const unsigned, const unsigned) const;
170 :
171 : public:
172 : explicit VesDeltaF(const ActionOptions&);
173 : void calculate() override;
174 : void update() override;
175 : static void registerKeywords(Keywords& keys);
176 : };
177 :
178 10427 : PLUMED_REGISTER_ACTION(VesDeltaF,"VES_DELTA_F")
179 :
180 5 : void VesDeltaF::registerKeywords(Keywords& keys) {
181 5 : Bias::registerKeywords(keys);
182 5 : keys.use("ARG");
183 10 : keys.add("optional","TEMP","temperature is compulsory, but it can be sometimes fetched from the MD engine");
184 : //local free energies
185 10 : keys.add("numbered","FILE_F","names of files containing local free energies and derivatives. "
186 : "The first one, FILE_F0, is used as reference for all the free energy differences.");
187 10 : keys.reset_style("FILE_F","compulsory");
188 10 : keys.addFlag("NORMALIZE",false,"normalize all local free energies so that alpha will be (approx) \\f$\\Delta F\\f$");
189 10 : keys.addFlag("NO_MINTOZERO",false,"leave local free energies as provided, without shifting them to zero min");
190 : //target distribution
191 10 : keys.add("compulsory","BIASFACTOR","0","the \\f$\\gamma\\f$ bias factor used for well-tempered target \\f$p(\\mathbf{s})\\f$."
192 : " Set to 0 for non-tempered flat target");
193 10 : keys.add("optional","TG_STRIDE","( default=1 ) number of AV_STRIDE between updates"
194 : " of target \\f$p(\\mathbf{s})\\f$ and reweighing factor \\f$c(t)\\f$");
195 : //optimization
196 10 : keys.add("compulsory","M_STEP","1.0","the \\f$\\mu\\f$ step used for the \\f$\\Omega\\f$ functional minimization");
197 10 : keys.add("compulsory","AV_STRIDE","500","number of simulation steps between alpha updates");
198 10 : keys.add("optional","TAU_MEAN","exponentially decaying average for alpha (rescaled using AV_STRIDE)."
199 : " Should be used only in very specific cases");
200 10 : keys.add("optional","INITIAL_ALPHA","( default=0 ) an initial guess for the bias potential parameter alpha");
201 10 : keys.addFlag("DAMPING_OFF",false,"do not use an AdaGrad-like term to rescale M_STEP");
202 : //output parameters file
203 10 : keys.add("compulsory","ALPHA_FILE","ALPHA","file name for output minimization parameters");
204 10 : keys.add("optional","PRINT_STRIDE","( default=10 ) stride for printing to ALPHA_FILE");
205 10 : keys.add("optional","FMT","specify format for ALPHA_FILE");
206 : //debug flags
207 10 : keys.addFlag("SERIAL",false,"perform the calculation in serial even if multiple tasks are available");
208 10 : keys.addFlag("MULTIPLE_WALKERS",false,"use multiple walkers connected via MPI for the optimization");
209 5 : keys.use("RESTART");
210 :
211 : //output components
212 5 : componentsAreNotOptional(keys);
213 10 : keys.addOutputComponent("rct","default","the reweighting factor \\f$c(t)\\f$");
214 10 : keys.addOutputComponent("work","default","the work done by the bias in one AV_STRIDE");
215 5 : }
216 :
217 4 : VesDeltaF::VesDeltaF(const ActionOptions&ao)
218 : : PLUMED_BIAS_INIT(ao)
219 4 : , isFirstStep_(true)
220 4 : , afterCalculate_(false)
221 4 : , mean_counter_(0)
222 4 : , av_counter_(0)
223 4 : , work_(0)
224 : {
225 : //set beta
226 4 : const double Kb=plumed.getAtoms().getKBoltzmann();
227 4 : double temp=0;
228 4 : parse("TEMP",temp);
229 4 : double KbT=Kb*temp;
230 4 : if(KbT==0)
231 : {
232 0 : KbT=plumed.getAtoms().getKbT();
233 0 : plumed_massert(KbT>0,"your MD engine does not pass the temperature to plumed, you must specify it using TEMP");
234 : }
235 4 : beta_=1.0/KbT;
236 :
237 : //initialize probability grids using local free energies
238 : bool spline=true;
239 : bool sparsegrid=false;
240 4 : std::string funcl="file.free"; //typical name given by sum_hills
241 :
242 : std::vector<std::string> fes_names;
243 8 : for(unsigned n=0;; n++)//NB: here we start from FILE_F0 not from FILE_F1
244 : {
245 : std::string filename;
246 24 : if(!parseNumbered("FILE_F",n,filename))
247 : break;
248 8 : fes_names.push_back(filename);
249 8 : IFile gridfile;
250 8 : gridfile.open(filename);
251 8 : auto g=GridBase::create(funcl,getArguments(),gridfile,sparsegrid,spline,true);
252 : // we assume this cannot be sparse. in case we want it to be sparse, some of the methods
253 : // that are available only in Grid should be ported to GridBase
254 8 : auto gg=dynamic_cast<Grid*>(g.get());
255 : // if this throws, g is deleted
256 8 : plumed_assert(gg);
257 : // release ownership in order to transfer it to emplaced pointer
258 : g.release();
259 8 : grid_p_.emplace_back(gg);
260 16 : }
261 4 : plumed_massert(grid_p_.size()>1,"at least 2 basins must be defined, starting from FILE_F0");
262 4 : alpha_size_=grid_p_.size()-1;
263 4 : sym_alpha_size_=alpha_size_*(alpha_size_+1)/2; //useful for symmetric matrix [alpha_size_]x[alpha_size_]
264 : //check for consistency with first local free energy
265 8 : for(unsigned n=1; n<grid_p_.size(); n++)
266 : {
267 8 : std::string error_tag="FILE_F"+std::to_string(n)+" '"+fes_names[n]+"' not compatible with reference one, FILE_F0";
268 4 : plumed_massert(grid_p_[n]->getSize()==grid_p_[0]->getSize(),error_tag);
269 4 : plumed_massert(grid_p_[n]->getMin()==grid_p_[0]->getMin(),error_tag);
270 4 : plumed_massert(grid_p_[n]->getMax()==grid_p_[0]->getMax(),error_tag);
271 4 : plumed_massert(grid_p_[n]->getBinVolume()==grid_p_[0]->getBinVolume(),error_tag);
272 : }
273 :
274 4 : bool no_mintozero=false;
275 4 : parseFlag("NO_MINTOZERO",no_mintozero);
276 4 : if(!no_mintozero)
277 : {
278 6 : for(unsigned n=0; n<grid_p_.size(); n++)
279 4 : grid_p_[n]->setMinToZero();
280 : }
281 4 : bool normalize=false;
282 4 : parseFlag("NORMALIZE",normalize);
283 4 : norm_.resize(grid_p_.size(),0);
284 4 : std::vector<double> c_norm(grid_p_.size());
285 : //convert the FESs to probability distributions
286 : //NB: the spline interpolation will be done on the probability distributions, not on the given FESs
287 : const unsigned ncv=getNumberOfArguments(); //just for ease
288 12 : for(unsigned n=0; n<grid_p_.size(); n++)
289 : {
290 808 : for(Grid::index_t t=0; t<grid_p_[n]->getSize(); t++)
291 : {
292 800 : std::vector<double> der(ncv);
293 800 : const double val=std::exp(-beta_*grid_p_[n]->getValueAndDerivatives(t,der));
294 1600 : for(unsigned s=0; s<ncv; s++)
295 800 : der[s]*=-beta_*val;
296 800 : grid_p_[n]->setValueAndDerivatives(t,val,der);
297 800 : norm_[n]+=val;
298 : }
299 8 : c_norm[n]=1./beta_*std::log(norm_[n]);
300 8 : if(normalize)
301 : {
302 4 : grid_p_[n]->scaleAllValuesAndDerivatives(1./norm_[n]);
303 4 : norm_[n]=1;
304 : }
305 : }
306 :
307 : //get target
308 4 : double biasfactor=0;
309 4 : parse("BIASFACTOR",biasfactor);
310 4 : plumed_massert(biasfactor==0 || biasfactor>1,"BIASFACTOR must be zero (for uniform target) or greater than one");
311 4 : if(biasfactor==0)
312 2 : inv_gamma_=0;
313 : else
314 2 : inv_gamma_=1./biasfactor;
315 4 : tg_counter_=0;
316 4 : tg_stride_=1;
317 4 : parse("TG_STRIDE",tg_stride_);
318 4 : tg_dV_dAlpha_.resize(alpha_size_,0);
319 4 : tg_d2V_dAlpha2_.resize(sym_alpha_size_,0);
320 :
321 : //setup optimization stuff
322 4 : minimization_step_=1;
323 4 : parse("M_STEP",minimization_step_);
324 :
325 4 : av_stride_=500;
326 4 : parse("AV_STRIDE",av_stride_);
327 4 : av_dV_dAlpha_.resize(alpha_size_,0);
328 4 : av_dV_dAlpha_prod_.resize(sym_alpha_size_,0);
329 4 : av_d2V_dAlpha2_.resize(sym_alpha_size_,0);
330 :
331 4 : mean_weight_tau_=0;
332 4 : parse("TAU_MEAN",mean_weight_tau_);
333 4 : if(mean_weight_tau_!=1) //set it to 1 for basic SGD
334 : {
335 4 : plumed_massert((mean_weight_tau_==0 || mean_weight_tau_>av_stride_),"TAU_MEAN is rescaled with AV_STRIDE, so it has to be greater");
336 4 : mean_weight_tau_/=av_stride_; //this way you can look at the number of simulation steps to choose TAU_MEAN
337 : }
338 :
339 8 : parseVector("INITIAL_ALPHA",mean_alpha_);
340 4 : if(mean_alpha_.size()>0)
341 : {
342 2 : plumed_massert(mean_alpha_.size()==alpha_size_,"provide one INITIAL_ALPHA for each basin beyond the first one");
343 : }
344 : else
345 2 : mean_alpha_.resize(alpha_size_,0);
346 4 : inst_alpha_=mean_alpha_;
347 4 : exp_alpha_.resize(alpha_size_);
348 8 : for(unsigned i=0; i<alpha_size_; i++)
349 4 : exp_alpha_[i]=std::exp(-beta_*mean_alpha_[i]);
350 4 : prev_exp_alpha_=exp_alpha_;
351 :
352 4 : damping_off_=false;
353 4 : parseFlag("DAMPING_OFF",damping_off_);
354 4 : if(damping_off_)
355 2 : past_increment2_.resize(alpha_size_,1);
356 : else
357 2 : past_increment2_.resize(alpha_size_,0);
358 :
359 : //file printing options
360 4 : std::string alphaFileName("ALPHA");
361 4 : parse("ALPHA_FILE",alphaFileName);
362 4 : print_stride_=10;
363 8 : parse("PRINT_STRIDE",print_stride_);
364 : std::string fmt;
365 4 : parse("FMT",fmt);
366 :
367 : //other flags, mainly for debugging
368 4 : NumParallel_=comm.Get_size();
369 4 : rank_=comm.Get_rank();
370 4 : bool serial=false;
371 4 : parseFlag("SERIAL",serial);
372 4 : if(serial)
373 : {
374 2 : log.printf(" -- SERIAL: running without loop parallelization\n");
375 2 : NumParallel_=1;
376 2 : rank_=0;
377 : }
378 :
379 4 : bool multiple_walkers=false;
380 4 : parseFlag("MULTIPLE_WALKERS",multiple_walkers);
381 4 : if(!multiple_walkers)
382 2 : NumWalkers_=1;
383 : else
384 : {
385 2 : if(comm.Get_rank()==0)//multi_sim_comm works well on first rank only
386 2 : NumWalkers_=multi_sim_comm.Get_size();
387 2 : if(comm.Get_size()>1) //if each walker has more than one processor update them all
388 0 : comm.Bcast(NumWalkers_,0);
389 : }
390 :
391 4 : checkRead();
392 :
393 : //restart if needed
394 4 : if(getRestart())
395 : {
396 2 : IFile ifile;
397 2 : ifile.link(*this);
398 2 : if(NumWalkers_>1)
399 4 : ifile.enforceSuffix("");
400 2 : if(ifile.FileExist(alphaFileName))
401 : {
402 2 : log.printf(" Restarting from: %s\n",alphaFileName.c_str());
403 2 : log.printf(" all options (also PRINT_STRIDE) must be consistent!\n");
404 2 : log.printf(" any INITIAL_ALPHA will be overwritten\n");
405 2 : ifile.open(alphaFileName);
406 : double time;
407 2 : std::vector<double> damping(alpha_size_);
408 20 : while(ifile.scanField("time",time)) //room for improvements: only last line is important
409 : {
410 16 : for(unsigned i=0; i<alpha_size_; i++)
411 : {
412 8 : const std::string index(std::to_string(i+1));
413 8 : prev_exp_alpha_[i]=std::exp(-beta_*mean_alpha_[i]);
414 16 : ifile.scanField("alpha_"+index,mean_alpha_[i]);
415 16 : ifile.scanField("auxiliary_"+index,inst_alpha_[i]);
416 16 : ifile.scanField("damping_"+index,damping[i]);
417 : }
418 8 : ifile.scanField();
419 8 : mean_counter_+=print_stride_;
420 : }
421 4 : for(unsigned i=0; i<alpha_size_; i++)
422 : {
423 2 : exp_alpha_[i]=std::exp(-beta_*mean_alpha_[i]);
424 2 : past_increment2_[i]=damping[i]*damping[i];
425 : }
426 : //sync all walkers and treads. Not sure is mandatory but is no harm
427 2 : comm.Barrier();
428 2 : if(comm.Get_rank()==0)
429 2 : multi_sim_comm.Barrier();
430 : }
431 : else
432 0 : log.printf(" -- WARNING: restart requested, but no '%s' file found!\n",alphaFileName.c_str());
433 2 : }
434 :
435 : //setup output file with Alpha values
436 4 : alphaOfile_.link(*this);
437 4 : if(NumWalkers_>1)
438 : {
439 2 : if(comm.Get_rank()==0 && multi_sim_comm.Get_rank()>0)
440 : alphaFileName="/dev/null"; //only first walker writes on file
441 4 : alphaOfile_.enforceSuffix("");
442 : }
443 4 : alphaOfile_.open(alphaFileName);
444 4 : if(fmt.length()>0)
445 8 : alphaOfile_.fmtField(" "+fmt);
446 :
447 : //add other output components
448 8 : addComponent("rct"); componentIsNotPeriodic("rct");
449 8 : addComponent("work"); componentIsNotPeriodic("work");
450 :
451 : //print some info
452 4 : log.printf(" Temperature T: %g\n",1./(Kb*beta_));
453 4 : log.printf(" Beta (1/Kb*T): %g\n",beta_);
454 4 : log.printf(" Local free energy basins files and normalization constants:\n");
455 12 : for(unsigned n=0; n<grid_p_.size(); n++)
456 8 : log.printf(" F_%d filename: %s c_%d=%g\n",n,fes_names[n].c_str(),n,c_norm[n]);
457 4 : if(no_mintozero)
458 2 : log.printf(" -- NO_MINTOZERO: local free energies are not shifted to be zero at minimum\n");
459 4 : if(normalize)
460 2 : log.printf(" -- NORMALIZE: F_n+=c_n, alpha=DeltaF\n");
461 4 : log.printf(" Using target distribution with 1/gamma = %g\n",inv_gamma_);
462 4 : log.printf(" and updated with stride %d\n",tg_stride_);
463 4 : log.printf(" Step for the minimization algorithm: %g\n",minimization_step_);
464 4 : log.printf(" Stride for the ensemble average: %d\n",av_stride_);
465 4 : if(mean_weight_tau_>1)
466 2 : log.printf(" Exponentially decaying average with weight=tau/av_stride=%d\n",mean_weight_tau_);
467 4 : if(mean_weight_tau_==1)
468 0 : log.printf(" +++ WARNING +++ setting TAU_MEAN=1 is equivalent to use simple SGD, without mean alpha nor hessian contribution\n");
469 4 : log.printf(" Initial guess for alpha:\n");
470 8 : for(unsigned i=0; i<alpha_size_; i++)
471 4 : log.printf(" alpha_%d = %g\n",i+1,mean_alpha_[i]);
472 4 : if(damping_off_)
473 2 : log.printf(" -- DAMPING_OFF: the minimization step will NOT become smaller as the simulation goes on\n");
474 4 : log.printf(" Printing on file %s with stride %d\n",alphaFileName.c_str(),print_stride_);
475 4 : if(serial)
476 2 : log.printf(" -- SERIAL: running without loop parallelization\n");
477 4 : if(NumParallel_>1)
478 2 : log.printf(" Using multiple threads per simulation: %d\n",NumParallel_);
479 4 : if(multiple_walkers)
480 : {
481 2 : log.printf(" -- MULTIPLE_WALKERS: multiple simulations will combine statistics for the optimization\n");
482 2 : if(NumWalkers_>1)
483 : {
484 2 : log.printf(" number of walkers: %d\n",NumWalkers_);
485 2 : log.printf(" walker rank: %d\n",multi_sim_comm.Get_rank()); //only comm.Get_rank()=0 prints, so this is fine
486 : }
487 : else
488 0 : log.printf(" +++ WARNING +++ only one replica found: are you sure you are running MPI-connected simulations?\n");
489 : }
490 4 : log.printf(" Bibliography ");
491 8 : log<<plumed.cite("Invernizzi and Parrinello, J. Chem. Theory Comput. 15, 2187-2194 (2019)");
492 8 : log<<plumed.cite("Valsson and Parrinello, Phys. Rev. Lett. 113, 090601 (2014)");
493 4 : if(inv_gamma_>0)
494 6 : log<<plumed.cite("Valsson and Parrinello, J. Chem. Theory Comput. 11, 1996-2002 (2015)");
495 :
496 : //last initializations
497 4 : prob_.resize(grid_p_.size());
498 4 : der_prob_.resize(grid_p_.size(),std::vector<double>(getNumberOfArguments()));
499 4 : update_tg_and_rct();
500 8 : }
501 :
502 804 : void VesDeltaF::calculate()
503 : {
504 : //get CVs
505 804 : const unsigned ncv=getNumberOfArguments(); //just for ease
506 804 : std::vector<double> cv(ncv);
507 1608 : for(unsigned s=0; s<ncv; s++)
508 804 : cv[s]=getArgument(s);
509 : //get probabilities for each basin, and total one
510 2412 : for(unsigned n=0; n<grid_p_.size(); n++)
511 1608 : prob_[n]=grid_p_[n]->getValueAndDerivatives(cv,der_prob_[n]);
512 804 : tot_prob_=prob_[0];
513 1608 : for(unsigned i=0; i<alpha_size_; i++)
514 804 : tot_prob_+=prob_[i+1]*exp_alpha_[i];
515 :
516 : //update bias and forces: V=-(1-inv_gamma_)*fes
517 804 : setBias((1-inv_gamma_)/beta_*std::log(tot_prob_));
518 1608 : for(unsigned s=0; s<ncv; s++)
519 : {
520 804 : double dProb_dCV_s=der_prob_[0][s];
521 1608 : for(unsigned i=0; i<alpha_size_; i++)
522 804 : dProb_dCV_s+=der_prob_[i+1][s]*exp_alpha_[i];
523 804 : setOutputForce(s,-(1-inv_gamma_)/beta_/tot_prob_*dProb_dCV_s);
524 : }
525 804 : afterCalculate_=true;
526 804 : }
527 :
528 804 : void VesDeltaF::update()
529 : {
530 : //skip first step to sync getTime() and av_counter_, as in METAD
531 804 : if(isFirstStep_)
532 : {
533 4 : isFirstStep_=false;
534 4 : return;
535 : }
536 800 : plumed_massert(afterCalculate_,"VesDeltaF::update() must be called after VesDeltaF::calculate() to work properly");
537 800 : afterCalculate_=false;
538 :
539 : //calculate derivatives for ensemble averages
540 800 : std::vector<double> dV_dAlpha(alpha_size_);
541 800 : std::vector<double> d2V_dAlpha2(sym_alpha_size_);
542 1600 : for(unsigned i=0; i<alpha_size_; i++)
543 800 : dV_dAlpha[i]=-(1-inv_gamma_)/tot_prob_*prob_[i+1]*exp_alpha_[i];
544 1600 : for(unsigned i=0; i<alpha_size_; i++)
545 : {
546 800 : d2V_dAlpha2[get_index(i,i)]=-beta_*dV_dAlpha[i];
547 1600 : for(unsigned j=i; j<alpha_size_; j++)
548 800 : d2V_dAlpha2[get_index(i,j)]-=beta_/(1-inv_gamma_)*dV_dAlpha[i]*dV_dAlpha[j];
549 : }
550 : //update ensemble averages
551 800 : av_counter_++;
552 1600 : for(unsigned i=0; i<alpha_size_; i++)
553 : {
554 800 : av_dV_dAlpha_[i]+=(dV_dAlpha[i]-av_dV_dAlpha_[i])/av_counter_;
555 1600 : for(unsigned j=i; j<alpha_size_; j++)
556 : {
557 800 : const unsigned ij=get_index(i,j);
558 800 : av_dV_dAlpha_prod_[ij]+=(dV_dAlpha[i]*dV_dAlpha[j]-av_dV_dAlpha_prod_[ij])/av_counter_;
559 800 : av_d2V_dAlpha2_[ij]+=(d2V_dAlpha2[ij]-av_d2V_dAlpha2_[ij])/av_counter_;
560 : }
561 : }
562 : //update work
563 800 : double prev_tot_prob=prob_[0];
564 1600 : for(unsigned i=0; i<alpha_size_; i++)
565 800 : prev_tot_prob+=prob_[i+1]*prev_exp_alpha_[i];
566 800 : work_+=(1-inv_gamma_)/beta_*std::log(tot_prob_/prev_tot_prob);
567 :
568 : //update coefficients
569 800 : if(av_counter_==av_stride_)
570 : {
571 16 : update_alpha();
572 16 : tg_counter_++;
573 16 : if(tg_counter_==tg_stride_)
574 : {
575 12 : update_tg_and_rct();
576 12 : tg_counter_=0;
577 : }
578 : //reset the ensemble averages
579 16 : av_counter_=0;
580 : std::fill(av_dV_dAlpha_.begin(),av_dV_dAlpha_.end(),0);
581 : std::fill(av_dV_dAlpha_prod_.begin(),av_dV_dAlpha_prod_.end(),0);
582 : std::fill(av_d2V_dAlpha2_.begin(),av_d2V_dAlpha2_.end(),0);
583 : }
584 : }
585 :
586 16 : void VesDeltaF::update_tg_and_rct()
587 : {
588 : //calculate target averages
589 16 : double Z_0=norm_[0];
590 32 : for(unsigned i=0; i<alpha_size_; i++)
591 16 : Z_0+=norm_[i+1]*exp_alpha_[i];
592 16 : double Z_tg=0;
593 : std::fill(tg_dV_dAlpha_.begin(),tg_dV_dAlpha_.end(),0);
594 : std::fill(tg_d2V_dAlpha2_.begin(),tg_d2V_dAlpha2_.end(),0);
595 1116 : for(Grid::index_t t=rank_; t<grid_p_[0]->getSize(); t+=NumParallel_)
596 : { //TODO can we recycle some code?
597 1100 : std::vector<double> prob(grid_p_.size());
598 3300 : for(unsigned n=0; n<grid_p_.size(); n++)
599 2200 : prob[n]=grid_p_[n]->getValue(t);
600 1100 : double tot_prob=prob[0];
601 2200 : for(unsigned i=0; i<alpha_size_; i++)
602 1100 : tot_prob+=prob[i+1]*exp_alpha_[i];
603 1100 : std::vector<double> dV_dAlpha(alpha_size_);
604 1100 : std::vector<double> d2V_dAlpha2(sym_alpha_size_);
605 2200 : for(unsigned i=0; i<alpha_size_; i++)
606 1100 : dV_dAlpha[i]=-(1-inv_gamma_)/tot_prob*prob[i+1]*exp_alpha_[i];
607 2200 : for(unsigned i=0; i<alpha_size_; i++)
608 : {
609 1100 : d2V_dAlpha2[get_index(i,i)]=-beta_*dV_dAlpha[i];
610 2200 : for(unsigned j=i; j<alpha_size_; j++)
611 1100 : d2V_dAlpha2[get_index(i,j)]-=beta_/(1-inv_gamma_)*dV_dAlpha[i]*dV_dAlpha[j];
612 : }
613 1100 : const double unnorm_tg_p=std::pow(tot_prob,inv_gamma_);
614 1100 : Z_tg+=unnorm_tg_p;
615 2200 : for(unsigned i=0; i<alpha_size_; i++)
616 1100 : tg_dV_dAlpha_[i]+=unnorm_tg_p*dV_dAlpha[i];
617 2200 : for(unsigned ij=0; ij<sym_alpha_size_; ij++)
618 1100 : tg_d2V_dAlpha2_[ij]+=unnorm_tg_p*d2V_dAlpha2[ij];
619 : }
620 16 : if(NumParallel_>1)
621 : {
622 10 : comm.Sum(Z_tg);
623 10 : comm.Sum(tg_dV_dAlpha_);
624 10 : comm.Sum(tg_d2V_dAlpha2_);
625 : }
626 32 : for(unsigned i=0; i<alpha_size_; i++)
627 16 : tg_dV_dAlpha_[i]/=Z_tg;
628 32 : for(unsigned ij=0; ij<sym_alpha_size_; ij++)
629 16 : tg_d2V_dAlpha2_[ij]/=Z_tg;
630 16 : getPntrToComponent("rct")->set(-1./beta_*std::log(Z_tg/Z_0)); //Z_tg is the best available estimate of Z_V
631 16 : }
632 :
633 16 : void VesDeltaF::update_alpha()
634 : {
635 : //combining the averages of multiple walkers
636 16 : if(NumWalkers_>1)
637 : {
638 8 : if(comm.Get_rank()==0) //sum only once: in the first rank of each walker
639 : {
640 8 : multi_sim_comm.Sum(av_dV_dAlpha_);
641 8 : multi_sim_comm.Sum(av_dV_dAlpha_prod_);
642 8 : multi_sim_comm.Sum(av_d2V_dAlpha2_);
643 16 : for(unsigned i=0; i<alpha_size_; i++)
644 8 : av_dV_dAlpha_[i]/=NumWalkers_;
645 16 : for(unsigned ij=0; ij<sym_alpha_size_; ij++)
646 : {
647 8 : av_dV_dAlpha_prod_[ij]/=NumWalkers_;
648 8 : av_d2V_dAlpha2_[ij]/=NumWalkers_;
649 : }
650 : }
651 8 : if(comm.Get_size()>1)//if there are more ranks for each walker, everybody has to know
652 : {
653 0 : comm.Bcast(av_dV_dAlpha_,0);
654 0 : comm.Bcast(av_dV_dAlpha_prod_,0);
655 0 : comm.Bcast(av_d2V_dAlpha2_,0);
656 : }
657 : }
658 : //set work and reset it
659 16 : getPntrToComponent("work")->set(work_);
660 16 : work_=0;
661 :
662 : //build the gradient and the Hessian of the functional
663 16 : std::vector<double> grad_omega(alpha_size_);
664 16 : std::vector<double> hess_omega(sym_alpha_size_);
665 32 : for(unsigned i=0; i<alpha_size_; i++)
666 : {
667 16 : grad_omega[i]=tg_dV_dAlpha_[i]-av_dV_dAlpha_[i];
668 32 : for(unsigned j=i; j<alpha_size_; j++)
669 : {
670 16 : const unsigned ij=get_index(i,j);
671 16 : hess_omega[ij]=beta_*(av_dV_dAlpha_prod_[ij]-av_dV_dAlpha_[i]*av_dV_dAlpha_[j])+tg_d2V_dAlpha2_[ij]-av_d2V_dAlpha2_[ij];
672 : }
673 : }
674 : //calculate the increment and update alpha
675 16 : mean_counter_++;
676 : long unsigned mean_weight=mean_counter_;
677 16 : if(mean_weight_tau_>0 && mean_weight_tau_<mean_counter_)
678 : mean_weight=mean_weight_tau_;
679 16 : std::vector<double> damping(alpha_size_);
680 32 : for(unsigned i=0; i<alpha_size_; i++)
681 : {
682 16 : double increment_i=grad_omega[i];
683 32 : for(unsigned j=0; j<alpha_size_; j++)
684 16 : increment_i+=hess_omega[get_index(i,j)]*(inst_alpha_[j]-mean_alpha_[j]);
685 16 : if(!damping_off_)
686 8 : past_increment2_[i]+=increment_i*increment_i;
687 16 : damping[i]=std::sqrt(past_increment2_[i]);
688 16 : prev_exp_alpha_[i]=std::exp(-beta_*mean_alpha_[i]);
689 16 : inst_alpha_[i]-=minimization_step_/damping[i]*increment_i;
690 16 : mean_alpha_[i]+=(inst_alpha_[i]-mean_alpha_[i])/mean_weight;
691 16 : exp_alpha_[i]=std::exp(-beta_*mean_alpha_[i]);
692 : }
693 :
694 : //update the Alpha file
695 16 : if(mean_counter_%print_stride_==0)
696 : {
697 16 : alphaOfile_.printField("time",getTime());
698 32 : for(unsigned i=0; i<alpha_size_; i++)
699 : {
700 16 : const std::string index(std::to_string(i+1));
701 32 : alphaOfile_.printField("alpha_"+index,mean_alpha_[i]);
702 32 : alphaOfile_.printField("auxiliary_"+index,inst_alpha_[i]);
703 32 : alphaOfile_.printField("damping_"+index,damping[i]);
704 : }
705 16 : alphaOfile_.printField();
706 : }
707 16 : }
708 :
709 : //mapping of a [alpha_size_]x[alpha_size_] symmetric matrix into a vector of size sym_alpha_size_, useful for the communicator
710 4632 : inline unsigned VesDeltaF::get_index(const unsigned i, const unsigned j) const
711 : {
712 4632 : if(i<=j)
713 4632 : return j+i*(alpha_size_-1)-i*(i-1)/2;
714 : else
715 0 : return get_index(j,i);
716 : }
717 :
718 : }
719 : }
|
