Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2018 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "Bias.h"
23 : #include "ActionRegister.h"
24 : #include "core/ActionSet.h"
25 : #include "tools/Grid.h"
26 : #include "core/PlumedMain.h"
27 : #include "core/Atoms.h"
28 : #include "tools/Exception.h"
29 : #include "tools/Random.h"
30 : #include <string>
31 : #include <cstring>
32 : #include "tools/File.h"
33 : #include "time.h"
34 : #include <iostream>
35 : #include <limits>
36 :
37 : #define DP2CUTOFF 6.25
38 :
39 : using namespace std;
40 :
41 :
42 : namespace PLMD {
43 : namespace bias {
44 :
45 : //+PLUMEDOC BIAS PBMETAD
46 : /*
47 : Used to performed Parallel Bias MetaDynamics.
48 :
49 : This action activate Parallel Bias MetaDynamics (PBMetaD) \cite pbmetad, a version of MetaDynamics \cite metad in which
50 : multiple low-dimensional bias potentials are applied in parallel.
51 : In the current implementation, these have the form of mono-dimensional MetaDynamics bias
52 : potentials:
53 :
54 : \f[
55 : {V(s_1,t), ..., V(s_N,t)}
56 : \f]
57 :
58 : where:
59 :
60 : \f[
61 : V(s_i,t) = \sum_{ k \tau < t} W_i(k \tau)
62 : \exp\left(
63 : - \frac{(s_i-s_i^{(0)}(k \tau))^2}{2\sigma_i^2}
64 : \right).
65 : \f]
66 :
67 : To ensure the convergence of each mono-dimensional bias potential to the corresponding free energy,
68 : at each deposition step the Gaussian heights are multiplied by the so-called conditional term:
69 :
70 : \f[
71 : W_i(k \tau)=W_0 \frac{\exp\left(
72 : - \frac{V(s_i,k \tau)}{k_B T}
73 : \right)}{\sum_{i=1}^N
74 : \exp\left(
75 : - \frac{V(s_i,k \tau)}{k_B T}
76 : \right)}
77 : \f]
78 :
79 : where \f$W_0\f$ is the initial Gaussian height.
80 :
81 : The PBMetaD bias potential is defined by:
82 :
83 : \f[
84 : V_{PB}(\vec{s},t) = -k_B T \log{\sum_{i=1}^N
85 : \exp\left(
86 : - \frac{V(s_i,t)}{k_B T}
87 : \right)}.
88 : \f]
89 :
90 :
91 : Information on the Gaussian functions that build each bias potential are printed to
92 : multiple HILLS files, which
93 : are used both to restart the calculation and to reconstruct the mono-dimensional
94 : free energies as a function of the corresponding CVs.
95 : These can be reconstructed using the \ref sum_hills utility because the final bias is given
96 : by:
97 :
98 : \f[
99 : V(s_i) = -F(s_i)
100 : \f]
101 :
102 : Currently, only a subset of the \ref METAD options are available in PBMetaD.
103 :
104 : The bias potentials can be stored on a grid to increase performances of long PBMetaD simulations.
105 : You should
106 : provide either the number of bins for every collective variable (GRID_BIN) or
107 : the desired grid spacing (GRID_SPACING). In case you provide both PLUMED will use
108 : the most conservative choice (highest number of bins) for each dimension.
109 : In case you do not provide any information about bin size (neither GRID_BIN nor GRID_SPACING)
110 : and if Gaussian width is fixed PLUMED will use 1/5 of the Gaussian width as grid spacing.
111 : This default choice should be reasonable for most applications.
112 :
113 : Another option that is available is well-tempered metadynamics \cite Barducci:2008. In this
114 : variant of PBMetaD the heights of the Gaussian hills are rescaled at each step by the
115 : additional well-tempered metadynamics term.
116 : This ensures that each bias converges more smoothly. It should be noted that, in the case of well-tempered metadynamics, in
117 : the output printed the Gaussian height is re-scaled using the bias factor.
118 : Also notice that with well-tempered metadynamics the HILLS files do not contain the bias,
119 : but the negative of the free-energy estimate. This choice has the advantage that
120 : one can restart a simulation using a different value for the \f$\Delta T\f$. The applied bias will be scaled accordingly.
121 :
122 : With the keyword INTERVAL one changes the metadynamics algorithm setting the bias force equal to zero
123 : outside boundary \cite baftizadeh2012protein. If, for example, metadynamics is performed on a CV s and one is interested only
124 : to the free energy for s > sw, the history dependent potential is still updated according to the above
125 : equations but the metadynamics force is set to zero for s < sw. Notice that Gaussians are added also
126 : if s < sw, as the tails of these Gaussians influence VG in the relevant region s > sw. In this way, the
127 : force on the system in the region s > sw comes from both metadynamics and the force field, in the region
128 : s < sw only from the latter. This approach allows obtaining a history-dependent bias potential VG that
129 : fluctuates around a stable estimator, equal to the negative of the free energy far enough from the
130 : boundaries. Note that:
131 : - It works only for one-dimensional biases;
132 : - It works both with and without GRID;
133 : - The interval limit sw in a region where the free energy derivative is not large;
134 : - If in the region outside the limit sw the system has a free energy minimum, the INTERVAL keyword should
135 : be used together with a \ref UPPER_WALLS or \ref LOWER_WALLS at sw.
136 :
137 : Multiple walkers \cite multiplewalkers can also be used, in the MPI implementation only. See below the examples.
138 :
139 : \par Examples
140 :
141 : The following input is for PBMetaD calculation using as
142 : collective variables the distance between atoms 3 and 5
143 : and the distance between atoms 2 and 4. The value of the CVs and
144 : the PBMetaD bias potential are written to the COLVAR file every 100 steps.
145 : \verbatim
146 : DISTANCE ATOMS=3,5 LABEL=d1
147 : DISTANCE ATOMS=2,4 LABEL=d2
148 : PBMETAD ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3 PACE=500 LABEL=pb FILE=HILLS_d1,HILLS_d2
149 : PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
150 : \endverbatim
151 : (See also \ref DISTANCE and \ref PRINT).
152 :
153 : \par
154 : If you use well-tempered metadynamics, you should specify a single biasfactor and initial
155 : Gaussian height.
156 : \verbatim
157 : DISTANCE ATOMS=3,5 LABEL=d1
158 : DISTANCE ATOMS=2,4 LABEL=d2
159 : PBMETAD ...
160 : ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3
161 : PACE=500 BIASFACTOR=8 LABEL=pb
162 : FILE=HILLS_d1,HILLS_d2
163 : ... PBMETAD
164 : PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
165 : \endverbatim
166 :
167 : \par
168 : Only the MPI version of multiple-walkers is currently implemented.
169 : \verbatim
170 : DISTANCE ATOMS=3,5 LABEL=d1
171 : DISTANCE ATOMS=2,4 LABEL=d2
172 : PBMETAD ...
173 : ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3
174 : PACE=500 BIASFACTOR=8 LABEL=pb
175 : FILE=HILLS_d1,HILLS_d2
176 : MULTIPLE_WALKERS
177 : ... PBMETAD
178 : PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
179 : \endverbatim
180 :
181 : */
182 : //+ENDPLUMEDOC
183 :
184 : class PBMetaD : public Bias {
185 :
186 : private:
187 1068 : struct Gaussian {
188 : vector<double> center;
189 : vector<double> sigma;
190 : double height;
191 192 : Gaussian(const vector<double> & center,const vector<double> & sigma, double height):
192 192 : center(center),sigma(sigma),height(height) {}
193 : };
194 : vector<double> sigma0_;
195 : vector< vector<Gaussian> > hills_;
196 : vector<OFile*> hillsOfiles_;
197 : vector<OFile*> gridfiles_;
198 : vector<Grid*> BiasGrids_;
199 : bool grid_;
200 : double height0_;
201 : double biasf_;
202 : double kbt_;
203 : int stride_;
204 : int wgridstride_;
205 : bool welltemp_;
206 : bool multiple_w;
207 : unsigned nw_;
208 : unsigned mw_;
209 : vector<double> uppI_;
210 : vector<double> lowI_;
211 : vector<bool> doInt_;
212 : bool isFirstStep;
213 :
214 : void readGaussians(unsigned iarg, IFile*);
215 : bool readChunkOfGaussians(unsigned iarg, IFile *ifile, unsigned n);
216 : void writeGaussian(unsigned iarg, const Gaussian&, OFile*);
217 : void addGaussian(unsigned iarg, const Gaussian&);
218 : double getBiasAndDerivatives(unsigned iarg, const vector<double>&, double* der=NULL);
219 : double evaluateGaussian(unsigned iarg, const vector<double>&, const Gaussian&,double* der=NULL);
220 : vector<unsigned> getGaussianSupport(unsigned iarg, const Gaussian&);
221 : bool scanOneHill(unsigned iarg, IFile *ifile, vector<Value> &v, vector<double> ¢er, vector<double> &sigma, double &height);
222 : std::string fmt;
223 :
224 : public:
225 : explicit PBMetaD(const ActionOptions&);
226 : ~PBMetaD();
227 : void calculate();
228 : void update();
229 : static void registerKeywords(Keywords& keys);
230 : };
231 :
232 2529 : PLUMED_REGISTER_ACTION(PBMetaD,"PBMETAD")
233 :
234 7 : void PBMetaD::registerKeywords(Keywords& keys) {
235 7 : Bias::registerKeywords(keys);
236 7 : keys.use("ARG");
237 7 : keys.add("compulsory","SIGMA","the widths of the Gaussian hills");
238 7 : keys.add("compulsory","PACE","the frequency for hill addition, one for all biases");
239 7 : keys.add("optional","FILE","files in which the lists of added hills are stored");
240 7 : keys.add("optional","HEIGHT","the height of the Gaussian hills, one for all biases. Compulsory unless TAU, TEMP and BIASFACTOR are given");
241 7 : keys.add("optional","FMT","specify format for HILLS files (useful for decrease the number of digits in regtests)");
242 7 : keys.add("optional","BIASFACTOR","use well tempered metadynamics with this biasfactor, one for all biases. Please note you must also specify temp");
243 7 : keys.add("optional","TEMP","the system temperature - this is only needed if you are doing well-tempered metadynamics");
244 7 : keys.add("optional","TAU","in well tempered metadynamics, sets height to (kb*DeltaT*pace*timestep)/tau");
245 7 : keys.add("optional","GRID_RFILES", "read grid for the bias");
246 7 : keys.add("optional","GRID_WFILES", "dump grid for the bias");
247 7 : keys.add("optional","GRID_WSTRIDE", "frequency for dumping the grid");
248 7 : keys.add("optional","GRID_MIN","the lower bounds for the grid");
249 7 : keys.add("optional","GRID_MAX","the upper bounds for the grid");
250 7 : keys.add("optional","GRID_BIN","the number of bins for the grid");
251 7 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
252 7 : keys.addFlag("GRID_SPARSE",false,"use a sparse grid to store hills");
253 7 : keys.addFlag("GRID_NOSPLINE",false,"don't use spline interpolation with grids");
254 7 : keys.add("optional","INTERVAL_MIN","monodimensional lower limits, outside the limits the system will not feel the biasing force.");
255 7 : keys.add("optional","INTERVAL_MAX","monodimensional upper limits, outside the limits the system will not feel the biasing force.");
256 7 : keys.addFlag("MULTIPLE_WALKERS",false,"Switch on MPI version of multiple walkers");
257 7 : keys.use("RESTART");
258 7 : keys.use("UPDATE_FROM");
259 7 : keys.use("UPDATE_UNTIL");
260 7 : }
261 :
262 24 : PBMetaD::~PBMetaD() {
263 6 : for(unsigned i=0; i<BiasGrids_.size(); ++i) delete BiasGrids_[i];
264 18 : for(unsigned i=0; i<hillsOfiles_.size(); ++i) {
265 12 : hillsOfiles_[i]->close();
266 12 : delete hillsOfiles_[i];
267 : }
268 6 : if(wgridstride_ > 0) {
269 3 : for(unsigned i=0; i<gridfiles_.size(); ++i) {
270 2 : gridfiles_[i]->close();
271 2 : delete gridfiles_[i];
272 : }
273 : }
274 18 : }
275 :
276 6 : PBMetaD::PBMetaD(const ActionOptions& ao):
277 : PLUMED_BIAS_INIT(ao),
278 6 : grid_(false), height0_(std::numeric_limits<double>::max()),
279 : biasf_(1.0), kbt_(0.0), stride_(0), wgridstride_(0), welltemp_(false),
280 12 : multiple_w(false), isFirstStep(true)
281 : {
282 :
283 6 : parse("FMT",fmt);
284 :
285 : // parse the sigma
286 6 : parseVector("SIGMA",sigma0_);
287 6 : if( sigma0_.size()!=getNumberOfArguments() ) error("number of arguments does not match number of SIGMA parameters");
288 :
289 : // note: HEIGHT is not compulsory, since one could use the TAU keyword, see below
290 6 : parse("HEIGHT",height0_);
291 6 : parse("PACE",stride_);
292 6 : if(stride_<=0) error("frequency for hill addition is nonsensical");
293 :
294 6 : vector<string> hillsfname;
295 6 : parseVector("FILE",hillsfname);
296 6 : if(hillsfname.size()==0) {
297 0 : for(unsigned i=0; i<getNumberOfArguments(); i++) hillsfname.push_back("HILLS."+getPntrToArgument(i)->getName());
298 : }
299 6 : if( hillsfname.size()!=getNumberOfArguments() ) {
300 0 : error("number of FILE arguments does not match number of HILLS files");
301 : }
302 :
303 6 : parse("BIASFACTOR",biasf_);
304 6 : if( biasf_<1.0 ) error("well tempered bias factor is nonsensical");
305 6 : double temp=0.0;
306 6 : parse("TEMP",temp);
307 6 : if(temp>0.0) kbt_=plumed.getAtoms().getKBoltzmann()*temp;
308 0 : else kbt_=plumed.getAtoms().getKbT();
309 6 : if(biasf_>1.0) {
310 5 : if(kbt_==0.0) error("Unless the MD engine passes the temperature to plumed, with well-tempered metad you must specify it using TEMP");
311 5 : welltemp_=true;
312 : }
313 6 : double tau=0.0;
314 6 : parse("TAU",tau);
315 6 : if(tau==0.0) {
316 6 : if(height0_==std::numeric_limits<double>::max()) error("At least one between HEIGHT and TAU should be specified");
317 : // if tau is not set, we compute it here from the other input parameters
318 6 : if(welltemp_) tau=(kbt_*(biasf_-1.0))/height0_*getTimeStep()*stride_;
319 : } else {
320 0 : if(!welltemp_)error("TAU only makes sense in well-tempered metadynamics");
321 0 : if(height0_!=std::numeric_limits<double>::max()) error("At most one between HEIGHT and TAU should be specified");
322 0 : height0_=(kbt_*(biasf_-1.0))/tau*getTimeStep()*stride_;
323 : }
324 :
325 : // MPI version
326 6 : parseFlag("MULTIPLE_WALKERS",multiple_w);
327 :
328 : // Grid file
329 12 : vector<string> gridfilenames_;
330 6 : parseVector("GRID_WFILES",gridfilenames_);
331 6 : parse("GRID_WSTRIDE",wgridstride_);
332 6 : if (wgridstride_ == 0 && gridfilenames_.size() > 0) {
333 0 : error("frequency with which to output grid not specified use GRID_WSTRIDE");
334 : }
335 6 : if(gridfilenames_.size() > 0 && hillsfname.size() > 0 && gridfilenames_.size() != hillsfname.size())
336 0 : error("number of GRID_WFILES arguments does not match number of HILLS files");
337 :
338 : // Read grid
339 12 : vector<string> gridreadfilenames_;
340 6 : parseVector("GRID_RFILES",gridreadfilenames_);
341 :
342 : // Grid Stuff
343 12 : vector<std::string> gmin(getNumberOfArguments());
344 6 : parseVector("GRID_MIN",gmin);
345 6 : if(gmin.size()!=getNumberOfArguments() && gmin.size()!=0) error("not enough values for GRID_MIN");
346 12 : vector<std::string> gmax(getNumberOfArguments());
347 6 : parseVector("GRID_MAX",gmax);
348 6 : if(gmax.size()!=getNumberOfArguments() && gmax.size()!=0) error("not enough values for GRID_MAX");
349 12 : vector<unsigned> gbin(getNumberOfArguments());
350 12 : vector<double> gspacing;
351 6 : parseVector("GRID_BIN",gbin);
352 6 : if(gbin.size()!=getNumberOfArguments() && gbin.size()!=0) error("not enough values for GRID_BIN");
353 6 : parseVector("GRID_SPACING",gspacing);
354 6 : if(gspacing.size()!=getNumberOfArguments() && gspacing.size()!=0) error("not enough values for GRID_SPACING");
355 6 : if(gmin.size()!=gmax.size()) error("GRID_MAX and GRID_MIN should be either present or absent");
356 6 : if(gspacing.size()!=0 && gmin.size()==0) error("If GRID_SPACING is present also GRID_MIN should be present");
357 6 : if(gbin.size()!=0 && gmin.size()==0) error("If GRID_SPACING is present also GRID_MIN should be present");
358 6 : if(gmin.size()!=0) {
359 1 : if(gbin.size()==0 && gspacing.size()==0) {
360 1 : log<<" Binsize not specified, 1/10 of sigma will be be used\n";
361 1 : plumed_assert(sigma0_.size()==getNumberOfArguments());
362 1 : gspacing.resize(getNumberOfArguments());
363 1 : for(unsigned i=0; i<gspacing.size(); i++) gspacing[i]=0.1*sigma0_[i];
364 0 : } else if(gspacing.size()!=0 && gbin.size()==0) {
365 0 : log<<" The number of bins will be estimated from GRID_SPACING\n";
366 0 : } else if(gspacing.size()!=0 && gbin.size()!=0) {
367 0 : log<<" You specified both GRID_BIN and GRID_SPACING\n";
368 0 : log<<" The more conservative (highest) number of bins will be used for each variable\n";
369 : }
370 1 : if(gbin.size()==0) gbin.assign(getNumberOfArguments(),1);
371 3 : if(gspacing.size()!=0) for(unsigned i=0; i<getNumberOfArguments(); i++) {
372 : double a,b;
373 2 : Tools::convert(gmin[i],a);
374 2 : Tools::convert(gmax[i],b);
375 2 : unsigned n=((b-a)/gspacing[i])+1;
376 2 : if(gbin[i]<n) gbin[i]=n;
377 : }
378 : }
379 6 : bool sparsegrid=false;
380 6 : parseFlag("GRID_SPARSE",sparsegrid);
381 6 : bool nospline=false;
382 6 : parseFlag("GRID_NOSPLINE",nospline);
383 6 : bool spline=!nospline;
384 6 : if(gbin.size()>0) {grid_=true;}
385 6 : if(!grid_&&gridfilenames_.size() > 0) error("To write a grid you need first to define it!");
386 6 : if(!grid_&&gridreadfilenames_.size() > 0) error("To read a grid you need first to define it!");
387 :
388 6 : doInt_.resize(getNumberOfArguments(),false);
389 : // Interval keyword
390 6 : parseVector("INTERVAL_MIN",lowI_);
391 6 : parseVector("INTERVAL_MAX",uppI_);
392 : // various checks
393 6 : if(lowI_.size()!=uppI_.size()) error("both a lower and an upper limits must be provided with INTERVAL");
394 6 : if(lowI_.size()!=0 && lowI_.size()!=getNumberOfArguments()) error("check number of argument of INTERVAL");
395 6 : for(unsigned i=0; i<lowI_.size(); ++i) {
396 0 : if(uppI_[i]<lowI_[i]) error("The Upper limit must be greater than the Lower limit!");
397 0 : if(getPntrToArgument(i)->isPeriodic()) warning("INTERVAL is not used for periodic variables");
398 0 : else doInt_[i]=true;
399 : }
400 :
401 6 : checkRead();
402 :
403 6 : log.printf(" Gaussian width ");
404 6 : for(unsigned i=0; i<sigma0_.size(); ++i) log.printf(" %f",sigma0_[i]);
405 6 : log.printf(" Gaussian height %f\n",height0_);
406 6 : log.printf(" Gaussian deposition pace %d\n",stride_);
407 :
408 6 : log.printf(" Gaussian files ");
409 6 : for(unsigned i=0; i<hillsfname.size(); ++i) log.printf("%s ",hillsfname[i].c_str());
410 6 : log.printf("\n");
411 6 : if(welltemp_) {
412 5 : log.printf(" Well-Tempered Bias Factor %f\n",biasf_);
413 5 : log.printf(" Hills relaxation time (tau) %f\n",tau);
414 5 : log.printf(" KbT %f\n",kbt_);
415 : }
416 6 : if(multiple_w) log.printf(" Multiple walkers active using MPI communnication\n");
417 18 : for(unsigned i=0; i<doInt_.size(); i++) {
418 12 : if(doInt_[i]) log.printf(" Upper and Lower limits boundaries for the bias of CV %u are activated\n", i);
419 : }
420 6 : if(grid_) {
421 1 : log.printf(" Grid min");
422 1 : for(unsigned i=0; i<gmin.size(); ++i) log.printf(" %s",gmin[i].c_str() );
423 1 : log.printf("\n");
424 1 : log.printf(" Grid max");
425 1 : for(unsigned i=0; i<gmax.size(); ++i) log.printf(" %s",gmax[i].c_str() );
426 1 : log.printf("\n");
427 1 : log.printf(" Grid bin");
428 1 : for(unsigned i=0; i<gbin.size(); ++i) log.printf(" %u",gbin[i]);
429 1 : log.printf("\n");
430 1 : if(spline) {log.printf(" Grid uses spline interpolation\n");}
431 1 : if(sparsegrid) {log.printf(" Grid uses sparse grid\n");}
432 1 : if(wgridstride_>0) {
433 3 : for(unsigned i=0; i<gridfilenames_.size(); ++i) {
434 2 : log.printf(" Grid is written on file %s with stride %d\n",gridfilenames_[i].c_str(),wgridstride_);
435 : }
436 : }
437 1 : if(gridreadfilenames_.size()>0) {
438 0 : for(unsigned i=0; i<gridreadfilenames_.size(); ++i) {
439 0 : log.printf(" Reading bias from grid in file %s \n",gridreadfilenames_[i].c_str());
440 : }
441 : }
442 : }
443 :
444 : // initializing vector of hills
445 6 : hills_.resize(getNumberOfArguments());
446 :
447 : // restart from external grid
448 6 : bool restartedFromGrid=false;
449 :
450 : // initializing and checking grid
451 6 : if(grid_) {
452 : // check for mesh and sigma size
453 3 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
454 : double a,b;
455 2 : Tools::convert(gmin[i],a);
456 2 : Tools::convert(gmax[i],b);
457 2 : double mesh=(b-a)/((double)gbin[i]);
458 2 : if(mesh>0.5*sigma0_[i]) log<<" WARNING: Using a METAD with a Grid Spacing larger than half of the Gaussians width can produce artifacts\n";
459 : }
460 1 : std::string funcl=getLabel() + ".bias";
461 3 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
462 2 : std::vector<Value*> args(1);
463 2 : args[0] = getPntrToArgument(i);
464 4 : vector<std::string> gmin_t(1);
465 4 : vector<std::string> gmax_t(1);
466 4 : vector<unsigned> gbin_t(1);
467 2 : gmin_t[0] = gmin[i];
468 2 : gmax_t[0] = gmax[i];
469 2 : gbin_t[0] = gbin[i];
470 : Grid* BiasGrid_;
471 : // Read grid from file
472 2 : if(gridreadfilenames_.size()>0) {
473 0 : IFile gridfile;
474 0 : gridfile.link(*this);
475 0 : if(gridfile.FileExist(gridreadfilenames_[i])) {
476 0 : gridfile.open(gridreadfilenames_[i]);
477 : } else {
478 0 : error("The GRID file you want to read: " + gridreadfilenames_[i] + ", cannot be found!");
479 : }
480 0 : string funcl = getLabel() + ".bias";
481 0 : BiasGrid_ = Grid::create(funcl, args, gridfile, gmin_t, gmax_t, gbin_t, sparsegrid, spline, true);
482 0 : gridfile.close();
483 0 : if(BiasGrid_->getDimension() != args.size()) {
484 0 : error("mismatch between dimensionality of input grid and number of arguments");
485 : }
486 0 : if(getPntrToArgument(i)->isPeriodic() != BiasGrid_->getIsPeriodic()[0]) {
487 0 : error("periodicity mismatch between arguments and input bias");
488 : }
489 0 : log.printf(" Restarting from %s:",gridreadfilenames_[i].c_str());
490 0 : if(getRestart()) restartedFromGrid=true;
491 : } else {
492 2 : if(!sparsegrid) {BiasGrid_=new Grid(funcl,args,gmin_t,gmax_t,gbin_t,spline,true);}
493 0 : else {BiasGrid_=new SparseGrid(funcl,args,gmin_t,gmax_t,gbin_t,spline,true);}
494 2 : std::vector<std::string> actualmin=BiasGrid_->getMin();
495 4 : std::vector<std::string> actualmax=BiasGrid_->getMax();
496 2 : if(gmin_t[0]!=actualmin[0]) log<<" WARNING: GRID_MIN["<<i<<"] has been adjusted to "<<actualmin[0]<<" to fit periodicity\n";
497 4 : if(gmax_t[0]!=actualmax[0]) log<<" WARNING: GRID_MAX["<<i<<"] has been adjusted to "<<actualmax[0]<<" to fit periodicity\n";
498 : }
499 2 : BiasGrids_.push_back(BiasGrid_);
500 3 : }
501 : }
502 :
503 : // read Gaussians if restarting
504 6 : if (!restartedFromGrid) {
505 18 : for(unsigned i=0; i<hillsfname.size(); ++i) {
506 12 : IFile *ifile = new IFile();
507 12 : ifile->link(*this);
508 12 : if(ifile->FileExist(hillsfname[i])) {
509 0 : ifile->open(hillsfname[i]);
510 0 : if(getRestart()) {
511 0 : log.printf(" Restarting from %s:",hillsfname[i].c_str());
512 0 : readGaussians(i, ifile);
513 : }
514 0 : ifile->reset(false);
515 0 : ifile->close();
516 0 : delete ifile;
517 : }
518 : }
519 : }
520 :
521 6 : comm.Barrier();
522 : // this barrier is needed when using walkers_mpi
523 : // to be sure that all files have been read before
524 : // backing them up
525 : // it should not be used when walkers_mpi is false otherwise
526 : // it would introduce troubles when using replicas without METAD
527 : // (e.g. in bias exchange with a neutral replica)
528 : // see issue #168 on github
529 6 : if(multiple_w) {
530 4 : if(comm.Get_rank()==0) {
531 2 : multi_sim_comm.Barrier();
532 2 : nw_ = multi_sim_comm.Get_size();
533 2 : mw_ = multi_sim_comm.Get_rank();
534 : }
535 4 : comm.Bcast(nw_,0);
536 4 : comm.Bcast(mw_,0);
537 : }
538 :
539 : // open hills files for writing
540 18 : for(unsigned i=0; i<hillsfname.size(); ++i) {
541 12 : OFile *ofile = new OFile();
542 12 : ofile->link(*this);
543 12 : string hillsfname_tmp = hillsfname[i];
544 : // if multiple walkers, only rank 0 will write to file
545 12 : if(multiple_w) {
546 8 : int r=0;
547 8 : if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
548 8 : comm.Bcast(r,0);
549 8 : if(r>0) hillsfname_tmp="/dev/null";
550 8 : ofile->enforceSuffix("");
551 : }
552 12 : ofile->open(hillsfname_tmp);
553 12 : if(fmt.length()>0) ofile->fmtField(fmt);
554 12 : ofile->addConstantField("multivariate");
555 12 : if(doInt_[i]) {
556 0 : ofile->addConstantField("lower_int").printField("lower_int",lowI_[i]);
557 0 : ofile->addConstantField("upper_int").printField("upper_int",uppI_[i]);
558 : }
559 12 : ofile->setHeavyFlush();
560 : // output periodicities of variables
561 12 : ofile->setupPrintValue( getPntrToArgument(i) );
562 : // push back
563 12 : hillsOfiles_.push_back(ofile);
564 12 : }
565 :
566 : // Dump grid to files
567 6 : if(wgridstride_ > 0) {
568 3 : for(unsigned i = 0; i < gridfilenames_.size(); ++i) {
569 2 : OFile *ofile = new OFile();
570 2 : ofile->link(*this);
571 2 : string gridfname_tmp = gridfilenames_[i];
572 2 : if(multiple_w) {
573 0 : int r = 0;
574 0 : if(comm.Get_rank() == 0) {
575 0 : r = multi_sim_comm.Get_rank();
576 : }
577 0 : comm.Bcast(r, 0);
578 0 : if(r>0) {
579 0 : gridfname_tmp = "/dev/null";
580 : }
581 0 : ofile->enforceSuffix("");
582 : }
583 2 : ofile->open(gridfname_tmp);
584 2 : ofile->setHeavyFlush();
585 2 : gridfiles_.push_back(ofile);
586 2 : }
587 : }
588 :
589 6 : log<<" Bibliography "<<plumed.cite("Pfaendtner and Bonomi. J. Chem. Theory Comput. 11, 5062 (2015)");
590 6 : if(doInt_[0]) log<<plumed.cite(
591 0 : "Baftizadeh, Cossio, Pietrucci, and Laio, Curr. Phys. Chem. 2, 79 (2012)");
592 6 : if(multiple_w) log<<plumed.cite(
593 4 : "Raiteri, Laio, Gervasio, Micheletti, and Parrinello, J. Phys. Chem. B 110, 3533 (2006)");
594 12 : log<<"\n";
595 6 : }
596 :
597 0 : void PBMetaD::readGaussians(unsigned iarg, IFile *ifile) {
598 0 : vector<double> center(1);
599 0 : vector<double> sigma(1);
600 : double height;
601 0 : int nhills=0;
602 :
603 0 : std::vector<Value> tmpvalues;
604 0 : tmpvalues.push_back( Value( this, getPntrToArgument(iarg)->getName(), false ) );
605 :
606 0 : while(scanOneHill(iarg,ifile,tmpvalues,center,sigma,height)) {
607 : ;
608 0 : nhills++;
609 0 : if(welltemp_) {height*=(biasf_-1.0)/biasf_;}
610 0 : addGaussian(iarg, Gaussian(center,sigma,height));
611 : }
612 0 : log.printf(" %d Gaussians read\n",nhills);
613 0 : }
614 :
615 0 : bool PBMetaD::readChunkOfGaussians(unsigned iarg, IFile *ifile, unsigned n) {
616 0 : vector<double> center(1);
617 0 : vector<double> sigma(1);
618 : double height;
619 0 : unsigned nhills=0;
620 0 : std::vector<Value> tmpvalues;
621 0 : tmpvalues.push_back( Value( this, getPntrToArgument(iarg)->getName(), false ) );
622 :
623 0 : while(scanOneHill(iarg,ifile,tmpvalues,center,sigma,height)) {
624 : ;
625 0 : if(welltemp_) {height*=(biasf_-1.0)/biasf_;}
626 0 : addGaussian(iarg, Gaussian(center,sigma,height));
627 0 : if(nhills==n) {
628 0 : log.printf(" %u Gaussians read\n",nhills);
629 0 : return true;
630 : }
631 0 : nhills++;
632 : }
633 0 : log.printf(" %u Gaussians read\n",nhills);
634 0 : return false;
635 : }
636 :
637 192 : void PBMetaD::writeGaussian(unsigned iarg, const Gaussian& hill, OFile *ofile) {
638 192 : ofile->printField("time",getTimeStep()*getStep());
639 192 : ofile->printField(getPntrToArgument(iarg),hill.center[0]);
640 192 : ofile->printField("multivariate","false");
641 192 : ofile->printField("sigma_"+getPntrToArgument(iarg)->getName(),hill.sigma[0]);
642 192 : double height=hill.height;
643 192 : if(welltemp_) {height *= biasf_/(biasf_-1.0);}
644 192 : ofile->printField("height",height);
645 192 : ofile->printField("biasf",biasf_);
646 192 : ofile->printField();
647 192 : }
648 :
649 192 : void PBMetaD::addGaussian(unsigned iarg, const Gaussian& hill) {
650 192 : if(!grid_) {hills_[iarg].push_back(hill);}
651 : else {
652 8 : vector<unsigned> nneighb=getGaussianSupport(iarg, hill);
653 16 : vector<Grid::index_t> neighbors=BiasGrids_[iarg]->getNeighbors(hill.center,nneighb);
654 16 : vector<double> der(1);
655 16 : vector<double> xx(1);
656 8 : if(comm.Get_size()==1) {
657 592 : for(unsigned i=0; i<neighbors.size(); ++i) {
658 584 : Grid::index_t ineigh=neighbors[i];
659 584 : der[0]=0.0;
660 584 : BiasGrids_[iarg]->getPoint(ineigh,xx);
661 584 : double bias=evaluateGaussian(iarg,xx,hill,&der[0]);
662 584 : BiasGrids_[iarg]->addValueAndDerivatives(ineigh,bias,der);
663 : }
664 : } else {
665 0 : unsigned stride=comm.Get_size();
666 0 : unsigned rank=comm.Get_rank();
667 0 : vector<double> allder(neighbors.size(),0.0);
668 0 : vector<double> allbias(neighbors.size(),0.0);
669 0 : for(unsigned i=rank; i<neighbors.size(); i+=stride) {
670 0 : Grid::index_t ineigh=neighbors[i];
671 0 : BiasGrids_[iarg]->getPoint(ineigh,xx);
672 0 : allbias[i]=evaluateGaussian(iarg,xx,hill,&allder[i]);
673 : }
674 0 : comm.Sum(allbias);
675 0 : comm.Sum(allder);
676 0 : for(unsigned i=0; i<neighbors.size(); ++i) {
677 0 : Grid::index_t ineigh=neighbors[i];
678 0 : der[0]=allder[i];
679 0 : BiasGrids_[iarg]->addValueAndDerivatives(ineigh,allbias[i],der);
680 0 : }
681 8 : }
682 : }
683 192 : }
684 :
685 8 : vector<unsigned> PBMetaD::getGaussianSupport(unsigned iarg, const Gaussian& hill) {
686 8 : vector<unsigned> nneigh;
687 8 : const double cutoff=sqrt(2.0*DP2CUTOFF)*hill.sigma[0];
688 8 : if(doInt_[iarg]) {
689 0 : if(hill.center[0]+cutoff > uppI_[iarg] || hill.center[0]-cutoff < lowI_[iarg]) {
690 : // in this case, we updated the entire grid to avoid problems
691 0 : return BiasGrids_[iarg]->getNbin();
692 : } else {
693 0 : nneigh.push_back( static_cast<unsigned>(ceil(cutoff/BiasGrids_[iarg]->getDx()[0])));
694 0 : return nneigh;
695 : }
696 : }
697 :
698 8 : nneigh.push_back( static_cast<unsigned>(ceil(cutoff/BiasGrids_[iarg]->getDx()[0])) );
699 :
700 8 : return nneigh;
701 : }
702 :
703 220 : double PBMetaD::getBiasAndDerivatives(unsigned iarg, const vector<double>& cv, double* der)
704 : {
705 220 : double bias=0.0;
706 220 : if(!grid_) {
707 202 : unsigned stride=comm.Get_size();
708 202 : unsigned rank=comm.Get_rank();
709 1106 : for(unsigned i=rank; i<hills_[iarg].size(); i+=stride) {
710 904 : bias += evaluateGaussian(iarg,cv,hills_[iarg][i],der);
711 : }
712 202 : comm.Sum(bias);
713 202 : if(der) comm.Sum(der,1);
714 : } else {
715 18 : if(der) {
716 10 : vector<double> vder(1);
717 10 : bias = BiasGrids_[iarg]->getValueAndDerivatives(cv,vder);
718 10 : der[0] = vder[0];
719 : } else {
720 8 : bias = BiasGrids_[iarg]->getValue(cv);
721 : }
722 : }
723 :
724 220 : return bias;
725 : }
726 :
727 1488 : double PBMetaD::evaluateGaussian(unsigned iarg, const vector<double>& cv, const Gaussian& hill, double* der)
728 : {
729 1488 : double bias=0.0;
730 : // I use a pointer here because cv is const (and should be const)
731 : // but when using doInt it is easier to locally replace cv[0] with
732 : // the upper/lower limit in case it is out of range
733 1488 : const double *pcv=NULL;
734 : double tmpcv[1]; // tmp array with cv (to be used with doInt_)
735 1488 : tmpcv[0]=cv[0];
736 1488 : bool isOutOfInt = false;
737 1488 : if(doInt_[iarg]) {
738 0 : if(cv[0]<lowI_[iarg]) { tmpcv[0]=lowI_[iarg]; isOutOfInt = true; }
739 0 : else if(cv[0]>uppI_[iarg]) { tmpcv[0]=uppI_[iarg]; isOutOfInt = true; }
740 : }
741 1488 : pcv=&(tmpcv[0]);
742 1488 : double dp = difference(iarg, hill.center[0],pcv[0]) / hill.sigma[0];
743 1488 : double dp2 = 0.5 * dp * dp;
744 1488 : if(dp2<DP2CUTOFF) {
745 1472 : bias = hill.height*exp(-dp2);
746 1472 : if(der &&!isOutOfInt) {
747 1020 : der[0] += -bias * dp / hill.sigma[0];
748 : }
749 : }
750 1488 : return bias;
751 : }
752 :
753 58 : void PBMetaD::calculate()
754 : {
755 58 : vector<double> cv(1);
756 : double der[1];
757 116 : vector<double> bias(getNumberOfArguments());
758 116 : vector<double> deriv(getNumberOfArguments());
759 :
760 58 : double ncv = (double) getNumberOfArguments();
761 58 : double bmin = 1.0e+19;
762 174 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
763 116 : cv[0] = getArgument(i);
764 116 : der[0] = 0.0;
765 116 : bias[i] = getBiasAndDerivatives(i, cv, der);
766 116 : deriv[i] = der[0];
767 116 : if(bias[i] < bmin) bmin = bias[i];
768 : }
769 58 : double ene = 0.;
770 174 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
771 116 : ene += exp((-bias[i]+bmin)/kbt_);
772 : }
773 :
774 : // set Forces
775 174 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
776 116 : const double f = - exp((-bias[i]+bmin)/kbt_) / (ene) * deriv[i];
777 116 : setOutputForce(i, f);
778 : }
779 :
780 : // set bias
781 58 : ene = -kbt_ * (std::log(ene) - std::log(ncv)) + bmin;
782 116 : setBias(ene);
783 58 : }
784 :
785 58 : void PBMetaD::update()
786 : {
787 : // adding hills criteria
788 : bool nowAddAHill;
789 58 : if(getStep()%stride_==0 && !isFirstStep ) {nowAddAHill=true;} else {nowAddAHill=false; isFirstStep=false;}
790 :
791 58 : if(nowAddAHill) {
792 : // get all CVs value
793 52 : vector<double> cv(getNumberOfArguments());
794 52 : for(unsigned i=0; i<getNumberOfArguments(); ++i) cv[i] = getArgument(i);
795 : // get all biases and heights
796 104 : vector<double> bias(getNumberOfArguments());
797 104 : vector<double> height(getNumberOfArguments());
798 104 : vector<double> cv_tmp(1);
799 104 : vector<double> sigma_tmp(1);
800 52 : double bmin = 1.0e+19;
801 156 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
802 104 : cv_tmp[0] = cv[i];
803 104 : bias[i] = getBiasAndDerivatives(i, cv_tmp);
804 104 : if(bias[i] < bmin) bmin = bias[i];
805 : }
806 52 : double norm = 0.0;
807 : // calculate heights and norm
808 156 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
809 104 : double h = exp((-bias[i]+bmin)/kbt_);
810 104 : norm += h;
811 104 : height[i] = h;
812 : }
813 : // normalize and apply welltemp correction
814 156 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
815 104 : height[i] *= height0_ / norm;
816 104 : if(welltemp_) height[i] *= exp(-bias[i]/(kbt_*(biasf_-1.0)));
817 : }
818 :
819 : // Multiple walkers: share hills and add them all
820 52 : if(multiple_w) {
821 : // Allocate arrays to store all walkers hills
822 44 : std::vector<double> all_cv(nw_*cv.size(), 0.0);
823 88 : std::vector<double> all_height(nw_*height.size(), 0.0);
824 44 : if(comm.Get_rank()==0) {
825 : // fill in value
826 66 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
827 44 : unsigned j = mw_ * getNumberOfArguments() + i;
828 44 : all_cv[j] = cv[i];
829 44 : all_height[j] = height[i];
830 : }
831 : // Communicate (only root)
832 22 : multi_sim_comm.Sum(&all_cv[0], all_cv.size());
833 22 : multi_sim_comm.Sum(&all_height[0], all_height.size());
834 : }
835 : // Share info with group members
836 44 : comm.Sum(&all_cv[0], all_cv.size());
837 44 : comm.Sum(&all_height[0], all_height.size());
838 : // now add hills one by one
839 132 : for(unsigned j=0; j<nw_; ++j) {
840 264 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
841 176 : cv_tmp[0] = all_cv[j*cv.size()+i];
842 176 : sigma_tmp[0] = sigma0_[i];
843 : // new Gaussian
844 176 : Gaussian newhill = Gaussian(cv_tmp, sigma_tmp, all_height[j*cv.size()+i]);
845 176 : addGaussian(i, newhill);
846 : // print on HILLS file
847 176 : writeGaussian(i, newhill, hillsOfiles_[i]);
848 176 : }
849 44 : }
850 : // just add your own hills
851 : } else {
852 24 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
853 16 : cv_tmp[0] = cv[i];
854 16 : sigma_tmp[0] = sigma0_[i];
855 : // new Gaussian
856 16 : Gaussian newhill = Gaussian(cv_tmp, sigma_tmp, height[i]);
857 16 : addGaussian(i, newhill);
858 : // print on HILLS file
859 16 : writeGaussian(i, newhill, hillsOfiles_[i]);
860 16 : }
861 52 : }
862 : }
863 : // write grid files
864 58 : if(wgridstride_>0 && (getStep()%wgridstride_==0 || getCPT())) {
865 5 : int r = 0;
866 5 : if(multiple_w) {
867 0 : if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
868 0 : comm.Bcast(r,0);
869 : }
870 5 : if(r==0) {
871 15 : for(unsigned i=0; i<gridfiles_.size(); ++i) {
872 10 : gridfiles_[i]->rewind();
873 10 : BiasGrids_[i]->writeToFile(*gridfiles_[i]);
874 10 : gridfiles_[i]->flush();
875 : }
876 : }
877 : }
878 58 : }
879 :
880 : /// takes a pointer to the file and a template string with values v and gives back the next center, sigma and height
881 0 : bool PBMetaD::scanOneHill(unsigned iarg, IFile *ifile, vector<Value> &tmpvalues, vector<double> ¢er, vector<double> &sigma, double &height) {
882 : double dummy;
883 :
884 0 : if(ifile->scanField("time",dummy)) {
885 0 : ifile->scanField( &tmpvalues[0] );
886 0 : if( tmpvalues[0].isPeriodic() && ! getPntrToArgument(iarg)->isPeriodic() ) {
887 0 : error("in hills file periodicity for variable " + tmpvalues[0].getName() + " does not match periodicity in input");
888 0 : } else if( tmpvalues[0].isPeriodic() ) {
889 0 : std::string imin, imax; tmpvalues[0].getDomain( imin, imax );
890 0 : std::string rmin, rmax; getPntrToArgument(iarg)->getDomain( rmin, rmax );
891 0 : if( imin!=rmin || imax!=rmax ) {
892 0 : error("in hills file periodicity for variable " + tmpvalues[0].getName() + " does not match periodicity in input");
893 0 : }
894 : }
895 0 : center[0]=tmpvalues[0].get();
896 : // scan for multivariate label: record the actual file position so to eventually rewind
897 0 : std::string sss;
898 0 : ifile->scanField("multivariate",sss);
899 0 : ifile->scanField("sigma_"+getPntrToArgument(iarg)->getName(),sigma[0]);
900 0 : ifile->scanField("height",height);
901 0 : ifile->scanField("biasf",dummy);
902 0 : if(ifile->FieldExist("clock")) ifile->scanField("clock",dummy);
903 0 : if(ifile->FieldExist("lower_int")) ifile->scanField("lower_int",dummy);
904 0 : if(ifile->FieldExist("upper_int")) ifile->scanField("upper_int",dummy);
905 0 : ifile->scanField();
906 0 : return true;
907 : } else {
908 0 : return false;
909 : };
910 : }
911 :
912 : }
913 2523 : }
914 :
|
