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