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