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