Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2019 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 "core/FlexibleBin.h"
30 : #include "tools/Matrix.h"
31 : #include "tools/Random.h"
32 : #include <string>
33 : #include <cstring>
34 : #include "tools/File.h"
35 : #include <iostream>
36 : #include <limits>
37 : #include <ctime>
38 :
39 : #define DP2CUTOFF 6.25
40 :
41 : using namespace std;
42 :
43 :
44 : namespace PLMD {
45 : namespace bias {
46 :
47 : //+PLUMEDOC BIAS PBMETAD
48 : /*
49 : Used to performed Parallel Bias MetaDynamics.
50 :
51 : This action activate Parallel Bias MetaDynamics (PBMetaD) \cite pbmetad, a version of MetaDynamics \cite metad in which
52 : multiple low-dimensional bias potentials are applied in parallel.
53 : In the current implementation, these have the form of mono-dimensional MetaDynamics bias
54 : potentials:
55 :
56 : \f[
57 : {V(s_1,t), ..., V(s_N,t)}
58 : \f]
59 :
60 : where:
61 :
62 : \f[
63 : V(s_i,t) = \sum_{ k \tau < t} W_i(k \tau)
64 : \exp\left(
65 : - \frac{(s_i-s_i^{(0)}(k \tau))^2}{2\sigma_i^2}
66 : \right).
67 : \f]
68 :
69 : To ensure the convergence of each mono-dimensional bias potential to the corresponding free energy,
70 : at each deposition step the Gaussian heights are multiplied by the so-called conditional term:
71 :
72 : \f[
73 : W_i(k \tau)=W_0 \frac{\exp\left(
74 : - \frac{V(s_i,k \tau)}{k_B T}
75 : \right)}{\sum_{i=1}^N
76 : \exp\left(
77 : - \frac{V(s_i,k \tau)}{k_B T}
78 : \right)}
79 : \f]
80 :
81 : where \f$W_0\f$ is the initial Gaussian height.
82 :
83 : The PBMetaD bias potential is defined by:
84 :
85 : \f[
86 : V_{PB}(\vec{s},t) = -k_B T \log{\sum_{i=1}^N
87 : \exp\left(
88 : - \frac{V(s_i,t)}{k_B T}
89 : \right)}.
90 : \f]
91 :
92 :
93 : Information on the Gaussian functions that build each bias potential are printed to
94 : multiple HILLS files, which
95 : are used both to restart the calculation and to reconstruct the mono-dimensional
96 : free energies as a function of the corresponding CVs.
97 : These can be reconstructed using the \ref sum_hills utility because the final bias is given
98 : by:
99 :
100 : \f[
101 : V(s_i) = -F(s_i)
102 : \f]
103 :
104 : Currently, only a subset of the \ref METAD options are available in PBMetaD.
105 :
106 : The bias potentials can be stored on a grid to increase performances of long PBMetaD simulations.
107 : You should
108 : provide either the number of bins for every collective variable (GRID_BIN) or
109 : the desired grid spacing (GRID_SPACING). In case you provide both PLUMED will use
110 : the most conservative choice (highest number of bins) for each dimension.
111 : In case you do not provide any information about bin size (neither GRID_BIN nor GRID_SPACING)
112 : and if Gaussian width is fixed PLUMED will use 1/5 of the Gaussian width as grid spacing.
113 : This default choice should be reasonable for most applications.
114 :
115 : Another option that is available is well-tempered metadynamics \cite Barducci:2008. In this
116 : variant of PBMetaD the heights of the Gaussian hills are rescaled at each step by the
117 : additional well-tempered metadynamics term.
118 : This ensures that each bias converges more smoothly. It should be noted that, in the case of well-tempered metadynamics, in
119 : the output printed the Gaussian height is re-scaled using the bias factor.
120 : Also notice that with well-tempered metadynamics the HILLS files do not contain the bias,
121 : but the negative of the free-energy estimate. This choice has the advantage that
122 : one can restart a simulation using a different value for the \f$\Delta T\f$. The applied bias will be scaled accordingly.
123 :
124 : Note that you can use here also the flexible gaussian approach \cite Branduardi:2012dl
125 : in which you can adapt the gaussian to the extent of Cartesian space covered by a variable or
126 : to the space in collective variable covered in a given time. In this case the width of the deposited
127 : gaussian potential is denoted by one value only that is a Cartesian space (ADAPTIVE=GEOM) or a time
128 : (ADAPTIVE=DIFF). Note that a specific integration technique for the deposited gaussians
129 : should be used in this case. Check the documentation for utility sum_hills.
130 :
131 : With the keyword INTERVAL one changes the metadynamics algorithm setting the bias force equal to zero
132 : outside boundary \cite baftizadeh2012protein. If, for example, metadynamics is performed on a CV s and one is interested only
133 : to the free energy for s > sw, the history dependent potential is still updated according to the above
134 : equations but the metadynamics force is set to zero for s < sw. Notice that Gaussians are added also
135 : if s < sw, as the tails of these Gaussians influence VG in the relevant region s > sw. In this way, the
136 : force on the system in the region s > sw comes from both metadynamics and the force field, in the region
137 : s < sw only from the latter. This approach allows obtaining a history-dependent bias potential VG that
138 : fluctuates around a stable estimator, equal to the negative of the free energy far enough from the
139 : boundaries. Note that:
140 : - It works only for one-dimensional biases;
141 : - It works both with and without GRID;
142 : - The interval limit sw in a region where the free energy derivative is not large;
143 : - If in the region outside the limit sw the system has a free energy minimum, the INTERVAL keyword should
144 : be used together with a \ref UPPER_WALLS or \ref LOWER_WALLS at sw.
145 :
146 : Multiple walkers \cite multiplewalkers can also be used. See below the examples.
147 :
148 : \par Examples
149 :
150 : The following input is for PBMetaD calculation using as
151 : collective variables the distance between atoms 3 and 5
152 : and the distance between atoms 2 and 4. The value of the CVs and
153 : the PBMetaD bias potential are written to the COLVAR file every 100 steps.
154 : \plumedfile
155 : DISTANCE ATOMS=3,5 LABEL=d1
156 : DISTANCE ATOMS=2,4 LABEL=d2
157 : PBMETAD ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3 PACE=500 LABEL=pb FILE=HILLS_d1,HILLS_d2
158 : PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
159 : \endplumedfile
160 : (See also \ref DISTANCE and \ref PRINT).
161 :
162 : \par
163 : If you use well-tempered metadynamics, you should specify a single biasfactor and initial
164 : Gaussian height.
165 : \plumedfile
166 : DISTANCE ATOMS=3,5 LABEL=d1
167 : DISTANCE ATOMS=2,4 LABEL=d2
168 : PBMETAD ...
169 : ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3
170 : PACE=500 BIASFACTOR=8 LABEL=pb
171 : FILE=HILLS_d1,HILLS_d2
172 : ... PBMETAD
173 : PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
174 : \endplumedfile
175 :
176 : \par
177 : The following input enables the MPI version of multiple-walkers.
178 : \plumedfile
179 : DISTANCE ATOMS=3,5 LABEL=d1
180 : DISTANCE ATOMS=2,4 LABEL=d2
181 : PBMETAD ...
182 : ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3
183 : PACE=500 BIASFACTOR=8 LABEL=pb
184 : FILE=HILLS_d1,HILLS_d2
185 : WALKERS_MPI
186 : ... PBMETAD
187 : PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
188 : \endplumedfile
189 :
190 : \par
191 : The disk version of multiple-walkers can be
192 : enabled by setting the number of walker used, the id of the
193 : current walker which interprets the input file, the directory where the
194 : hills containing files resides, and the frequency to read the other walkers.
195 : Here is an example
196 : \plumedfile
197 : DISTANCE ATOMS=3,5 LABEL=d1
198 : DISTANCE ATOMS=2,4 LABEL=d2
199 : PBMETAD ...
200 : ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3
201 : PACE=500 BIASFACTOR=8 LABEL=pb
202 : FILE=HILLS_d1,HILLS_d2
203 : WALKERS_N=10
204 : WALKERS_ID=3
205 : WALKERS_DIR=../
206 : WALKERS_RSTRIDE=100
207 : ... PBMETAD
208 : PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
209 : \endplumedfile
210 : where WALKERS_N is the total number of walkers, WALKERS_ID is the
211 : id of the present walker (starting from 0 ) and the WALKERS_DIR is the directory
212 : where all the walkers are located. WALKERS_RSTRIDE is the number of step between
213 : one update and the other.
214 :
215 : */
216 : //+ENDPLUMEDOC
217 :
218 : class PBMetaD : public Bias {
219 :
220 : private:
221 10686 : struct Gaussian {
222 : vector<double> center;
223 : vector<double> sigma;
224 : double height;
225 : bool multivariate; // this is required to discriminate the one dimensional case
226 : vector<double> invsigma;
227 896 : Gaussian(const vector<double> & center,const vector<double> & sigma, double height, bool multivariate):
228 896 : center(center),sigma(sigma),height(height),multivariate(multivariate),invsigma(sigma) {
229 : // to avoid troubles from zero element in flexible hills
230 5376 : for(unsigned i=0; i<invsigma.size(); ++i) abs(invsigma[i])>1.e-20?invsigma[i]=1.0/invsigma[i]:0.;
231 896 : }
232 : };
233 : vector<double> sigma0_;
234 : vector<double> sigma0min_;
235 : vector<double> sigma0max_;
236 : vector< vector<Gaussian> > hills_;
237 : vector<OFile*> hillsOfiles_;
238 : vector<OFile*> gridfiles_;
239 : vector<Grid*> BiasGrids_;
240 : bool grid_;
241 : double height0_;
242 : double biasf_;
243 : double kbt_;
244 : int stride_;
245 : int wgridstride_;
246 : bool welltemp_;
247 : int mw_n_;
248 : string mw_dir_;
249 : int mw_id_;
250 : int mw_rstride_;
251 : bool walkers_mpi;
252 : unsigned mpi_nw_;
253 : unsigned mpi_id_;
254 : vector<string> hillsfname;
255 : vector<IFile*> ifiles;
256 : vector<string> ifilesnames;
257 : vector<double> uppI_;
258 : vector<double> lowI_;
259 : vector<bool> doInt_;
260 : int adaptive_;
261 : vector<FlexibleBin> flexbin;
262 : bool isFirstStep;
263 : // variable for selector
264 : string selector_;
265 : bool do_select_;
266 : unsigned select_value_;
267 : unsigned current_value_;
268 :
269 : void readGaussians(unsigned iarg, IFile*);
270 : bool readChunkOfGaussians(unsigned iarg, IFile *ifile, unsigned n);
271 : void writeGaussian(unsigned iarg, const Gaussian&, OFile*);
272 : void addGaussian(unsigned iarg, const Gaussian&);
273 : double getBiasAndDerivatives(unsigned iarg, const vector<double>&, double* der=NULL);
274 : double evaluateGaussian(unsigned iarg, const vector<double>&, const Gaussian&,double* der=NULL);
275 : vector<unsigned> getGaussianSupport(unsigned iarg, const Gaussian&);
276 : bool scanOneHill(unsigned iarg, IFile *ifile, vector<Value> &v, vector<double> ¢er, vector<double> &sigma, double &height, bool &multivariate);
277 : std::string fmt;
278 :
279 : public:
280 : explicit PBMetaD(const ActionOptions&);
281 : ~PBMetaD();
282 : void calculate();
283 : void update();
284 : static void registerKeywords(Keywords& keys);
285 : bool checkNeedsGradients()const;
286 : };
287 :
288 6482 : PLUMED_REGISTER_ACTION(PBMetaD,"PBMETAD")
289 :
290 31 : void PBMetaD::registerKeywords(Keywords& keys) {
291 31 : Bias::registerKeywords(keys);
292 62 : keys.use("ARG");
293 124 : keys.add("compulsory","SIGMA","the widths of the Gaussian hills");
294 124 : keys.add("compulsory","PACE","the frequency for hill addition, one for all biases");
295 124 : keys.add("optional","FILE","files in which the lists of added hills are stored, default names are assigned using arguments if FILE is not found");
296 124 : keys.add("optional","HEIGHT","the height of the Gaussian hills, one for all biases. Compulsory unless TAU, TEMP and BIASFACTOR are given");
297 124 : keys.add("optional","FMT","specify format for HILLS files (useful for decrease the number of digits in regtests)");
298 124 : keys.add("optional","BIASFACTOR","use well tempered metadynamics with this biasfactor, one for all biases. Please note you must also specify temp");
299 124 : keys.add("optional","TEMP","the system temperature - this is only needed if you are doing well-tempered metadynamics");
300 124 : keys.add("optional","TAU","in well tempered metadynamics, sets height to (kb*DeltaT*pace*timestep)/tau");
301 124 : keys.add("optional","GRID_RFILES", "read grid for the bias");
302 124 : keys.add("optional","GRID_WSTRIDE", "frequency for dumping the grid");
303 124 : keys.add("optional","GRID_WFILES", "dump grid for the bias, default names are used if GRID_WSTRIDE is used without GRID_WFILES.");
304 124 : keys.add("optional","GRID_MIN","the lower bounds for the grid");
305 124 : keys.add("optional","GRID_MAX","the upper bounds for the grid");
306 124 : keys.add("optional","GRID_BIN","the number of bins for the grid");
307 124 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
308 93 : keys.addFlag("GRID_SPARSE",false,"use a sparse grid to store hills");
309 93 : keys.addFlag("GRID_NOSPLINE",false,"don't use spline interpolation with grids");
310 124 : keys.add("optional","SELECTOR", "add forces and do update based on the value of SELECTOR");
311 124 : keys.add("optional","SELECTOR_ID", "value of SELECTOR");
312 124 : keys.add("optional","WALKERS_ID", "walker id");
313 124 : keys.add("optional","WALKERS_N", "number of walkers");
314 124 : keys.add("optional","WALKERS_DIR", "shared directory with the hills files from all the walkers");
315 124 : keys.add("optional","WALKERS_RSTRIDE","stride for reading hills files");
316 124 : keys.add("optional","INTERVAL_MIN","monodimensional lower limits, outside the limits the system will not feel the biasing force.");
317 124 : keys.add("optional","INTERVAL_MAX","monodimensional upper limits, outside the limits the system will not feel the biasing force.");
318 124 : keys.add("optional","ADAPTIVE","use a geometric (=GEOM) or diffusion (=DIFF) based hills width scheme. Sigma is one number that has distance units or timestep dimensions");
319 124 : keys.add("optional","SIGMA_MAX","the upper bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
320 124 : keys.add("optional","SIGMA_MIN","the lower bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
321 93 : keys.addFlag("WALKERS_MPI",false,"Switch on MPI version of multiple walkers - not compatible with WALKERS_* options other than WALKERS_DIR");
322 62 : keys.use("RESTART");
323 62 : keys.use("UPDATE_FROM");
324 62 : keys.use("UPDATE_UNTIL");
325 31 : }
326 :
327 210 : PBMetaD::~PBMetaD() {
328 66 : for(unsigned i=0; i<BiasGrids_.size(); ++i) delete BiasGrids_[i];
329 240 : for(unsigned i=0; i<hillsOfiles_.size(); ++i) {
330 60 : hillsOfiles_[i]->close();
331 60 : delete hillsOfiles_[i];
332 : }
333 30 : if(wgridstride_ > 0) {
334 8 : for(unsigned i=0; i<gridfiles_.size(); ++i) {
335 2 : gridfiles_[i]->close();
336 2 : delete gridfiles_[i];
337 : }
338 : }
339 : // close files
340 240 : for(unsigned i=0; i<ifiles.size(); ++i) {
341 60 : if(ifiles[i]->isOpen()) ifiles[i]->close();
342 60 : delete ifiles[i];
343 : }
344 60 : }
345 :
346 30 : PBMetaD::PBMetaD(const ActionOptions& ao):
347 : PLUMED_BIAS_INIT(ao),
348 : grid_(false), height0_(std::numeric_limits<double>::max()),
349 : biasf_(1.0), kbt_(0.0), stride_(0), wgridstride_(0), welltemp_(false),
350 : mw_n_(1), mw_dir_(""), mw_id_(0), mw_rstride_(1),
351 : walkers_mpi(false), mpi_nw_(0),
352 : adaptive_(FlexibleBin::none),
353 : isFirstStep(true),
354 480 : do_select_(false)
355 : {
356 : // parse the flexible hills
357 : string adaptiveoption;
358 : adaptiveoption="NONE";
359 60 : parse("ADAPTIVE",adaptiveoption);
360 30 : if(adaptiveoption=="GEOM") {
361 0 : log.printf(" Uses Geometry-based hills width: sigma must be in distance units and only one sigma is needed\n");
362 0 : adaptive_=FlexibleBin::geometry;
363 30 : } else if(adaptiveoption=="DIFF") {
364 0 : log.printf(" Uses Diffusion-based hills width: sigma must be in timesteps and only one sigma is needed\n");
365 0 : adaptive_=FlexibleBin::diffusion;
366 30 : } else if(adaptiveoption=="NONE") {
367 30 : adaptive_=FlexibleBin::none;
368 : } else {
369 0 : error("I do not know this type of adaptive scheme");
370 : }
371 :
372 60 : parse("FMT",fmt);
373 :
374 : // parse the sigma
375 60 : parseVector("SIGMA",sigma0_);
376 30 : if(adaptive_==FlexibleBin::none) {
377 : // if you use normal sigma you need one sigma per argument
378 30 : if( sigma0_.size()!=getNumberOfArguments() ) error("number of arguments does not match number of SIGMA parameters");
379 : } else {
380 : // if you use flexible hills you need one sigma
381 0 : if(sigma0_.size()!=1) {
382 0 : error("If you choose ADAPTIVE you need only one sigma according to your choice of type (GEOM/DIFF)");
383 : }
384 : // if adaptive then the number must be an integer
385 0 : if(adaptive_==FlexibleBin::diffusion) {
386 0 : if(int(sigma0_[0])-sigma0_[0]>1.e-9 || int(sigma0_[0])-sigma0_[0] <-1.e-9 || int(sigma0_[0])<1 ) {
387 0 : error("In case of adaptive hills with diffusion, the sigma must be an integer which is the number of timesteps\n");
388 : }
389 : }
390 : // here evtl parse the sigma min and max values
391 0 : parseVector("SIGMA_MIN",sigma0min_);
392 0 : if(sigma0min_.size()>0 && sigma0min_.size()!=getNumberOfArguments()) {
393 0 : error("the number of SIGMA_MIN values be the same of the number of the arguments");
394 0 : } else if(sigma0min_.size()==0) {
395 0 : sigma0min_.resize(getNumberOfArguments());
396 0 : for(unsigned i=0; i<getNumberOfArguments(); i++) {sigma0min_[i]=-1.;}
397 : }
398 :
399 0 : parseVector("SIGMA_MAX",sigma0max_);
400 0 : if(sigma0max_.size()>0 && sigma0max_.size()!=getNumberOfArguments()) {
401 0 : error("the number of SIGMA_MAX values be the same of the number of the arguments");
402 0 : } else if(sigma0max_.size()==0) {
403 0 : sigma0max_.resize(getNumberOfArguments());
404 0 : for(unsigned i=0; i<getNumberOfArguments(); i++) {sigma0max_[i]=-1.;}
405 : }
406 :
407 0 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
408 : vector<double> tmp_smin, tmp_smax;
409 0 : tmp_smin.resize(1,sigma0min_[i]);
410 0 : tmp_smax.resize(1,sigma0max_[i]);
411 0 : flexbin.push_back(FlexibleBin(adaptive_,this,i,sigma0_[0],tmp_smin,tmp_smax));
412 : }
413 : }
414 :
415 : // note: HEIGHT is not compulsory, since one could use the TAU keyword, see below
416 60 : parse("HEIGHT",height0_);
417 60 : parse("PACE",stride_);
418 30 : if(stride_<=0) error("frequency for hill addition is nonsensical");
419 :
420 60 : parseVector("FILE",hillsfname);
421 30 : if(hillsfname.size()==0) {
422 0 : for(unsigned i=0; i<getNumberOfArguments(); i++) hillsfname.push_back("HILLS."+getPntrToArgument(i)->getName());
423 : }
424 30 : if( hillsfname.size()!=getNumberOfArguments() ) {
425 0 : error("number of FILE arguments does not match number of HILLS files");
426 : }
427 :
428 60 : parse("BIASFACTOR",biasf_);
429 30 : if( biasf_<1.0 ) error("well tempered bias factor is nonsensical");
430 30 : double temp=0.0;
431 60 : parse("TEMP",temp);
432 60 : if(temp>0.0) kbt_=plumed.getAtoms().getKBoltzmann()*temp;
433 0 : else kbt_=plumed.getAtoms().getKbT();
434 30 : if(biasf_>1.0) {
435 29 : if(kbt_==0.0) error("Unless the MD engine passes the temperature to plumed, with well-tempered metad you must specify it using TEMP");
436 29 : welltemp_=true;
437 : }
438 30 : double tau=0.0;
439 60 : parse("TAU",tau);
440 30 : if(tau==0.0) {
441 30 : if(height0_==std::numeric_limits<double>::max()) error("At least one between HEIGHT and TAU should be specified");
442 : // if tau is not set, we compute it here from the other input parameters
443 30 : if(welltemp_) tau=(kbt_*(biasf_-1.0))/height0_*getTimeStep()*stride_;
444 : } else {
445 0 : if(!welltemp_)error("TAU only makes sense in well-tempered metadynamics");
446 0 : if(height0_!=std::numeric_limits<double>::max()) error("At most one between HEIGHT and TAU should be specified");
447 0 : height0_=(kbt_*(biasf_-1.0))/tau*getTimeStep()*stride_;
448 : }
449 :
450 : // Multiple walkers
451 60 : parse("WALKERS_N",mw_n_);
452 60 : parse("WALKERS_ID",mw_id_);
453 30 : if(mw_n_<=mw_id_) error("walker ID should be a numerical value less than the total number of walkers");
454 60 : parse("WALKERS_DIR",mw_dir_);
455 60 : parse("WALKERS_RSTRIDE",mw_rstride_);
456 :
457 : // MPI version
458 60 : parseFlag("WALKERS_MPI",walkers_mpi);
459 :
460 : // Grid file
461 60 : parse("GRID_WSTRIDE",wgridstride_);
462 30 : vector<string> gridfilenames_;
463 60 : parseVector("GRID_WFILES",gridfilenames_);
464 59 : if (wgridstride_ == 0 && gridfilenames_.size() > 0) {
465 0 : error("frequency with which to output grid not specified use GRID_WSTRIDE");
466 : }
467 30 : if(gridfilenames_.size()==0 && wgridstride_ > 0) {
468 0 : for(unsigned i=0; i<getNumberOfArguments(); i++) gridfilenames_.push_back("GRID."+getPntrToArgument(i)->getName());
469 : }
470 31 : if(gridfilenames_.size() > 0 && hillsfname.size() > 0 && gridfilenames_.size() != hillsfname.size())
471 0 : error("number of GRID_WFILES arguments does not match number of HILLS files");
472 :
473 : // Read grid
474 30 : vector<string> gridreadfilenames_;
475 60 : parseVector("GRID_RFILES",gridreadfilenames_);
476 :
477 : // Grid Stuff
478 60 : vector<std::string> gmin(getNumberOfArguments());
479 60 : parseVector("GRID_MIN",gmin);
480 30 : if(gmin.size()!=getNumberOfArguments() && gmin.size()!=0) error("not enough values for GRID_MIN");
481 60 : vector<std::string> gmax(getNumberOfArguments());
482 60 : parseVector("GRID_MAX",gmax);
483 30 : if(gmax.size()!=getNumberOfArguments() && gmax.size()!=0) error("not enough values for GRID_MAX");
484 30 : vector<unsigned> gbin(getNumberOfArguments());
485 : vector<double> gspacing;
486 60 : parseVector("GRID_BIN",gbin);
487 30 : if(gbin.size()!=getNumberOfArguments() && gbin.size()!=0) error("not enough values for GRID_BIN");
488 60 : parseVector("GRID_SPACING",gspacing);
489 30 : if(gspacing.size()!=getNumberOfArguments() && gspacing.size()!=0) error("not enough values for GRID_SPACING");
490 30 : if(gmin.size()!=gmax.size()) error("GRID_MAX and GRID_MIN should be either present or absent");
491 30 : if(gspacing.size()!=0 && gmin.size()==0) error("If GRID_SPACING is present also GRID_MIN should be present");
492 30 : if(gbin.size()!=0 && gmin.size()==0) error("If GRID_SPACING is present also GRID_MIN should be present");
493 30 : if(gmin.size()!=0) {
494 2 : if(gbin.size()==0 && gspacing.size()==0) {
495 1 : if(adaptive_==FlexibleBin::none) {
496 1 : log<<" Binsize not specified, 1/10 of sigma will be be used\n";
497 1 : plumed_assert(sigma0_.size()==getNumberOfArguments());
498 1 : gspacing.resize(getNumberOfArguments());
499 10 : for(unsigned i=0; i<gspacing.size(); i++) gspacing[i]=0.1*sigma0_[i];
500 : } else {
501 : // with adaptive hills and grid a sigma min must be specified
502 0 : for(unsigned i=0; i<sigma0min_.size(); i++) if(sigma0min_[i]<=0) error("When using Adaptive Gaussians on a grid SIGMA_MIN must be specified");
503 0 : log<<" Binsize not specified, 1/5 of sigma_min will be be used\n";
504 0 : gspacing.resize(getNumberOfArguments());
505 0 : for(unsigned i=0; i<gspacing.size(); i++) gspacing[i]=0.2*sigma0min_[i];
506 : }
507 0 : } else if(gspacing.size()!=0 && gbin.size()==0) {
508 0 : log<<" The number of bins will be estimated from GRID_SPACING\n";
509 0 : } else if(gspacing.size()!=0 && gbin.size()!=0) {
510 0 : log<<" You specified both GRID_BIN and GRID_SPACING\n";
511 0 : log<<" The more conservative (highest) number of bins will be used for each variable\n";
512 : }
513 3 : if(gbin.size()==0) gbin.assign(getNumberOfArguments(),1);
514 6 : if(gspacing.size()!=0) for(unsigned i=0; i<getNumberOfArguments(); i++) {
515 : double a,b;
516 4 : Tools::convert(gmin[i],a);
517 2 : Tools::convert(gmax[i],b);
518 4 : unsigned n=((b-a)/gspacing[i])+1;
519 2 : if(gbin[i]<n) gbin[i]=n;
520 : }
521 : }
522 30 : bool sparsegrid=false;
523 60 : parseFlag("GRID_SPARSE",sparsegrid);
524 30 : bool nospline=false;
525 60 : parseFlag("GRID_NOSPLINE",nospline);
526 30 : bool spline=!nospline;
527 30 : if(gbin.size()>0) {grid_=true;}
528 59 : if(!grid_&&gridfilenames_.size() > 0) error("To write a grid you need first to define it!");
529 59 : if(!grid_&&gridreadfilenames_.size() > 0) error("To read a grid you need first to define it!");
530 :
531 30 : doInt_.resize(getNumberOfArguments(),false);
532 : // Interval keyword
533 60 : parseVector("INTERVAL_MIN",lowI_);
534 60 : parseVector("INTERVAL_MAX",uppI_);
535 : // various checks
536 30 : if(lowI_.size()!=uppI_.size()) error("both a lower and an upper limits must be provided with INTERVAL");
537 30 : if(lowI_.size()!=0 && lowI_.size()!=getNumberOfArguments()) error("check number of argument of INTERVAL");
538 60 : for(unsigned i=0; i<lowI_.size(); ++i) {
539 0 : if(uppI_[i]<lowI_[i]) error("The Upper limit must be greater than the Lower limit!");
540 0 : if(getPntrToArgument(i)->isPeriodic()) warning("INTERVAL is not used for periodic variables");
541 : else doInt_[i]=true;
542 : }
543 :
544 : // parse selector stuff
545 60 : parse("SELECTOR", selector_);
546 30 : if(selector_.length()>0) {
547 0 : do_select_ = true;
548 0 : parse("SELECTOR_ID", select_value_);
549 : }
550 :
551 30 : checkRead();
552 :
553 30 : log.printf(" Gaussian width ");
554 30 : if (adaptive_==FlexibleBin::diffusion)log.printf(" (Note: The units of sigma are in timesteps) ");
555 30 : if (adaptive_==FlexibleBin::geometry)log.printf(" (Note: The units of sigma are in dist units) ");
556 240 : for(unsigned i=0; i<sigma0_.size(); ++i) log.printf(" %f",sigma0_[i]);
557 30 : log.printf(" Gaussian height %f\n",height0_);
558 30 : log.printf(" Gaussian deposition pace %d\n",stride_);
559 :
560 30 : log.printf(" Gaussian files ");
561 240 : for(unsigned i=0; i<hillsfname.size(); ++i) log.printf("%s ",hillsfname[i].c_str());
562 30 : log.printf("\n");
563 30 : if(welltemp_) {
564 29 : log.printf(" Well-Tempered Bias Factor %f\n",biasf_);
565 29 : log.printf(" Hills relaxation time (tau) %f\n",tau);
566 29 : log.printf(" KbT %f\n",kbt_);
567 : }
568 :
569 30 : if(do_select_) {
570 0 : log.printf(" Add forces and update bias based on the value of SELECTOR %s\n",selector_.c_str());
571 0 : log.printf(" Id of the SELECTOR for this action %u\n", select_value_);
572 : }
573 :
574 30 : if(mw_n_>1) {
575 0 : if(walkers_mpi) error("MPI version of multiple walkers is not compatible with filesystem version of multiple walkers");
576 0 : log.printf(" %d multiple walkers active\n",mw_n_);
577 0 : log.printf(" walker id %d\n",mw_id_);
578 0 : log.printf(" reading stride %d\n",mw_rstride_);
579 0 : if(mw_dir_!="")log.printf(" directory with hills files %s\n",mw_dir_.c_str());
580 : } else {
581 30 : if(walkers_mpi) {
582 28 : log.printf(" Multiple walkers active using MPI communnication\n");
583 28 : if(mw_dir_!="")log.printf(" directory with hills files %s\n",mw_dir_.c_str());
584 28 : if(comm.Get_rank()==0) {
585 : // Only root of group can communicate with other walkers
586 14 : mpi_nw_ = multi_sim_comm.Get_size();
587 14 : mpi_id_ = multi_sim_comm.Get_rank();
588 : }
589 : // Communicate to the other members of the same group
590 : // info abount number of walkers and walker index
591 28 : comm.Bcast(mpi_nw_,0);
592 28 : comm.Bcast(mpi_id_,0);
593 : }
594 : }
595 :
596 240 : for(unsigned i=0; i<doInt_.size(); i++) {
597 60 : if(doInt_[i]) log.printf(" Upper and Lower limits boundaries for the bias of CV %u are activated\n", i);
598 : }
599 30 : if(grid_) {
600 1 : log.printf(" Grid min");
601 8 : for(unsigned i=0; i<gmin.size(); ++i) log.printf(" %s",gmin[i].c_str() );
602 1 : log.printf("\n");
603 1 : log.printf(" Grid max");
604 8 : for(unsigned i=0; i<gmax.size(); ++i) log.printf(" %s",gmax[i].c_str() );
605 1 : log.printf("\n");
606 1 : log.printf(" Grid bin");
607 8 : for(unsigned i=0; i<gbin.size(); ++i) log.printf(" %u",gbin[i]);
608 1 : log.printf("\n");
609 1 : if(spline) {log.printf(" Grid uses spline interpolation\n");}
610 1 : if(sparsegrid) {log.printf(" Grid uses sparse grid\n");}
611 1 : if(wgridstride_>0) {
612 8 : for(unsigned i=0; i<gridfilenames_.size(); ++i) {
613 4 : log.printf(" Grid is written on file %s with stride %d\n",gridfilenames_[i].c_str(),wgridstride_);
614 : }
615 : }
616 1 : if(gridreadfilenames_.size()>0) {
617 0 : for(unsigned i=0; i<gridreadfilenames_.size(); ++i) {
618 0 : log.printf(" Reading bias from grid in file %s \n",gridreadfilenames_[i].c_str());
619 : }
620 : }
621 : }
622 :
623 : // initializing vector of hills
624 30 : hills_.resize(getNumberOfArguments());
625 :
626 : // restart from external grid
627 : bool restartedFromGrid=false;
628 :
629 : // initializing and checking grid
630 30 : if(grid_) {
631 : // check for mesh and sigma size
632 5 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
633 : double a,b;
634 4 : Tools::convert(gmin[i],a);
635 2 : Tools::convert(gmax[i],b);
636 4 : double mesh=(b-a)/((double)gbin[i]);
637 2 : if(adaptive_==FlexibleBin::none) {
638 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";
639 : } else {
640 0 : if(mesh>0.5*sigma0min_[i]||sigma0min_[i]<0.) log<<" WARNING: to use a METAD with a GRID and ADAPTIVE you need to set a Grid Spacing larger than half of the Gaussians \n";
641 : }
642 : }
643 1 : std::string funcl=getLabel() + ".bias";
644 5 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
645 2 : std::vector<Value*> args(1);
646 2 : args[0] = getPntrToArgument(i);
647 4 : vector<std::string> gmin_t(1);
648 4 : vector<std::string> gmax_t(1);
649 2 : vector<unsigned> gbin_t(1);
650 : gmin_t[0] = gmin[i];
651 : gmax_t[0] = gmax[i];
652 2 : gbin_t[0] = gbin[i];
653 : Grid* BiasGrid_;
654 : // Read grid from file
655 2 : if(gridreadfilenames_.size()>0) {
656 0 : IFile gridfile;
657 0 : gridfile.link(*this);
658 0 : if(gridfile.FileExist(gridreadfilenames_[i])) {
659 0 : gridfile.open(gridreadfilenames_[i]);
660 : } else {
661 0 : error("The GRID file you want to read: " + gridreadfilenames_[i] + ", cannot be found!");
662 : }
663 0 : string funcl = getLabel() + ".bias";
664 0 : BiasGrid_ = Grid::create(funcl, args, gridfile, gmin_t, gmax_t, gbin_t, sparsegrid, spline, true);
665 0 : gridfile.close();
666 0 : if(BiasGrid_->getDimension() != args.size()) {
667 0 : error("mismatch between dimensionality of input grid and number of arguments");
668 : }
669 0 : if(getPntrToArgument(i)->isPeriodic() != BiasGrid_->getIsPeriodic()[0]) {
670 0 : error("periodicity mismatch between arguments and input bias");
671 : }
672 0 : log.printf(" Restarting from %s:",gridreadfilenames_[i].c_str());
673 0 : if(getRestart()) restartedFromGrid=true;
674 : } else {
675 2 : if(!sparsegrid) {BiasGrid_=new Grid(funcl,args,gmin_t,gmax_t,gbin_t,spline,true);}
676 0 : else {BiasGrid_=new SparseGrid(funcl,args,gmin_t,gmax_t,gbin_t,spline,true);}
677 4 : std::vector<std::string> actualmin=BiasGrid_->getMin();
678 4 : std::vector<std::string> actualmax=BiasGrid_->getMax();
679 : std::string is;
680 2 : Tools::convert(i,is);
681 2 : if(gmin_t[0]!=actualmin[0]) error("GRID_MIN["+is+"] must be adjusted to "+actualmin[0]+" to fit periodicity");
682 2 : if(gmax_t[0]!=actualmax[0]) error("GRID_MAX["+is+"] must be adjusted to "+actualmax[0]+" to fit periodicity");
683 : }
684 2 : BiasGrids_.push_back(BiasGrid_);
685 : }
686 : }
687 :
688 :
689 : // creating vector of ifile* for hills reading
690 : // open all files at the beginning and read Gaussians if restarting
691 90 : for(int j=0; j<mw_n_; ++j) {
692 240 : for(unsigned i=0; i<hillsfname.size(); ++i) {
693 60 : unsigned k=j*hillsfname.size()+i;
694 : string fname;
695 60 : if(mw_dir_!="") {
696 0 : if(mw_n_>1) {
697 0 : stringstream out; out << j;
698 0 : fname = mw_dir_+"/"+hillsfname[i]+"."+out.str();
699 0 : } else if(walkers_mpi) {
700 0 : fname = mw_dir_+"/"+hillsfname[i];
701 : } else {
702 : fname = hillsfname[i];
703 : }
704 : } else {
705 60 : if(mw_n_>1) {
706 0 : stringstream out; out << j;
707 0 : fname = hillsfname[i]+"."+out.str();
708 : } else {
709 : fname = hillsfname[i];
710 : }
711 : }
712 60 : IFile *ifile = new IFile();
713 60 : ifile->link(*this);
714 60 : ifiles.push_back(ifile);
715 60 : ifilesnames.push_back(fname);
716 60 : if(ifile->FileExist(fname)) {
717 16 : ifile->open(fname);
718 16 : if(getRestart()&&!restartedFromGrid) {
719 0 : log.printf(" Restarting from %s:",ifilesnames[k].c_str());
720 0 : readGaussians(i,ifiles[k]);
721 : }
722 32 : ifiles[k]->reset(false);
723 : // close only the walker own hills file for later writing
724 32 : if(j==mw_id_) ifiles[k]->close();
725 : } else {
726 : // in case a file does not exist and we are restarting, complain that the file was not found
727 44 : if(getRestart()) log<<" WARNING: restart file "<<fname<<" not found\n";
728 : }
729 : }
730 : }
731 :
732 30 : comm.Barrier();
733 30 : if(comm.Get_rank()==0 && walkers_mpi) multi_sim_comm.Barrier();
734 :
735 : // open hills files for writing
736 240 : for(unsigned i=0; i<hillsfname.size(); ++i) {
737 60 : OFile *ofile = new OFile();
738 60 : ofile->link(*this);
739 : // if MPI multiple walkers, only rank 0 will write to file
740 60 : if(walkers_mpi) {
741 56 : int r=0;
742 56 : if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
743 56 : comm.Bcast(r,0);
744 84 : if(r>0) ifilesnames[mw_id_*hillsfname.size()+i]="/dev/null";
745 112 : ofile->enforceSuffix("");
746 : }
747 60 : if(mw_n_>1) ofile->enforceSuffix("");
748 180 : ofile->open(ifilesnames[mw_id_*hillsfname.size()+i]);
749 60 : if(fmt.length()>0) ofile->fmtField(fmt);
750 120 : ofile->addConstantField("multivariate");
751 60 : if(doInt_[i]) {
752 0 : ofile->addConstantField("lower_int").printField("lower_int",lowI_[i]);
753 0 : ofile->addConstantField("upper_int").printField("upper_int",uppI_[i]);
754 : }
755 60 : ofile->setHeavyFlush();
756 : // output periodicities of variables
757 60 : ofile->setupPrintValue( getPntrToArgument(i) );
758 : // push back
759 60 : hillsOfiles_.push_back(ofile);
760 : }
761 :
762 : // Dump grid to files
763 30 : if(wgridstride_ > 0) {
764 8 : for(unsigned i = 0; i < gridfilenames_.size(); ++i) {
765 2 : OFile *ofile = new OFile();
766 2 : ofile->link(*this);
767 : string gridfname_tmp = gridfilenames_[i];
768 2 : if(walkers_mpi) {
769 0 : int r = 0;
770 0 : if(comm.Get_rank() == 0) {
771 0 : r = multi_sim_comm.Get_rank();
772 : }
773 0 : comm.Bcast(r, 0);
774 0 : if(r>0) {
775 : gridfname_tmp = "/dev/null";
776 : }
777 0 : ofile->enforceSuffix("");
778 : }
779 2 : if(mw_n_>1) ofile->enforceSuffix("");
780 2 : ofile->open(gridfname_tmp);
781 2 : ofile->setHeavyFlush();
782 2 : gridfiles_.push_back(ofile);
783 : }
784 : }
785 :
786 90 : log<<" Bibliography "<<plumed.cite("Pfaendtner and Bonomi. J. Chem. Theory Comput. 11, 5062 (2015)");
787 30 : if(doInt_[0]) log<<plumed.cite(
788 0 : "Baftizadeh, Cossio, Pietrucci, and Laio, Curr. Phys. Chem. 2, 79 (2012)");
789 114 : if(mw_n_>1||walkers_mpi) log<<plumed.cite(
790 28 : "Raiteri, Laio, Gervasio, Micheletti, and Parrinello, J. Phys. Chem. B 110, 3533 (2006)");
791 30 : if(adaptive_!=FlexibleBin::none) log<<plumed.cite(
792 0 : "Branduardi, Bussi, and Parrinello, J. Chem. Theory Comput. 8, 2247 (2012)");
793 30 : log<<"\n";
794 30 : }
795 :
796 0 : void PBMetaD::readGaussians(unsigned iarg, IFile *ifile)
797 : {
798 0 : vector<double> center(1);
799 0 : vector<double> sigma(1);
800 : double height;
801 : int nhills=0;
802 0 : bool multivariate=false;
803 :
804 0 : std::vector<Value> tmpvalues;
805 0 : tmpvalues.push_back( Value( this, getPntrToArgument(iarg)->getName(), false ) );
806 :
807 0 : while(scanOneHill(iarg,ifile,tmpvalues,center,sigma,height,multivariate)) {
808 : ;
809 0 : nhills++;
810 0 : if(welltemp_) {height*=(biasf_-1.0)/biasf_;}
811 0 : addGaussian(iarg, Gaussian(center,sigma,height,multivariate));
812 : }
813 0 : log.printf(" %d Gaussians read\n",nhills);
814 0 : }
815 :
816 0 : bool PBMetaD::readChunkOfGaussians(unsigned iarg, IFile *ifile, unsigned n)
817 : {
818 0 : vector<double> center(1);
819 0 : vector<double> sigma(1);
820 : double height;
821 : unsigned nhills=0;
822 0 : bool multivariate=false;
823 0 : std::vector<Value> tmpvalues;
824 0 : tmpvalues.push_back( Value( this, getPntrToArgument(iarg)->getName(), false ) );
825 :
826 0 : while(scanOneHill(iarg,ifile,tmpvalues,center,sigma,height,multivariate)) {
827 : ;
828 0 : if(welltemp_) {height*=(biasf_-1.0)/biasf_;}
829 0 : addGaussian(iarg, Gaussian(center,sigma,height,multivariate));
830 0 : if(nhills==n) {
831 0 : log.printf(" %u Gaussians read\n",nhills);
832 : return true;
833 : }
834 0 : nhills++;
835 : }
836 0 : log.printf(" %u Gaussians read\n",nhills);
837 : return false;
838 : }
839 :
840 896 : void PBMetaD::writeGaussian(unsigned iarg, const Gaussian& hill, OFile *ofile)
841 : {
842 1792 : ofile->printField("time",getTimeStep()*getStep());
843 896 : ofile->printField(getPntrToArgument(iarg),hill.center[0]);
844 :
845 896 : if(hill.multivariate) {
846 0 : ofile->printField("multivariate","true");
847 0 : double lower = sqrt(1./hill.sigma[0]);
848 0 : ofile->printField("sigma_"+getPntrToArgument(iarg)->getName()+"_"+
849 : getPntrToArgument(iarg)->getName(),lower);
850 : } else {
851 2688 : ofile->printField("multivariate","false");
852 2688 : ofile->printField("sigma_"+getPntrToArgument(iarg)->getName(),hill.sigma[0]);
853 : }
854 896 : double height=hill.height;
855 896 : if(welltemp_) height *= biasf_/(biasf_-1.0);
856 1792 : ofile->printField("height",height);
857 1792 : ofile->printField("biasf",biasf_);
858 896 : if(mw_n_>1) ofile->printField("clock",int(std::time(0)));
859 896 : ofile->printField();
860 896 : }
861 :
862 896 : void PBMetaD::addGaussian(unsigned iarg, const Gaussian& hill)
863 : {
864 1784 : if(!grid_) {hills_[iarg].push_back(hill);}
865 : else {
866 8 : vector<unsigned> nneighb=getGaussianSupport(iarg, hill);
867 16 : vector<Grid::index_t> neighbors=BiasGrids_[iarg]->getNeighbors(hill.center,nneighb);
868 8 : vector<double> der(1);
869 8 : vector<double> xx(1);
870 8 : if(comm.Get_size()==1) {
871 1768 : for(unsigned i=0; i<neighbors.size(); ++i) {
872 584 : Grid::index_t ineigh=neighbors[i];
873 584 : der[0]=0.0;
874 584 : BiasGrids_[iarg]->getPoint(ineigh,xx);
875 584 : double bias=evaluateGaussian(iarg,xx,hill,&der[0]);
876 584 : BiasGrids_[iarg]->addValueAndDerivatives(ineigh,bias,der);
877 : }
878 : } else {
879 0 : unsigned stride=comm.Get_size();
880 0 : unsigned rank=comm.Get_rank();
881 0 : vector<double> allder(neighbors.size(),0.0);
882 0 : vector<double> allbias(neighbors.size(),0.0);
883 0 : for(unsigned i=rank; i<neighbors.size(); i+=stride) {
884 0 : Grid::index_t ineigh=neighbors[i];
885 0 : BiasGrids_[iarg]->getPoint(ineigh,xx);
886 0 : allbias[i]=evaluateGaussian(iarg,xx,hill,&allder[i]);
887 : }
888 0 : comm.Sum(allbias);
889 0 : comm.Sum(allder);
890 0 : for(unsigned i=0; i<neighbors.size(); ++i) {
891 0 : Grid::index_t ineigh=neighbors[i];
892 0 : der[0]=allder[i];
893 0 : BiasGrids_[iarg]->addValueAndDerivatives(ineigh,allbias[i],der);
894 : }
895 : }
896 : }
897 896 : }
898 :
899 8 : vector<unsigned> PBMetaD::getGaussianSupport(unsigned iarg, const Gaussian& hill)
900 : {
901 : vector<unsigned> nneigh;
902 : double cutoff;
903 8 : if(hill.multivariate) {
904 0 : double maxautoval=1./hill.sigma[0];
905 0 : cutoff=sqrt(2.0*DP2CUTOFF*maxautoval);
906 : } else {
907 8 : cutoff=sqrt(2.0*DP2CUTOFF)*hill.sigma[0];
908 : }
909 :
910 16 : if(doInt_[iarg]) {
911 0 : if(hill.center[0]+cutoff > uppI_[iarg] || hill.center[0]-cutoff < lowI_[iarg]) {
912 : // in this case, we updated the entire grid to avoid problems
913 0 : return BiasGrids_[iarg]->getNbin();
914 : } else {
915 0 : nneigh.push_back( static_cast<unsigned>(ceil(cutoff/BiasGrids_[iarg]->getDx()[0])));
916 : return nneigh;
917 : }
918 : }
919 :
920 32 : nneigh.push_back( static_cast<unsigned>(ceil(cutoff/BiasGrids_[iarg]->getDx()[0])) );
921 :
922 : return nneigh;
923 : }
924 :
925 972 : double PBMetaD::getBiasAndDerivatives(unsigned iarg, const vector<double>& cv, double* der)
926 : {
927 972 : double bias=0.0;
928 972 : if(!grid_) {
929 954 : unsigned stride=comm.Get_size();
930 954 : unsigned rank=comm.Get_rank();
931 15180 : for(unsigned i=rank; i<hills_[iarg].size(); i+=stride) {
932 4424 : bias += evaluateGaussian(iarg,cv,hills_[iarg][i],der);
933 : }
934 954 : comm.Sum(bias);
935 954 : if(der) comm.Sum(der,1);
936 : } else {
937 18 : if(der) {
938 10 : vector<double> vder(1);
939 20 : bias = BiasGrids_[iarg]->getValueAndDerivatives(cv,vder);
940 10 : der[0] = vder[0];
941 : } else {
942 16 : bias = BiasGrids_[iarg]->getValue(cv);
943 : }
944 : }
945 :
946 972 : return bias;
947 : }
948 :
949 5008 : double PBMetaD::evaluateGaussian(unsigned iarg, const vector<double>& cv, const Gaussian& hill, double* der)
950 : {
951 : double bias=0.0;
952 : // I use a pointer here because cv is const (and should be const)
953 : // but when using doInt it is easier to locally replace cv[0] with
954 : // the upper/lower limit in case it is out of range
955 : const double *pcv=NULL;
956 : double tmpcv[1]; // tmp array with cv (to be used with doInt_)
957 5008 : tmpcv[0]=cv[0];
958 : bool isOutOfInt = false;
959 10016 : if(doInt_[iarg]) {
960 0 : if(cv[0]<lowI_[iarg]) { tmpcv[0]=lowI_[iarg]; isOutOfInt = true; }
961 0 : else if(cv[0]>uppI_[iarg]) { tmpcv[0]=uppI_[iarg]; isOutOfInt = true; }
962 : }
963 : pcv=&(tmpcv[0]);
964 :
965 5008 : if(hill.multivariate) {
966 0 : double dp = difference(iarg, hill.center[0], pcv[0]);
967 0 : double dp2 = 0.5 * dp * dp * hill.sigma[0];
968 0 : if(dp2<DP2CUTOFF) {
969 0 : bias = hill.height*exp(-dp2);
970 0 : if(der && !isOutOfInt) {
971 0 : der[0] += -bias * dp * hill.sigma[0];
972 : }
973 : }
974 : } else {
975 10016 : double dp = difference(iarg, hill.center[0], pcv[0]) * hill.invsigma[0];
976 5008 : double dp2 = 0.5 * dp * dp;
977 5008 : if(dp2<DP2CUTOFF) {
978 4992 : bias = hill.height*exp(-dp2);
979 4992 : if(der && !isOutOfInt) {
980 5560 : der[0] += -bias * dp * hill.invsigma[0];
981 : }
982 : }
983 : }
984 :
985 5008 : return bias;
986 : }
987 :
988 258 : void PBMetaD::calculate()
989 : {
990 : // this is because presently there is no way to properly pass information
991 : // on adaptive hills (diff) after exchanges:
992 258 : if(adaptive_==FlexibleBin::diffusion && getExchangeStep()) error("ADAPTIVE=DIFF is not compatible with replica exchange");
993 :
994 258 : vector<double> cv(1);
995 : double der[1];
996 258 : vector<double> bias(getNumberOfArguments());
997 258 : vector<double> deriv(getNumberOfArguments());
998 :
999 258 : double ncv = (double) getNumberOfArguments();
1000 : double bmin = 1.0e+19;
1001 1290 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1002 516 : cv[0] = getArgument(i);
1003 516 : der[0] = 0.0;
1004 516 : bias[i] = getBiasAndDerivatives(i, cv, der);
1005 516 : deriv[i] = der[0];
1006 516 : if(bias[i] < bmin) bmin = bias[i];
1007 : }
1008 : double ene = 0.;
1009 1290 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1010 1032 : ene += exp((-bias[i]+bmin)/kbt_);
1011 : }
1012 :
1013 : // set Forces - set them to zero if SELECTOR is active
1014 258 : if(do_select_) current_value_ = static_cast<unsigned>(plumed.passMap[selector_]);
1015 :
1016 258 : if(!do_select_ || (do_select_ && select_value_==current_value_)) {
1017 1290 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1018 1548 : const double f = - exp((-bias[i]+bmin)/kbt_) / (ene) * deriv[i];
1019 516 : setOutputForce(i, f);
1020 : }
1021 : }
1022 :
1023 258 : if(do_select_ && select_value_!=current_value_) {
1024 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) setOutputForce(i, 0.0);
1025 : }
1026 :
1027 : // set bias
1028 258 : ene = -kbt_ * (std::log(ene) - std::log(ncv)) + bmin;
1029 : setBias(ene);
1030 258 : }
1031 :
1032 258 : void PBMetaD::update()
1033 : {
1034 : bool multivariate;
1035 : // adding hills criteria
1036 : bool nowAddAHill;
1037 258 : if(getStep()%stride_==0 && !isFirstStep) nowAddAHill=true;
1038 : else {
1039 : nowAddAHill=false;
1040 30 : isFirstStep=false;
1041 : }
1042 :
1043 : // if you use adaptive, call the FlexibleBin
1044 258 : if(adaptive_!=FlexibleBin::none) {
1045 0 : for(unsigned i=0; i<getNumberOfArguments(); i++) flexbin[i].update(nowAddAHill,i);
1046 : multivariate=true;
1047 : } else {
1048 : multivariate=false;
1049 : }
1050 :
1051 258 : if(nowAddAHill && (!do_select_ || (do_select_ && select_value_==current_value_))) {
1052 : // get all biases and heights
1053 228 : vector<double> cv(getNumberOfArguments());
1054 228 : vector<double> bias(getNumberOfArguments());
1055 228 : vector<double> thissigma(getNumberOfArguments());
1056 228 : vector<double> height(getNumberOfArguments());
1057 228 : vector<double> cv_tmp(1);
1058 228 : vector<double> sigma_tmp(1);
1059 : double norm = 0.0;
1060 : double bmin = 1.0e+19;
1061 1140 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1062 456 : if(adaptive_!=FlexibleBin::none) thissigma[i]=flexbin[i].getInverseMatrix(i)[0];
1063 912 : else thissigma[i]=sigma0_[i];
1064 912 : cv[i] = getArgument(i);
1065 456 : cv_tmp[0] = getArgument(i);
1066 456 : bias[i] = getBiasAndDerivatives(i, cv_tmp);
1067 456 : if(bias[i] < bmin) bmin = bias[i];
1068 : }
1069 : // calculate heights and norm
1070 1140 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1071 912 : double h = exp((-bias[i]+bmin)/kbt_);
1072 456 : norm += h;
1073 456 : height[i] = h;
1074 : }
1075 : // normalize and apply welltemp correction
1076 1140 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1077 912 : height[i] *= height0_ / norm;
1078 1352 : if(welltemp_) height[i] *= exp(-bias[i]/(kbt_*(biasf_-1.0)));
1079 : }
1080 :
1081 : // MPI Multiple walkers: share hills and add them all
1082 228 : if(walkers_mpi) {
1083 : // Allocate arrays to store all walkers hills
1084 440 : std::vector<double> all_cv(mpi_nw_*cv.size(), 0.0);
1085 220 : std::vector<double> all_sigma(mpi_nw_*getNumberOfArguments(), 0.0);
1086 440 : std::vector<double> all_height(mpi_nw_*height.size(), 0.0);
1087 220 : if(comm.Get_rank()==0) {
1088 : // fill in value
1089 550 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1090 220 : unsigned j = mpi_id_ * getNumberOfArguments() + i;
1091 660 : all_cv[j] = cv[i];
1092 220 : all_sigma[j] = thissigma[i];
1093 220 : all_height[j] = height[i];
1094 : }
1095 : // Communicate (only root)
1096 220 : multi_sim_comm.Sum(&all_cv[0], all_cv.size());
1097 220 : multi_sim_comm.Sum(&all_sigma[0], all_sigma.size());
1098 220 : multi_sim_comm.Sum(&all_height[0], all_height.size());
1099 : }
1100 : // Share info with group members
1101 440 : comm.Sum(&all_cv[0], all_cv.size());
1102 440 : comm.Sum(&all_sigma[0], all_sigma.size());
1103 440 : comm.Sum(&all_height[0], all_height.size());
1104 : // now add hills one by one
1105 1100 : for(unsigned j=0; j<mpi_nw_; ++j) {
1106 2200 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1107 2640 : cv_tmp[0] = all_cv[j*cv.size()+i];
1108 1760 : double height_tmp = all_height[j*cv.size()+i];
1109 880 : sigma_tmp[0] = all_sigma[j*cv.size()+i];
1110 1760 : Gaussian newhill = Gaussian(cv_tmp, sigma_tmp, height_tmp, multivariate);
1111 880 : addGaussian(i, newhill);
1112 880 : writeGaussian(i, newhill, hillsOfiles_[i]);
1113 : }
1114 : }
1115 : // just add your own hills
1116 : } else {
1117 40 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1118 32 : cv_tmp[0] = cv[i];
1119 16 : if(adaptive_!=FlexibleBin::none) sigma_tmp[0]=thissigma[i];
1120 16 : else sigma_tmp[0] = sigma0_[i];
1121 48 : Gaussian newhill = Gaussian(cv_tmp, sigma_tmp, height[i], multivariate);
1122 16 : addGaussian(i, newhill);
1123 16 : writeGaussian(i, newhill, hillsOfiles_[i]);
1124 : }
1125 : }
1126 : }
1127 :
1128 : // write grid files
1129 258 : if(wgridstride_>0 && (getStep()%wgridstride_==0 || getCPT())) {
1130 5 : int r = 0;
1131 5 : if(walkers_mpi) {
1132 0 : if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
1133 0 : comm.Bcast(r,0);
1134 : }
1135 5 : if(r==0) {
1136 40 : for(unsigned i=0; i<gridfiles_.size(); ++i) {
1137 10 : gridfiles_[i]->rewind();
1138 20 : BiasGrids_[i]->writeToFile(*gridfiles_[i]);
1139 10 : gridfiles_[i]->flush();
1140 : }
1141 : }
1142 : }
1143 :
1144 : // if multiple walkers and time to read Gaussians
1145 258 : if(mw_n_>1 && getStep()%mw_rstride_==0) {
1146 0 : for(int j=0; j<mw_n_; ++j) {
1147 0 : for(unsigned i=0; i<hillsfname.size(); ++i) {
1148 0 : unsigned k=j*hillsfname.size()+i;
1149 : // don't read your own Gaussians
1150 0 : if(j==mw_id_) continue;
1151 : // if the file is not open yet
1152 0 : if(!(ifiles[k]->isOpen())) {
1153 : // check if it exists now and open it!
1154 0 : if(ifiles[k]->FileExist(ifilesnames[k])) {
1155 0 : ifiles[k]->open(ifilesnames[k]);
1156 0 : ifiles[k]->reset(false);
1157 : }
1158 : // otherwise read the new Gaussians
1159 : } else {
1160 0 : log.printf(" Reading hills from %s:",ifilesnames[k].c_str());
1161 0 : readGaussians(i,ifiles[k]);
1162 0 : ifiles[k]->reset(false);
1163 : }
1164 : }
1165 : }
1166 : }
1167 :
1168 258 : }
1169 :
1170 : /// takes a pointer to the file and a template string with values v and gives back the next center, sigma and height
1171 0 : bool PBMetaD::scanOneHill(unsigned iarg, IFile *ifile, vector<Value> &tmpvalues, vector<double> ¢er, vector<double> &sigma, double &height, bool &multivariate)
1172 : {
1173 : double dummy;
1174 0 : multivariate=false;
1175 0 : if(ifile->scanField("time",dummy)) {
1176 0 : ifile->scanField( &tmpvalues[0] );
1177 0 : if( tmpvalues[0].isPeriodic() && ! getPntrToArgument(iarg)->isPeriodic() ) {
1178 0 : error("in hills file periodicity for variable " + tmpvalues[0].getName() + " does not match periodicity in input");
1179 0 : } else if( tmpvalues[0].isPeriodic() ) {
1180 0 : std::string imin, imax; tmpvalues[0].getDomain( imin, imax );
1181 0 : std::string rmin, rmax; getPntrToArgument(iarg)->getDomain( rmin, rmax );
1182 0 : if( imin!=rmin || imax!=rmax ) {
1183 0 : error("in hills file periodicity for variable " + tmpvalues[0].getName() + " does not match periodicity in input");
1184 : }
1185 : }
1186 0 : center[0]=tmpvalues[0].get();
1187 : std::string sss;
1188 0 : ifile->scanField("multivariate",sss);
1189 0 : if(sss=="true") multivariate=true;
1190 0 : else if(sss=="false") multivariate=false;
1191 0 : else plumed_merror("cannot parse multivariate = "+ sss);
1192 :
1193 0 : if(multivariate) {
1194 0 : ifile->scanField("sigma_"+getPntrToArgument(iarg)->getName()+"_"+
1195 : getPntrToArgument(iarg)->getName(),sigma[0]);
1196 0 : sigma[0] = 1./(sigma[0]*sigma[0]);
1197 : } else {
1198 0 : ifile->scanField("sigma_"+getPntrToArgument(iarg)->getName(),sigma[0]);
1199 : }
1200 0 : ifile->scanField("height",height);
1201 0 : ifile->scanField("biasf",dummy);
1202 0 : if(ifile->FieldExist("clock")) ifile->scanField("clock",dummy);
1203 0 : if(ifile->FieldExist("lower_int")) ifile->scanField("lower_int",dummy);
1204 0 : if(ifile->FieldExist("upper_int")) ifile->scanField("upper_int",dummy);
1205 0 : ifile->scanField();
1206 : return true;
1207 : } else {
1208 : return false;
1209 : }
1210 :
1211 : }
1212 :
1213 258 : bool PBMetaD::checkNeedsGradients()const
1214 : {
1215 258 : if(adaptive_==FlexibleBin::geometry) {
1216 0 : if(getStep()%stride_==0 && !isFirstStep) return true;
1217 : else return false;
1218 : } else return false;
1219 : }
1220 :
1221 : }
1222 4839 : }
1223 :
|
