Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 METAD
49 : /*
50 : Used to performed metadynamics on one or more collective variables.
51 :
52 : In a metadynamics simulations a history dependent bias composed of
53 : intermittently added Gaussian functions is added to the potential \cite metad.
54 :
55 : \f[
56 : V(\vec{s},t) = \sum_{ k \tau < t} W(k \tau)
57 : \exp\left(
58 : -\sum_{i=1}^{d} \frac{(s_i-s_i^{(0)}(k \tau))^2}{2\sigma_i^2}
59 : \right).
60 : \f]
61 :
62 : This potential forces the system away from the kinetic traps in the potential energy surface
63 : and out into the unexplored parts of the energy landscape. Information on the Gaussian
64 : functions from which this potential is composed is output to a file called HILLS, which
65 : is used both the restart the calculation and to reconstruct the free energy as a function of the CVs.
66 : The free energy can be reconstructed from a metadynamics calculation because the final bias is given
67 : by:
68 :
69 : \f[
70 : V(\vec{s}) = -F(\vec{s})
71 : \f]
72 :
73 : During post processing the free energy can be calculated in this way using the \ref sum_hills
74 : utility.
75 :
76 : In the simplest possible implementation of a metadynamics calculation the expense of a metadynamics
77 : calculation increases with the length of the simulation as one has to, at every step, evaluate
78 : the values of a larger and larger number of Gaussian kernels. To avoid this issue you can
79 : store the bias on a grid. This approach is similar to that proposed in \cite babi08jcp but has the
80 : advantage that the grid spacing is independent on the Gaussian width.
81 : Notice that you should
82 : provide either the number of bins for every collective variable (GRID_BIN) or
83 : the desired grid spacing (GRID_SPACING). In case you provide both PLUMED will use
84 : the most conservative choice (highest number of bins) for each dimension.
85 : In case you do not provide any information about bin size (neither GRID_BIN nor GRID_SPACING)
86 : and if Gaussian width is fixed PLUMED will use 1/5 of the Gaussian width as grid spacing.
87 : This default choice should be reasonable for most applications.
88 :
89 : Metadynamics can be restarted either from a HILLS file as well as from a GRID, in this second
90 : case one can first save a GRID using GRID_WFILE (and GRID_WSTRIDE) and at a later stage read
91 : it using GRID_RFILE.
92 :
93 : Another option that is available in plumed is well-tempered metadynamics \cite Barducci:2008. In this
94 : variant of metadynamics the heights of the Gaussian hills are scaled at each step so the bias is now
95 : given by:
96 :
97 : \f[
98 : V({s},t)= \sum_{t'=0,\tau_G,2\tau_G,\dots}^{t'<t} W e^{-V({s}({q}(t'),t')/\Delta T} \exp\left(
99 : -\sum_{i=1}^{d} \frac{(s_i({q})-s_i({q}(t'))^2}{2\sigma_i^2}
100 : \right),
101 : \f]
102 :
103 : This method ensures that the bias converges more smoothly. It should be noted that, in the case of well-tempered metadynamics, in
104 : the output printed the Gaussian height is re-scaled using the bias factor.
105 : Also notice that with well-tempered metadynamics the HILLS file does not contain the bias,
106 : but the negative of the free-energy estimate. This choice has the advantage that
107 : one can restart a simulation using a different value for the \f$\Delta T\f$. The applied bias will be scaled accordingly.
108 :
109 : Note that you can use here also the flexible Gaussian approach \cite Branduardi:2012dl
110 : in which you can adapt the Gaussian to the extent of Cartesian space covered by a variable or
111 : to the space in collective variable covered in a given time. In this case the width of the deposited
112 : Gaussian potential is denoted by one value only that is a Cartesian space (ADAPTIVE=GEOM) or a time
113 : (ADAPTIVE=DIFF). Note that a specific integration technique for the deposited Gaussian kernels
114 : should be used in this case. Check the documentation for utility sum_hills.
115 :
116 : With the keyword INTERVAL one changes the metadynamics algorithm setting the bias force equal to zero
117 : outside boundary \cite baftizadeh2012protein. If, for example, metadynamics is performed on a CV s and one is interested only
118 : to the free energy for s > boundary, the history dependent potential is still updated according to the above
119 : equations but the metadynamics force is set to zero for s < boundary. Notice that Gaussian kernels are added also
120 : if s < boundary, as the tails of these Gaussian kernels influence VG in the relevant region s > boundary. In this way, the
121 : force on the system in the region s > boundary comes from both metadynamics and the force field, in the region
122 : s < boundary only from the latter. This approach allows obtaining a history-dependent bias potential VG that
123 : fluctuates around a stable estimator, equal to the negative of the free energy far enough from the
124 : boundaries. Note that:
125 : - It works only for one-dimensional biases;
126 : - It works both with and without GRID;
127 : - The interval limit boundary in a region where the free energy derivative is not large;
128 : - If in the region outside the limit boundary the system has a free energy minimum, the INTERVAL keyword should
129 : be used together with a \ref UPPER_WALLS or \ref LOWER_WALLS at boundary.
130 :
131 : As a final note, since version 2.0.2 when the system is outside of the selected interval the force
132 : is set to zero and the bias value to the value at the corresponding boundary. This allows acceptances
133 : for replica exchange methods to be computed correctly.
134 :
135 : Multiple walkers \cite multiplewalkers can also be used. See below the examples.
136 :
137 :
138 : The \f$c(t)\f$ reweighting factor can also be calculated on the fly using the equations
139 : presented in \cite Tiwary_jp504920s.
140 : The expression used to calculate \f$c(t)\f$ follows directly from Eq. 3 in \cite Tiwary_jp504920s,
141 : where \f$F(\vec{s})=-\gamma/(\gamma-1) V(\vec{s})\f$.
142 : This gives smoother results than equivalent Eqs. 13 and Eqs. 14 in that paper.
143 : The \f$c(t)\f$ is given by the rct component while the bias
144 : normalized by \f$c(t)\f$ is given by the rbias component (rbias=bias-rct) which can be used
145 : to obtain a reweighted histogram.
146 : The calculation of \f$c(t)\f$ is enabled by using the keyword CALC_RCT.
147 : By default \f$c(t)\f$ is updated every time the bias changes, but if this slows down the simulation
148 : the keyword RCT_USTRIDE can be set to a value higher than 1.
149 : This option requires that a grid is used.
150 :
151 : Additional material and examples can be also found in the tutorials:
152 :
153 : - \ref lugano-3
154 :
155 : Concurrent metadynamics
156 : as done e.g. in Ref. \cite gil2015enhanced . This indeed can be obtained by using the METAD
157 : action multiple times in the same input file.
158 :
159 : \par Examples
160 :
161 : The following input is for a standard metadynamics calculation using as
162 : collective variables the distance between atoms 3 and 5
163 : and the distance between atoms 2 and 4. The value of the CVs and
164 : the metadynamics bias potential are written to the COLVAR file every 100 steps.
165 : \plumedfile
166 : DISTANCE ATOMS=3,5 LABEL=d1
167 : DISTANCE ATOMS=2,4 LABEL=d2
168 : METAD ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3 PACE=500 LABEL=restraint
169 : PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
170 : \endplumedfile
171 : (See also \ref DISTANCE \ref PRINT).
172 :
173 : \par
174 : If you use adaptive Gaussian kernels, with diffusion scheme where you use
175 : a Gaussian that should cover the space of 20 time steps in collective variables.
176 : Note that in this case the histogram correction is needed when summing up hills.
177 : \plumedfile
178 : DISTANCE ATOMS=3,5 LABEL=d1
179 : DISTANCE ATOMS=2,4 LABEL=d2
180 : METAD ARG=d1,d2 SIGMA=20 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=DIFF
181 : PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
182 : \endplumedfile
183 :
184 : \par
185 : If you use adaptive Gaussian kernels, with geometrical scheme where you use
186 : a Gaussian that should cover the space of 0.05 nm in Cartesian space.
187 : Note that in this case the histogram correction is needed when summing up hills.
188 : \plumedfile
189 : DISTANCE ATOMS=3,5 LABEL=d1
190 : DISTANCE ATOMS=2,4 LABEL=d2
191 : METAD ARG=d1,d2 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=GEOM
192 : PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
193 : \endplumedfile
194 :
195 : \par
196 : When using adaptive Gaussian kernels you might want to limit how the hills width can change.
197 : You can use SIGMA_MIN and SIGMA_MAX keywords.
198 : The sigmas should specified in terms of CV so you should use the CV units.
199 : Note that if you use a negative number, this means that the limit is not set.
200 : Note also that in this case the histogram correction is needed when summing up hills.
201 : \plumedfile
202 : DISTANCE ATOMS=3,5 LABEL=d1
203 : DISTANCE ATOMS=2,4 LABEL=d2
204 : METAD ...
205 : ARG=d1,d2 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=GEOM
206 : SIGMA_MIN=0.2,0.1 SIGMA_MAX=0.5,1.0
207 : ... METAD
208 : PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
209 : \endplumedfile
210 :
211 : \par
212 : Multiple walkers can be also use as in \cite multiplewalkers
213 : These are enabled by setting the number of walker used, the id of the
214 : current walker which interprets the input file, the directory where the
215 : hills containing files resides, and the frequency to read the other walkers.
216 : Here is an example
217 : \plumedfile
218 : DISTANCE ATOMS=3,5 LABEL=d1
219 : METAD ...
220 : ARG=d1 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint
221 : WALKERS_N=10
222 : WALKERS_ID=3
223 : WALKERS_DIR=../
224 : WALKERS_RSTRIDE=100
225 : ... METAD
226 : \endplumedfile
227 : where WALKERS_N is the total number of walkers, WALKERS_ID is the
228 : id of the present walker (starting from 0 ) and the WALKERS_DIR is the directory
229 : where all the walkers are located. WALKERS_RSTRIDE is the number of step between
230 : one update and the other. Since version 2.2.5, hills files are automatically
231 : flushed every WALKERS_RSTRIDE steps.
232 :
233 : \par
234 : The \f$c(t)\f$ reweighting factor can be calculated on the fly using the equations
235 : presented in \cite Tiwary_jp504920s as described above.
236 : This is enabled by using the keyword CALC_RCT,
237 : and can be done only if the bias is defined on a grid.
238 : \plumedfile
239 : phi: TORSION ATOMS=1,2,3,4
240 : psi: TORSION ATOMS=5,6,7,8
241 :
242 : METAD ...
243 : LABEL=metad
244 : ARG=phi,psi SIGMA=0.20,0.20 HEIGHT=1.20 BIASFACTOR=5 TEMP=300.0 PACE=500
245 : GRID_MIN=-pi,-pi GRID_MAX=pi,pi GRID_BIN=150,150
246 : CALC_RCT
247 : RCT_USTRIDE=10
248 : ... METAD
249 : \endplumedfile
250 : Here we have asked that the calculation is performed every 10 hills deposition by using
251 : RCT_USTRIDE keyword. If this keyword is not given, the calculation will
252 : by default be performed every time the bias changes. The \f$c(t)\f$ reweighting factor will be given
253 : in the rct component while the instantaneous value of the bias potential
254 : normalized using the \f$c(t)\f$ reweighting factor is given in the rbias component
255 : [rbias=bias-rct] which can be used to obtain a reweighted histogram or
256 : free energy surface using the \ref HISTOGRAM analysis.
257 :
258 : \par
259 : The kinetics of the transitions between basins can also be analyzed on the fly as
260 : in \cite PRL230602. The flag ACCELERATION turn on accumulation of the acceleration
261 : factor that can then be used to determine the rate. This method can be used together
262 : with \ref COMMITTOR analysis to stop the simulation when the system get to the target basin.
263 : It must be used together with Well-Tempered Metadynamics. If restarting from a previous
264 : metadynamics you need to use the ACCELERATION_RFILE keyword to give the name of the
265 : data file from which the previous value of the acceleration factor should be read, otherwise the
266 : calculation of the acceleration factor will be wrong.
267 :
268 : \par
269 : By using the flag FREQUENCY_ADAPTIVE the frequency adaptive scheme introduced in \cite Wang-JCP-2018
270 : is turned on. The frequency for hill addition then changes dynamically based on the acceleration factor
271 : according to the following equation
272 : \f[
273 : \tau_{\mathrm{dep}}(t) =
274 : \min\left[
275 : \tau_0 \cdot
276 : \max\left[\frac{\alpha(t)}{\theta},1\right]
277 : ,\tau_{c}
278 : \right]
279 : \f]
280 : where \f$\tau_0\f$ is the initial hill addition frequency given by the PACE keyword,
281 : \f$\tau_{c}\f$ is the maximum allowed frequency given by the FA_MAX_PACE keyword,
282 : \f$\alpha(t)\f$ is the instantaneous acceleration factor at time \f$t\f$,
283 : and \f$\theta\f$ is a threshold value that acceleration factor has to reach before
284 : triggering a change in the hill addition frequency given by the FA_MIN_ACCELERATION keyword.
285 : The frequency for updating the hill addition frequency according to this equation is
286 : given by the FA_UPDATE_FREQUENCY keyword, by default it is the same as the value given
287 : in PACE. The hill hill addition frequency increase monotonously such that if the
288 : instantaneous acceleration factor is lower than in the previous updating step the
289 : previous \f$\tau_{\mathrm{dep}}\f$ is kept rather than updating it to a lower value.
290 : The instantaneous hill addition frequency \f$\tau_{\mathrm{dep}}(t)\f$ is outputted
291 : to pace component. Note that if restarting from a previous metadynamics run you need to
292 : use the ACCELERATION_RFILE keyword to read in the acceleration factors from the
293 : previous run, otherwise the hill addition frequency will start from the initial
294 : frequency.
295 :
296 :
297 : \par
298 : You can also provide a target distribution using the keyword TARGET
299 : \cite white2015designing
300 : \cite marinelli2015ensemble
301 : \cite gil2016empirical
302 : The TARGET should be a grid containing a free-energy (i.e. the -\f$k_B\f$T*log of the desired target distribution).
303 : Gaussian kernels will then be scaled by a factor
304 : \f[
305 : e^{\beta(\tilde{F}(s)-\tilde{F}_{max})}
306 : \f]
307 : Here \f$\tilde{F}(s)\f$ is the free energy defined on the grid and \f$\tilde{F}_{max}\f$ its maximum value.
308 : Notice that we here used the maximum value as in ref \cite gil2016empirical
309 : This choice allows to avoid exceedingly large Gaussian kernels to be added. However,
310 : it could make the Gaussian too small. You should always choose carefully the HEIGHT parameter
311 : in this case.
312 : The grid file should be similar to other PLUMED grid files in that it should contain
313 : both the target free-energy and its derivatives.
314 :
315 : Notice that if you wish your simulation to converge to the target free energy you should use
316 : the DAMPFACTOR command to provide a global tempering \cite dama2014well
317 : Alternatively, if you use a BIASFACTOR your simulation will converge to a free
318 : energy that is a linear combination of the target free energy and of the intrinsic free energy
319 : determined by the original force field.
320 :
321 : \plumedfile
322 : DISTANCE ATOMS=3,5 LABEL=d1
323 : METAD ...
324 : LABEL=t1
325 : ARG=d1 SIGMA=0.05 TAU=200 DAMPFACTOR=100 PACE=250
326 : GRID_MIN=1.14 GRID_MAX=1.32 GRID_BIN=6
327 : TARGET=dist.grid
328 : ... METAD
329 :
330 : PRINT ARG=d1,t1.bias STRIDE=100 FILE=COLVAR
331 : \endplumedfile
332 :
333 : The file dist.dat for this calculation would read:
334 :
335 : \auxfile{dist.grid}
336 : #! FIELDS d1 t1.target der_d1
337 : #! SET min_d1 1.14
338 : #! SET max_d1 1.32
339 : #! SET nbins_d1 6
340 : #! SET periodic_d1 false
341 : 1.1400 0.0031 0.1101
342 : 1.1700 0.0086 0.2842
343 : 1.2000 0.0222 0.6648
344 : 1.2300 0.0521 1.4068
345 : 1.2600 0.1120 2.6873
346 : 1.2900 0.2199 4.6183
347 : 1.3200 0.3948 7.1055
348 : \endauxfile
349 :
350 : Notice that BIASFACTOR can also be chosen as equal to 1. In this case one will perform
351 : unbiased sampling. Instead of using HEIGHT, one should provide the TAU parameter.
352 : \plumedfile
353 : d: DISTANCE ATOMS=3,5
354 : METAD ARG=d SIGMA=0.1 TAU=4.0 TEMP=300 PACE=100 BIASFACTOR=1.0
355 : \endplumedfile
356 : The HILLS file obtained will still work with `plumed sum_hills` so as to plot a free-energy.
357 : The case where this makes sense is probably that of RECT simulations.
358 :
359 : Regarding RECT simulations, you can also use the RECT keyword so as to avoid using multiple input files.
360 : For instance, a single input file will be
361 : \plumedfile
362 : d: DISTANCE ATOMS=3,5
363 : METAD ARG=d SIGMA=0.1 TAU=4.0 TEMP=300 PACE=100 RECT=1.0,1.5,2.0,3.0
364 : \endplumedfile
365 : The number of elements in the RECT array should be equal to the number of replicas.
366 :
367 :
368 :
369 :
370 :
371 : */
372 : //+ENDPLUMEDOC
373 :
374 966 : class MetaD : public Bias {
375 :
376 : private:
377 65202 : struct Gaussian {
378 : vector<double> center;
379 : vector<double> sigma;
380 : double height;
381 : bool multivariate; // this is required to discriminate the one dimensional case
382 : vector<double> invsigma;
383 5925 : Gaussian(const vector<double> & center,const vector<double> & sigma,double height, bool multivariate ):
384 5925 : center(center),sigma(sigma),height(height),multivariate(multivariate),invsigma(sigma) {
385 : // to avoid troubles from zero element in flexible hills
386 57414 : 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;
387 5925 : }
388 : };
389 288 : struct TemperingSpecs {
390 : bool is_active;
391 : std::string name_stem;
392 : std::string name;
393 : double biasf;
394 : double threshold;
395 : double alpha;
396 144 : inline TemperingSpecs(bool is_active, const std::string &name_stem, const std::string &name, double biasf, double threshold, double alpha) :
397 432 : is_active(is_active), name_stem(name_stem), name(name), biasf(biasf), threshold(threshold), alpha(alpha)
398 144 : {}
399 : };
400 : vector<double> sigma0_;
401 : vector<double> sigma0min_;
402 : vector<double> sigma0max_;
403 : vector<Gaussian> hills_;
404 : OFile hillsOfile_;
405 : OFile gridfile_;
406 : std::unique_ptr<Grid> BiasGrid_;
407 : bool storeOldGrids_;
408 : int wgridstride_;
409 : bool grid_;
410 : double height0_;
411 : double biasf_;
412 : static const size_t n_tempering_options_ = 1;
413 : static const string tempering_names_[1][2];
414 : double dampfactor_;
415 : struct TemperingSpecs tt_specs_;
416 : std::string targetfilename_;
417 : std::unique_ptr<Grid> TargetGrid_;
418 : double kbt_;
419 : int stride_;
420 : bool welltemp_;
421 : //
422 : int current_stride;
423 : bool freq_adaptive_;
424 : int fa_update_frequency_;
425 : int fa_max_stride_;
426 : double fa_min_acceleration_;
427 : //
428 : std::unique_ptr<double[]> dp_;
429 : int adaptive_;
430 : std::unique_ptr<FlexibleBin> flexbin;
431 : int mw_n_;
432 : string mw_dir_;
433 : int mw_id_;
434 : int mw_rstride_;
435 : bool walkers_mpi;
436 : unsigned mpi_nw_;
437 : unsigned mpi_mw_;
438 : bool flying;
439 : bool acceleration;
440 : double acc;
441 : double acc_restart_mean_;
442 : bool calc_max_bias_;
443 : double max_bias_;
444 : bool calc_transition_bias_;
445 : double transition_bias_;
446 : vector<vector<double> > transitionwells_;
447 : vector<std::unique_ptr<IFile>> ifiles;
448 : vector<string> ifilesnames;
449 : double uppI_;
450 : double lowI_;
451 : bool doInt_;
452 : bool isFirstStep;
453 : bool calc_rct_;
454 : double reweight_factor_;
455 : unsigned rct_ustride_;
456 : double work_;
457 : long int last_step_warn_grid;
458 :
459 : static void registerTemperingKeywords(const std::string &name_stem, const std::string &name, Keywords &keys);
460 : void readTemperingSpecs(TemperingSpecs &t_specs);
461 : void logTemperingSpecs(const TemperingSpecs &t_specs);
462 : void readGaussians(IFile*);
463 : void writeGaussian(const Gaussian&,OFile&);
464 : void addGaussian(const Gaussian&);
465 : double getHeight(const vector<double>&);
466 : void temperHeight(double &height, const TemperingSpecs &t_specs, const double tempering_bias);
467 : double getBiasAndDerivatives(const vector<double>&,double* der=NULL);
468 : double evaluateGaussian(const vector<double>&, const Gaussian&,double* der=NULL);
469 : double getGaussianNormalization( const Gaussian& );
470 : vector<unsigned> getGaussianSupport(const Gaussian&);
471 : bool scanOneHill(IFile *ifile, vector<Value> &v, vector<double> ¢er, vector<double> &sigma, double &height, bool &multivariate);
472 : void computeReweightingFactor();
473 : double getTransitionBarrierBias();
474 : void updateFrequencyAdaptiveStride();
475 : string fmt;
476 :
477 : public:
478 : explicit MetaD(const ActionOptions&);
479 : void calculate();
480 : void update();
481 : static void registerKeywords(Keywords& keys);
482 : bool checkNeedsGradients()const;
483 : };
484 :
485 7639 : PLUMED_REGISTER_ACTION(MetaD,"METAD")
486 :
487 146 : void MetaD::registerKeywords(Keywords& keys) {
488 146 : Bias::registerKeywords(keys);
489 584 : keys.addOutputComponent("rbias","CALC_RCT","the instantaneous value of the bias normalized using the \\f$c(t)\\f$ reweighting factor [rbias=bias-rct]."
490 : "This component can be used to obtain a reweighted histogram.");
491 584 : keys.addOutputComponent("rct","CALC_RCT","the reweighting factor \\f$c(t)\\f$.");
492 584 : keys.addOutputComponent("work","default","accumulator for work");
493 584 : keys.addOutputComponent("acc","ACCELERATION","the metadynamics acceleration factor");
494 584 : keys.addOutputComponent("maxbias", "CALC_MAX_BIAS", "the maximum of the metadynamics V(s, t)");
495 584 : keys.addOutputComponent("transbias", "CALC_TRANSITION_BIAS", "the metadynamics transition bias V*(t)");
496 584 : keys.addOutputComponent("pace","FREQUENCY_ADAPTIVE","the hill addition frequency when employing frequency adaptive metadynamics");
497 292 : keys.use("ARG");
498 584 : keys.add("compulsory","SIGMA","the widths of the Gaussian hills");
499 584 : keys.add("compulsory","PACE","the frequency for hill addition");
500 730 : keys.add("compulsory","FILE","HILLS","a file in which the list of added hills is stored");
501 584 : keys.add("optional","HEIGHT","the heights of the Gaussian hills. Compulsory unless TAU and either BIASFACTOR or DAMPFACTOR are given");
502 584 : keys.add("optional","FMT","specify format for HILLS files (useful for decrease the number of digits in regtests)");
503 584 : keys.add("optional","BIASFACTOR","use well tempered metadynamics and use this bias factor. Please note you must also specify temp");
504 584 : keys.add("optional","RECT","list of bias factors for all the replicas");
505 584 : keys.add("optional","DAMPFACTOR","damp hills with exp(-max(V)/(\\f$k_B\\f$T*DAMPFACTOR)");
506 438 : for (size_t i = 0; i < n_tempering_options_; i++) {
507 146 : registerTemperingKeywords(tempering_names_[i][0], tempering_names_[i][1], keys);
508 : }
509 584 : keys.add("optional","TARGET","target to a predefined distribution");
510 584 : keys.add("optional","TEMP","the system temperature - this is only needed if you are doing well-tempered metadynamics");
511 584 : keys.add("optional","TAU","in well tempered metadynamics, sets height to (\\f$k_B \\Delta T\\f$*pace*timestep)/tau");
512 584 : keys.add("optional","GRID_MIN","the lower bounds for the grid");
513 584 : keys.add("optional","GRID_MAX","the upper bounds for the grid");
514 584 : keys.add("optional","GRID_BIN","the number of bins for the grid");
515 584 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
516 438 : keys.addFlag("CALC_RCT",false,"calculate the \\f$c(t)\\f$ reweighting factor and use that to obtain the normalized bias [rbias=bias-rct]."
517 : "This method is not compatible with metadynamics not on a grid.");
518 584 : keys.add("optional","RCT_USTRIDE","the update stride for calculating the \\f$c(t)\\f$ reweighting factor."
519 : "The default 1, so \\f$c(t)\\f$ is updated every time the bias is updated.");
520 438 : keys.addFlag("GRID_SPARSE",false,"use a sparse grid to store hills");
521 438 : keys.addFlag("GRID_NOSPLINE",false,"don't use spline interpolation with grids");
522 584 : keys.add("optional","GRID_WSTRIDE","write the grid to a file every N steps");
523 584 : keys.add("optional","GRID_WFILE","the file on which to write the grid");
524 584 : keys.add("optional","GRID_RFILE","a grid file from which the bias should be read at the initial step of the simulation");
525 438 : keys.addFlag("STORE_GRIDS",false,"store all the grid files the calculation generates. They will be deleted if this keyword is not present");
526 584 : keys.add("optional","ADAPTIVE","use a geometric (=GEOM) or diffusion (=DIFF) based hills width scheme. Sigma is one number that has distance units or time step dimensions");
527 584 : keys.add("optional","WALKERS_ID", "walker id");
528 584 : keys.add("optional","WALKERS_N", "number of walkers");
529 584 : keys.add("optional","WALKERS_DIR", "shared directory with the hills files from all the walkers");
530 584 : keys.add("optional","WALKERS_RSTRIDE","stride for reading hills files");
531 584 : keys.add("optional","INTERVAL","one dimensional lower and upper limits, outside the limits the system will not feel the biasing force.");
532 584 : keys.add("optional","SIGMA_MAX","the upper bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
533 584 : keys.add("optional","SIGMA_MIN","the lower bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
534 438 : keys.addFlag("WALKERS_MPI",false,"Switch on MPI version of multiple walkers - not compatible with WALKERS_* options other than WALKERS_DIR");
535 438 : keys.addFlag("FLYING_GAUSSIAN",false,"Switch on flying Gaussian method, must be used with WALKERS_MPI");
536 438 : keys.addFlag("ACCELERATION",false,"Set to TRUE if you want to compute the metadynamics acceleration factor.");
537 584 : keys.add("optional","ACCELERATION_RFILE","a data file from which the acceleration should be read at the initial step of the simulation");
538 438 : keys.addFlag("CALC_MAX_BIAS", false, "Set to TRUE if you want to compute the maximum of the metadynamics V(s, t)");
539 438 : keys.addFlag("CALC_TRANSITION_BIAS", false, "Set to TRUE if you want to compute a metadynamics transition bias V*(t)");
540 584 : keys.add("numbered", "TRANSITIONWELL", "This keyword appears multiple times as TRANSITIONWELL followed by an integer. Each specifies the coordinates for one well as in transition-tempered metadynamics. At least one must be provided.");
541 438 : keys.addFlag("FREQUENCY_ADAPTIVE",false,"Set to TRUE if you want to enable frequency adaptive metadynamics such that the frequency for hill addition to change dynamically based on the acceleration factor.");
542 584 : keys.add("optional","FA_UPDATE_FREQUENCY","the frequency for updating the hill addition pace in frequency adaptive metadynamics, by default this is equal to the value given in PACE");
543 584 : keys.add("optional","FA_MAX_PACE","the maximum hill addition frequency allowed in frequency adaptive metadynamics. By default there is no maximum value.");
544 584 : keys.add("optional","FA_MIN_ACCELERATION","only update the hill addition pace in frequency adaptive metadynamics after reaching the minimum acceleration factor given here. By default it is 1.0.");
545 292 : keys.use("RESTART");
546 292 : keys.use("UPDATE_FROM");
547 292 : keys.use("UPDATE_UNTIL");
548 146 : }
549 :
550 5517 : const std::string MetaD::tempering_names_[1][2] = {{"TT", "transition tempered"}};
551 :
552 146 : void MetaD::registerTemperingKeywords(const std::string &name_stem, const std::string &name, Keywords &keys) {
553 730 : keys.add("optional", name_stem + "BIASFACTOR", "use " + name + " metadynamics with this bias factor. Please note you must also specify temp");
554 1022 : keys.add("optional", name_stem + "BIASTHRESHOLD", "use " + name + " metadynamics with this bias threshold. Please note you must also specify " + name_stem + "BIASFACTOR");
555 1022 : keys.add("optional", name_stem + "ALPHA", "use " + name + " metadynamics with this hill size decay exponent parameter. Please note you must also specify " + name_stem + "BIASFACTOR");
556 146 : }
557 :
558 145 : MetaD::MetaD(const ActionOptions& ao):
559 : PLUMED_BIAS_INIT(ao),
560 : // Grid stuff initialization
561 : wgridstride_(0), grid_(false),
562 : // Metadynamics basic parameters
563 : height0_(std::numeric_limits<double>::max()), biasf_(-1.0), dampfactor_(0.0),
564 : tt_specs_(false, "TT", "Transition Tempered", -1.0, 0.0, 1.0),
565 : kbt_(0.0),
566 : stride_(0), welltemp_(false),
567 : // frequency adaptive
568 : current_stride(0),
569 : freq_adaptive_(false),
570 : fa_update_frequency_(0),
571 : fa_max_stride_(0),
572 : fa_min_acceleration_(1.0),
573 : // Other stuff
574 : adaptive_(FlexibleBin::none),
575 : // Multiple walkers initialization
576 : mw_n_(1), mw_dir_(""), mw_id_(0), mw_rstride_(1),
577 : walkers_mpi(false), mpi_nw_(0), mpi_mw_(0),
578 : // Flying Gaussian
579 : flying(false),
580 : acceleration(false), acc(0.0), acc_restart_mean_(0.0),
581 : calc_max_bias_(false), max_bias_(0.0),
582 : calc_transition_bias_(false), transition_bias_(0.0),
583 : // Interval initialization
584 : uppI_(-1), lowI_(-1), doInt_(false),
585 : isFirstStep(true),
586 : calc_rct_(false),
587 : reweight_factor_(0.0),
588 : rct_ustride_(1),
589 : work_(0),
590 2047 : last_step_warn_grid(0)
591 : {
592 : // parse the flexible hills
593 : string adaptiveoption;
594 : adaptiveoption="NONE";
595 288 : parse("ADAPTIVE",adaptiveoption);
596 144 : if(adaptiveoption=="GEOM") {
597 22 : log.printf(" Uses Geometry-based hills width: sigma must be in distance units and only one sigma is needed\n");
598 22 : adaptive_=FlexibleBin::geometry;
599 122 : } else if(adaptiveoption=="DIFF") {
600 3 : log.printf(" Uses Diffusion-based hills width: sigma must be in time steps and only one sigma is needed\n");
601 3 : adaptive_=FlexibleBin::diffusion;
602 119 : } else if(adaptiveoption=="NONE") {
603 118 : adaptive_=FlexibleBin::none;
604 : } else {
605 2 : error("I do not know this type of adaptive scheme");
606 : }
607 :
608 286 : parse("FMT",fmt);
609 :
610 : // parse the sigma
611 286 : parseVector("SIGMA",sigma0_);
612 143 : if(adaptive_==FlexibleBin::none) {
613 : // if you use normal sigma you need one sigma per argument
614 118 : if( sigma0_.size()!=getNumberOfArguments() ) error("number of arguments does not match number of SIGMA parameters");
615 : } else {
616 : // if you use flexible hills you need one sigma
617 25 : if(sigma0_.size()!=1) {
618 2 : error("If you choose ADAPTIVE you need only one sigma according to your choice of type (GEOM/DIFF)");
619 : }
620 : // if adaptive then the number must be an integer
621 24 : if(adaptive_==FlexibleBin::diffusion) {
622 3 : if(int(sigma0_[0])-sigma0_[0]>1.e-9 || int(sigma0_[0])-sigma0_[0] <-1.e-9 || int(sigma0_[0])<1 ) {
623 0 : error("In case of adaptive hills with diffusion, the sigma must be an integer which is the number of time steps\n");
624 : }
625 : }
626 : // here evtl parse the sigma min and max values
627 48 : parseVector("SIGMA_MIN",sigma0min_);
628 25 : if(sigma0min_.size()>0 && sigma0min_.size()!=getNumberOfArguments()) {
629 2 : error("the number of SIGMA_MIN values be the same of the number of the arguments");
630 23 : } else if(sigma0min_.size()==0) {
631 23 : sigma0min_.resize(getNumberOfArguments());
632 111 : for(unsigned i=0; i<getNumberOfArguments(); i++) {sigma0min_[i]=-1.;}
633 : }
634 :
635 46 : parseVector("SIGMA_MAX",sigma0max_);
636 24 : if(sigma0max_.size()>0 && sigma0max_.size()!=getNumberOfArguments()) {
637 2 : error("the number of SIGMA_MAX values be the same of the number of the arguments");
638 22 : } else if(sigma0max_.size()==0) {
639 22 : sigma0max_.resize(getNumberOfArguments());
640 106 : for(unsigned i=0; i<getNumberOfArguments(); i++) {sigma0max_[i]=-1.;}
641 : }
642 :
643 22 : flexbin.reset(new FlexibleBin(adaptive_,this,sigma0_[0],sigma0min_,sigma0max_));
644 : }
645 : // note: HEIGHT is not compulsory, since one could use the TAU keyword, see below
646 280 : parse("HEIGHT",height0_);
647 280 : parse("PACE",stride_);
648 139 : if(stride_<=0 ) error("frequency for hill addition is nonsensical");
649 139 : current_stride = stride_;
650 145 : string hillsfname="HILLS";
651 278 : parse("FILE",hillsfname);
652 :
653 : // Manually set to calculate special bias quantities
654 : // throughout the course of simulation. (These are chosen due to
655 : // relevance for tempering and event-driven logic as well.)
656 278 : parseFlag("CALC_MAX_BIAS", calc_max_bias_);
657 278 : parseFlag("CALC_TRANSITION_BIAS", calc_transition_bias_);
658 :
659 : std::vector<double> rect_biasf_;
660 278 : parseVector("RECT",rect_biasf_);
661 139 : if(rect_biasf_.size()>0) {
662 18 : int r=0;
663 18 : if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
664 18 : comm.Bcast(r,0);
665 36 : biasf_=rect_biasf_[r];
666 18 : log<<" You are using RECT\n";
667 : } else {
668 242 : parse("BIASFACTOR",biasf_);
669 : }
670 139 : if( biasf_<1.0 && biasf_!=-1.0) error("well tempered bias factor is nonsensical");
671 278 : parse("DAMPFACTOR",dampfactor_);
672 139 : double temp=0.0;
673 278 : parse("TEMP",temp);
674 188 : if(temp>0.0) kbt_=plumed.getAtoms().getKBoltzmann()*temp;
675 180 : else kbt_=plumed.getAtoms().getKbT();
676 139 : if(biasf_>=1.0) {
677 32 : if(kbt_==0.0) error("Unless the MD engine passes the temperature to plumed, with well-tempered metad you must specify it using TEMP");
678 32 : welltemp_=true;
679 : }
680 139 : if(dampfactor_>0.0) {
681 2 : if(kbt_==0.0) error("Unless the MD engine passes the temperature to plumed, with damped metad you must specify it using TEMP");
682 : }
683 :
684 : // Set transition tempering parameters.
685 : // Transition wells are read later via calc_transition_bias_.
686 139 : readTemperingSpecs(tt_specs_);
687 139 : if (tt_specs_.is_active) calc_transition_bias_ = true;
688 :
689 : // If any previous option specified to calculate a transition bias,
690 : // now read the transition wells for that quantity.
691 139 : if (calc_transition_bias_) {
692 14 : vector<double> tempcoords(getNumberOfArguments());
693 26 : for (unsigned i = 0; ; i++) {
694 104 : if (!parseNumberedVector("TRANSITIONWELL", i, tempcoords) ) break;
695 26 : if (tempcoords.size() != getNumberOfArguments()) {
696 0 : error("incorrect number of coordinates for transition tempering well");
697 : }
698 26 : transitionwells_.push_back(tempcoords);
699 : }
700 : }
701 :
702 278 : parse("TARGET",targetfilename_);
703 139 : if(targetfilename_.length()>0 && kbt_==0.0) error("with TARGET temperature must be specified");
704 139 : double tau=0.0;
705 278 : parse("TAU",tau);
706 139 : if(tau==0.0) {
707 117 : if(height0_==std::numeric_limits<double>::max()) error("At least one between HEIGHT and TAU should be specified");
708 : // if tau is not set, we compute it here from the other input parameters
709 117 : if(welltemp_) tau=(kbt_*(biasf_-1.0))/height0_*getTimeStep()*stride_;
710 104 : else if(dampfactor_>0.0) tau=(kbt_*dampfactor_)/height0_*getTimeStep()*stride_;
711 : } else {
712 22 : if(height0_!=std::numeric_limits<double>::max()) error("At most one between HEIGHT and TAU should be specified");
713 22 : if(welltemp_) {
714 19 : if(biasf_!=1.0) height0_=(kbt_*(biasf_-1.0))/tau*getTimeStep()*stride_;
715 4 : else height0_=kbt_/tau*getTimeStep()*stride_; // special case for gamma=1
716 : }
717 3 : else if(dampfactor_>0.0) height0_=(kbt_*dampfactor_)/tau*getTimeStep()*stride_;
718 2 : else error("TAU only makes sense in well-tempered or damped metadynamics");
719 : }
720 :
721 : // Grid Stuff
722 276 : vector<std::string> gmin(getNumberOfArguments());
723 276 : parseVector("GRID_MIN",gmin);
724 138 : if(gmin.size()!=getNumberOfArguments() && gmin.size()!=0) error("not enough values for GRID_MIN");
725 276 : vector<std::string> gmax(getNumberOfArguments());
726 276 : parseVector("GRID_MAX",gmax);
727 138 : if(gmax.size()!=getNumberOfArguments() && gmax.size()!=0) error("not enough values for GRID_MAX");
728 138 : vector<unsigned> gbin(getNumberOfArguments());
729 : vector<double> gspacing;
730 276 : parseVector("GRID_BIN",gbin);
731 138 : if(gbin.size()!=getNumberOfArguments() && gbin.size()!=0) error("not enough values for GRID_BIN");
732 276 : parseVector("GRID_SPACING",gspacing);
733 138 : if(gspacing.size()!=getNumberOfArguments() && gspacing.size()!=0) error("not enough values for GRID_SPACING");
734 138 : if(gmin.size()!=gmax.size()) error("GRID_MAX and GRID_MIN should be either present or absent");
735 140 : if(gspacing.size()!=0 && gmin.size()==0) error("If GRID_SPACING is present also GRID_MIN should be present");
736 188 : if(gbin.size()!=0 && gmin.size()==0) error("If GRID_SPACING is present also GRID_MIN should be present");
737 138 : if(gmin.size()!=0) {
738 54 : if(gbin.size()==0 && gspacing.size()==0) {
739 1 : if(adaptive_==FlexibleBin::none) {
740 1 : log<<" Binsize not specified, 1/5 of sigma will be be used\n";
741 1 : plumed_assert(sigma0_.size()==getNumberOfArguments());
742 1 : gspacing.resize(getNumberOfArguments());
743 10 : for(unsigned i=0; i<gspacing.size(); i++) gspacing[i]=0.2*sigma0_[i];
744 : } else {
745 : // with adaptive hills and grid a sigma min must be specified
746 0 : for(unsigned i=0; i<sigma0min_.size(); i++) if(sigma0min_[i]<=0) error("When using Adaptive Gaussians on a grid SIGMA_MIN must be specified");
747 0 : log<<" Binsize not specified, 1/5 of sigma_min will be be used\n";
748 0 : gspacing.resize(getNumberOfArguments());
749 0 : for(unsigned i=0; i<gspacing.size(); i++) gspacing[i]=0.2*sigma0min_[i];
750 : }
751 53 : } else if(gspacing.size()!=0 && gbin.size()==0) {
752 1 : log<<" The number of bins will be estimated from GRID_SPACING\n";
753 51 : } else if(gspacing.size()!=0 && gbin.size()!=0) {
754 1 : log<<" You specified both GRID_BIN and GRID_SPACING\n";
755 1 : log<<" The more conservative (highest) number of bins will be used for each variable\n";
756 : }
757 56 : if(gbin.size()==0) gbin.assign(getNumberOfArguments(),1);
758 67 : if(gspacing.size()!=0) for(unsigned i=0; i<getNumberOfArguments(); i++) {
759 : double a,b;
760 12 : Tools::convert(gmin[i],a);
761 6 : Tools::convert(gmax[i],b);
762 12 : unsigned n=((b-a)/gspacing[i])+1;
763 6 : if(gbin[i]<n) gbin[i]=n;
764 : }
765 : }
766 138 : bool sparsegrid=false;
767 276 : parseFlag("GRID_SPARSE",sparsegrid);
768 138 : bool nospline=false;
769 276 : parseFlag("GRID_NOSPLINE",nospline);
770 138 : bool spline=!nospline;
771 138 : if(gbin.size()>0) {grid_=true;}
772 276 : parse("GRID_WSTRIDE",wgridstride_);
773 : string gridfilename_;
774 276 : parse("GRID_WFILE",gridfilename_);
775 276 : parseFlag("STORE_GRIDS",storeOldGrids_);
776 190 : if(grid_ && gridfilename_.length()>0) {
777 16 : if(wgridstride_==0 ) error("frequency with which to output grid not specified use GRID_WSTRIDE");
778 : }
779 :
780 138 : if(grid_ && wgridstride_>0) {
781 16 : if(gridfilename_.length()==0) error("grid filename not specified use GRID_WFILE");
782 : }
783 : string gridreadfilename_;
784 276 : parse("GRID_RFILE",gridreadfilename_);
785 :
786 224 : if(!grid_&&gridfilename_.length()> 0) error("To write a grid you need first to define it!");
787 224 : if(!grid_&&gridreadfilename_.length()>0) error("To read a grid you need first to define it!");
788 :
789 : // Reweighting factor rct
790 276 : parseFlag("CALC_RCT",calc_rct_);
791 138 : if (calc_rct_)
792 5 : plumed_massert(grid_,"CALC_RCT is supported only if bias is on a grid");
793 276 : parse("RCT_USTRIDE",rct_ustride_);
794 :
795 138 : if(dampfactor_>0.0) {
796 2 : if(!grid_) error("With DAMPFACTOR you should use grids");
797 : }
798 :
799 : // Multiple walkers
800 276 : parse("WALKERS_N",mw_n_);
801 276 : parse("WALKERS_ID",mw_id_);
802 138 : if(mw_n_<=mw_id_) error("walker ID should be a numerical value less than the total number of walkers");
803 276 : parse("WALKERS_DIR",mw_dir_);
804 276 : parse("WALKERS_RSTRIDE",mw_rstride_);
805 :
806 : // MPI version
807 276 : parseFlag("WALKERS_MPI",walkers_mpi);
808 :
809 : // Flying Gaussian
810 276 : parseFlag("FLYING_GAUSSIAN", flying);
811 :
812 : // Inteval keyword
813 138 : vector<double> tmpI(2);
814 276 : parseVector("INTERVAL",tmpI);
815 138 : if(tmpI.size()!=2&&tmpI.size()!=0) error("both a lower and an upper limits must be provided with INTERVAL");
816 138 : else if(tmpI.size()==2) {
817 2 : lowI_=tmpI.at(0);
818 2 : uppI_=tmpI.at(1);
819 2 : if(getNumberOfArguments()!=1) error("INTERVAL limits correction works only for monodimensional metadynamics!");
820 2 : if(uppI_<lowI_) error("The Upper limit must be greater than the Lower limit!");
821 2 : if(getPntrToArgument(0)->isPeriodic()) error("INTERVAL cannot be used with periodic variables!");
822 2 : doInt_=true;
823 : }
824 :
825 138 : acceleration=false;
826 276 : parseFlag("ACCELERATION",acceleration);
827 : // Check for a restart acceleration if acceleration is active.
828 : string acc_rfilename;
829 138 : if (acceleration) {
830 8 : parse("ACCELERATION_RFILE", acc_rfilename);
831 : }
832 :
833 138 : freq_adaptive_=false;
834 276 : parseFlag("FREQUENCY_ADAPTIVE",freq_adaptive_);
835 : //
836 138 : fa_update_frequency_=0;
837 276 : parse("FA_UPDATE_FREQUENCY",fa_update_frequency_);
838 138 : if(fa_update_frequency_!=0 && !freq_adaptive_) {
839 0 : plumed_merror("It doesn't make sense to use the FA_MAX_PACE keyword if frequency adaptive MetaD hasn't been activated by using the FREQUENCY_ADAPTIVE flag");
840 : }
841 138 : if(fa_update_frequency_==0 && freq_adaptive_) {
842 0 : fa_update_frequency_=stride_;
843 : }
844 : //
845 138 : fa_max_stride_=0;
846 276 : parse("FA_MAX_PACE",fa_max_stride_);
847 138 : if(fa_max_stride_!=0 && !freq_adaptive_) {
848 0 : plumed_merror("It doesn't make sense to use the FA_MAX_PACE keyword if frequency adaptive MetaD hasn't been activated by using the FREQUENCY_ADAPTIVE flag");
849 : }
850 : //
851 138 : fa_min_acceleration_=1.0;
852 276 : parse("FA_MIN_ACCELERATION",fa_min_acceleration_);
853 138 : if(fa_min_acceleration_!=1.0 && !freq_adaptive_) {
854 0 : plumed_merror("It doesn't make sense to use the FA_MIN_ACCELERATION keyword if frequency adaptive MetaD hasn't been activated by using the FREQUENCY_ADAPTIVE flag");
855 : }
856 :
857 138 : checkRead();
858 :
859 138 : log.printf(" Gaussian width ");
860 138 : if (adaptive_==FlexibleBin::diffusion)log.printf(" (Note: The units of sigma are in timesteps) ");
861 138 : if (adaptive_==FlexibleBin::geometry)log.printf(" (Note: The units of sigma are in dist units) ");
862 975 : for(unsigned i=0; i<sigma0_.size(); ++i) log.printf(" %f",sigma0_[i]);
863 138 : log.printf(" Gaussian height %f\n",height0_);
864 138 : log.printf(" Gaussian deposition pace %d\n",stride_);
865 276 : log.printf(" Gaussian file %s\n",hillsfname.c_str());
866 138 : if(welltemp_) {
867 32 : log.printf(" Well-Tempered Bias Factor %f\n",biasf_);
868 32 : log.printf(" Hills relaxation time (tau) %f\n",tau);
869 32 : log.printf(" KbT %f\n",kbt_);
870 : }
871 : // Transition tempered metadynamics options
872 138 : if (tt_specs_.is_active) {
873 3 : logTemperingSpecs(tt_specs_);
874 : // Check that the appropriate transition bias quantity is calculated.
875 : // (Should never trip, given that the flag is automatically set.)
876 3 : if (!calc_transition_bias_) {
877 0 : error(" transition tempering requires calculation of a transition bias");
878 : }
879 : }
880 :
881 : // Overall tempering sanity check (this gets tricky when multiple are active).
882 : // When multiple temperings are active, it's fine to have one tempering attempt
883 : // to increase hill size with increasing bias, so long as the others can shrink
884 : // the hills faster than it increases their size in the long-time limit.
885 : // This set of checks ensures that the hill sizes eventually decay to zero as c(t)
886 : // diverges to infinity.
887 : // The alpha parameter allows hills to decay as 1/t^alpha instead of 1/t,
888 : // a slower decay, so as t -> infinity, only the temperings with the largest
889 : // alphas govern the final asymptotic decay. (Alpha helps prevent false convergence.)
890 138 : if (welltemp_ || dampfactor_ > 0.0 || tt_specs_.is_active) {
891 : // Determine the number of active temperings.
892 : int n_active = 0;
893 37 : if (welltemp_) n_active++;
894 37 : if (dampfactor_ > 0.0) n_active++;
895 37 : if (tt_specs_.is_active) n_active++;
896 : // Find the greatest alpha.
897 37 : double greatest_alpha = 0.0;
898 37 : if (welltemp_) greatest_alpha = max(greatest_alpha, 1.0);
899 39 : if (dampfactor_ > 0.0) greatest_alpha = max(greatest_alpha, 1.0);
900 40 : if (tt_specs_.is_active) greatest_alpha = max(greatest_alpha, tt_specs_.alpha);
901 : // Find the least alpha.
902 37 : double least_alpha = 1.0;
903 : if (welltemp_) least_alpha = min(least_alpha, 1.0);
904 39 : if (dampfactor_ > 0.0) least_alpha = min(least_alpha, 1.0);
905 40 : if (tt_specs_.is_active) least_alpha = min(least_alpha, tt_specs_.alpha);
906 : // Find the inverse harmonic average of the delta T parameters for all
907 : // of the temperings with the greatest alpha values.
908 : double total_governing_deltaT_inv = 0.0;
909 37 : if (welltemp_ && 1.0 == greatest_alpha && biasf_ != 1.0) total_governing_deltaT_inv += 1.0 / (biasf_ - 1.0);
910 37 : if (dampfactor_ > 0.0 && 1.0 == greatest_alpha) total_governing_deltaT_inv += 1.0 / (dampfactor_);
911 37 : if (tt_specs_.is_active && tt_specs_.alpha == greatest_alpha) total_governing_deltaT_inv += 1.0 / (tt_specs_.biasf - 1.0);
912 : // Give a newbie-friendly error message for people using one tempering if
913 : // only one is active.
914 37 : if (n_active == 1 && total_governing_deltaT_inv < 0.0) {
915 0 : error("for stable tempering, the bias factor must be greater than one");
916 : // Give a slightly more complex error message to users stacking multiple
917 : // tempering options at a time, but all with uniform alpha values.
918 37 : } else if (total_governing_deltaT_inv < 0.0 && greatest_alpha == least_alpha) {
919 0 : error("for stable tempering, the sum of the inverse Delta T parameters must be greater than zero!");
920 : // Give the most technical error message to users stacking multiple tempering
921 : // options with different alpha parameters.
922 37 : } else if (total_governing_deltaT_inv < 0.0 && greatest_alpha != least_alpha) {
923 0 : error("for stable tempering, the sum of the inverse Delta T parameters for the greatest asymptotic hill decay exponents must be greater than zero!");
924 : }
925 : }
926 :
927 138 : if(doInt_) log.printf(" Upper and Lower limits boundaries for the bias are activated at %f - %f\n", lowI_, uppI_);
928 138 : if(grid_) {
929 52 : log.printf(" Grid min");
930 371 : for(unsigned i=0; i<gmin.size(); ++i) log.printf(" %s",gmin[i].c_str() );
931 52 : log.printf("\n");
932 52 : log.printf(" Grid max");
933 371 : for(unsigned i=0; i<gmax.size(); ++i) log.printf(" %s",gmax[i].c_str() );
934 52 : log.printf("\n");
935 52 : log.printf(" Grid bin");
936 371 : for(unsigned i=0; i<gbin.size(); ++i) log.printf(" %u",gbin[i]);
937 52 : log.printf("\n");
938 52 : if(spline) {log.printf(" Grid uses spline interpolation\n");}
939 52 : if(sparsegrid) {log.printf(" Grid uses sparse grid\n");}
940 68 : if(wgridstride_>0) {log.printf(" Grid is written on file %s with stride %d\n",gridfilename_.c_str(),wgridstride_);}
941 : }
942 :
943 138 : if(mw_n_>1) {
944 6 : if(walkers_mpi) error("MPI version of multiple walkers is not compatible with filesystem version of multiple walkers");
945 6 : log.printf(" %d multiple walkers active\n",mw_n_);
946 6 : log.printf(" walker id %d\n",mw_id_);
947 6 : log.printf(" reading stride %d\n",mw_rstride_);
948 9 : if(mw_dir_!="")log.printf(" directory with hills files %s\n",mw_dir_.c_str());
949 : } else {
950 132 : if(walkers_mpi) {
951 36 : log.printf(" Multiple walkers active using MPI communnication\n");
952 36 : if(mw_dir_!="")log.printf(" directory with hills files %s\n",mw_dir_.c_str());
953 36 : if(comm.Get_rank()==0) {
954 : // Only root of group can communicate with other walkers
955 21 : mpi_nw_=multi_sim_comm.Get_size();
956 21 : mpi_mw_=multi_sim_comm.Get_rank();
957 : }
958 : // Communicate to the other members of the same group
959 : // info abount number of walkers and walker index
960 36 : comm.Bcast(mpi_nw_,0);
961 36 : comm.Bcast(mpi_mw_,0);
962 : }
963 : }
964 :
965 138 : if(flying) {
966 6 : if(!walkers_mpi) error("Flying Gaussian method must be used with MPI version of multiple walkers");
967 6 : log.printf(" Flying Gaussian method with %d walkers active\n",mpi_nw_);
968 : }
969 :
970 138 : if(calc_rct_) {
971 15 : addComponent("rbias"); componentIsNotPeriodic("rbias");
972 15 : addComponent("rct"); componentIsNotPeriodic("rct");
973 5 : log.printf(" The c(t) reweighting factor will be calculated every %u hills\n",rct_ustride_);
974 10 : getPntrToComponent("rct")->set(reweight_factor_);
975 : }
976 414 : addComponent("work"); componentIsNotPeriodic("work");
977 :
978 138 : if(acceleration) {
979 4 : if (kbt_ == 0.0) {
980 0 : error("The calculation of the acceleration works only if simulation temperature has been defined");
981 : }
982 4 : log.printf(" calculation on the fly of the acceleration factor\n");
983 12 : addComponent("acc"); componentIsNotPeriodic("acc");
984 : // Set the initial value of the the acceleration.
985 : // If this is not a restart, set to 1.0.
986 4 : if (acc_rfilename.length() == 0) {
987 4 : getPntrToComponent("acc")->set(1.0);
988 2 : if(getRestart()) {
989 1 : log.printf(" WARNING: calculating the acceleration factor in a restarted run without reading in the previous value will most likely lead to incorrect results. You should use the ACCELERATION_RFILE keyword.\n");
990 : }
991 : // Otherwise, read and set the restart value.
992 : } else {
993 : // Restart of acceleration does not make sense if the restart timestep is zero.
994 : //if (getStep() == 0) {
995 : // error("Restarting calculation of acceleration factors works only if simulation timestep is restarted correctly");
996 : //}
997 : // Open the ACCELERATION_RFILE.
998 4 : IFile acc_rfile;
999 2 : acc_rfile.link(*this);
1000 2 : if(acc_rfile.FileExist(acc_rfilename)) {
1001 2 : acc_rfile.open(acc_rfilename);
1002 : } else {
1003 0 : error("The ACCELERATION_RFILE file you want to read: " + acc_rfilename + ", cannot be found!");
1004 : }
1005 : // Read the file to find the restart acceleration.
1006 : double acc_rmean;
1007 : double acc_rtime;
1008 2 : std::string acclabel = getLabel() + ".acc";
1009 2 : acc_rfile.allowIgnoredFields();
1010 248 : while(acc_rfile.scanField("time", acc_rtime)) {
1011 122 : acc_rfile.scanField(acclabel, acc_rmean);
1012 122 : acc_rfile.scanField();
1013 : }
1014 2 : acc_restart_mean_ = acc_rmean;
1015 : // Set component based on the read values.
1016 4 : getPntrToComponent("acc")->set(acc_rmean);
1017 4 : log.printf(" initial acceleration factor read from file %s: value of %f at time %f\n",acc_rfilename.c_str(),acc_rmean,acc_rtime);
1018 : }
1019 : }
1020 138 : if (calc_max_bias_) {
1021 0 : if (!grid_) error("Calculating the maximum bias on the fly works only with a grid");
1022 0 : log.printf(" calculation on the fly of the maximum bias max(V(s,t)) \n");
1023 0 : addComponent("maxbias");
1024 0 : componentIsNotPeriodic("maxbias");
1025 : }
1026 138 : if (calc_transition_bias_) {
1027 13 : if (!grid_) error("Calculating the transition bias on the fly works only with a grid");
1028 13 : log.printf(" calculation on the fly of the transition bias V*(t)\n");
1029 26 : addComponent("transbias");
1030 26 : componentIsNotPeriodic("transbias");
1031 26 : log<<" Number of transition wells "<<transitionwells_.size()<<"\n";
1032 13 : if (transitionwells_.size() == 0) error("Calculating the transition bias on the fly requires definition of at least one transition well");
1033 : // Check that a grid is in use.
1034 13 : if (!grid_) error(" transition barrier finding requires a grid for the bias");
1035 : // Log the wells and check that they are in the grid.
1036 104 : for (unsigned i = 0; i < transitionwells_.size(); i++) {
1037 : // Log the coordinate.
1038 26 : log.printf(" Transition well %d at coordinate ", i);
1039 102 : for (unsigned j = 0; j < getNumberOfArguments(); j++) log.printf("%f ", transitionwells_[i][j]);
1040 26 : log.printf("\n");
1041 : // Check that the coordinate is in the grid.
1042 102 : for (unsigned j = 0; j < getNumberOfArguments(); j++) {
1043 : double max, min;
1044 76 : Tools::convert(gmin[j], min);
1045 38 : Tools::convert(gmax[j], max);
1046 38 : if (transitionwells_[i][j] < min || transitionwells_[i][j] > max) error(" transition well is not in grid");
1047 : }
1048 : }
1049 : }
1050 :
1051 138 : if(freq_adaptive_) {
1052 2 : if(!acceleration) {
1053 0 : plumed_merror("Frequency adaptive metadynamics only works if the calculation of the acceleration factor is enabled with the ACCELERATION keyword\n");
1054 : }
1055 2 : if(walkers_mpi) {
1056 0 : plumed_merror("Combining frequency adaptive metadynamics with MPI multiple walkers is not allowed");
1057 : }
1058 :
1059 2 : log.printf(" Frequency adaptive metadynamics enabled\n");
1060 3 : if(getRestart() && acc_rfilename.length() == 0) {
1061 0 : log.printf(" WARNING: using the frequency adaptive scheme in a restarted run without reading in the previous value of the acceleration factor will most likely lead to incorrect results. You should use the ACCELERATION_RFILE keyword.\n");
1062 : }
1063 2 : log.printf(" The frequency for hill addition will change dynamically based on the metadynamics acceleration factor\n");
1064 2 : log.printf(" The hill addition frequency will be updated every %d steps\n",fa_update_frequency_);
1065 2 : if(fa_min_acceleration_>1.0) {
1066 2 : log.printf(" The hill addition frequency will only be updated once the metadynamics acceleration factor becomes larger than %.1f \n",fa_min_acceleration_);
1067 : }
1068 2 : if(fa_max_stride_!=0) {
1069 2 : log.printf(" The hill addition frequency will not become larger than %d steps\n",fa_max_stride_);
1070 : }
1071 6 : addComponent("pace"); componentIsNotPeriodic("pace");
1072 2 : updateFrequencyAdaptiveStride();
1073 : }
1074 :
1075 : // for performance
1076 138 : dp_.reset( new double[getNumberOfArguments()] );
1077 :
1078 : // initializing and checking grid
1079 138 : if(grid_) {
1080 : // check for mesh and sigma size
1081 230 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
1082 : double a,b;
1083 178 : Tools::convert(gmin[i],a);
1084 89 : Tools::convert(gmax[i],b);
1085 178 : double mesh=(b-a)/((double)gbin[i]);
1086 89 : if(adaptive_==FlexibleBin::none) {
1087 89 : 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";
1088 : } else {
1089 0 : if(mesh>0.5*sigma0min_[i]||sigma0min_[i]<0.) log<<" WARNING: to use a METAD with a GRID and ADAPTIVE you need to set a Grid Spacing larger than half of the Gaussians \n";
1090 : }
1091 : }
1092 52 : std::string funcl=getLabel() + ".bias";
1093 98 : if(!sparsegrid) {BiasGrid_.reset(new Grid(funcl,getArguments(),gmin,gmax,gbin,spline,true));}
1094 6 : else {BiasGrid_.reset(new SparseGrid(funcl,getArguments(),gmin,gmax,gbin,spline,true));}
1095 104 : std::vector<std::string> actualmin=BiasGrid_->getMin();
1096 104 : std::vector<std::string> actualmax=BiasGrid_->getMax();
1097 230 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
1098 : std::string is;
1099 89 : Tools::convert(i,is);
1100 178 : if(gmin[i]!=actualmin[i]) error("GRID_MIN["+is+"] must be adjusted to "+actualmin[i]+" to fit periodicity");
1101 89 : if(gmax[i]!=actualmax[i]) error("GRID_MAX["+is+"] must be adjusted to "+actualmax[i]+" to fit periodicity");
1102 : }
1103 : }
1104 :
1105 : // restart from external grid
1106 : bool restartedFromGrid=false;
1107 138 : if(gridreadfilename_.length()>0) {
1108 : // read the grid in input, find the keys
1109 36 : IFile gridfile;
1110 18 : gridfile.link(*this);
1111 18 : if(gridfile.FileExist(gridreadfilename_)) {
1112 18 : gridfile.open(gridreadfilename_);
1113 : } else {
1114 0 : error("The GRID file you want to read: " + gridreadfilename_ + ", cannot be found!");
1115 : }
1116 18 : std::string funcl=getLabel() + ".bias";
1117 54 : BiasGrid_=Grid::create(funcl, getArguments(), gridfile, gmin, gmax, gbin, sparsegrid, spline, true);
1118 36 : if(BiasGrid_->getDimension()!=getNumberOfArguments()) error("mismatch between dimensionality of input grid and number of arguments");
1119 72 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1120 81 : if( getPntrToArgument(i)->isPeriodic()!=BiasGrid_->getIsPeriodic()[i] ) error("periodicity mismatch between arguments and input bias");
1121 : double a, b;
1122 27 : Tools::convert(gmin[i],a);
1123 27 : Tools::convert(gmax[i],b);
1124 54 : double mesh=(b-a)/((double)gbin[i]);
1125 27 : 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";
1126 : }
1127 36 : log.printf(" Restarting from %s:",gridreadfilename_.c_str());
1128 18 : if(getRestart()) restartedFromGrid=true;
1129 : }
1130 :
1131 : // initializing and checking grid
1132 190 : if(grid_&&!(gridreadfilename_.length()>0)) {
1133 : // check for adaptive and sigma_min
1134 34 : if(sigma0min_.size()==0&&adaptive_!=FlexibleBin::none) error("When using Adaptive Gaussians on a grid SIGMA_MIN must be specified");
1135 : // check for mesh and sigma size
1136 158 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
1137 : double a,b;
1138 124 : Tools::convert(gmin[i],a);
1139 62 : Tools::convert(gmax[i],b);
1140 124 : double mesh=(b-a)/((double)gbin[i]);
1141 62 : 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";
1142 : }
1143 34 : std::string funcl=getLabel() + ".bias";
1144 62 : if(!sparsegrid) {BiasGrid_.reset(new Grid(funcl,getArguments(),gmin,gmax,gbin,spline,true));}
1145 6 : else {BiasGrid_.reset(new SparseGrid(funcl,getArguments(),gmin,gmax,gbin,spline,true));}
1146 68 : std::vector<std::string> actualmin=BiasGrid_->getMin();
1147 68 : std::vector<std::string> actualmax=BiasGrid_->getMax();
1148 192 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
1149 124 : if(gmin[i]!=actualmin[i]) log<<" WARNING: GRID_MIN["<<i<<"] has been adjusted to "<<actualmin[i]<<" to fit periodicity\n";
1150 124 : if(gmax[i]!=actualmax[i]) log<<" WARNING: GRID_MAX["<<i<<"] has been adjusted to "<<actualmax[i]<<" to fit periodicity\n";
1151 : }
1152 : }
1153 :
1154 : // creating vector of ifile* for hills reading
1155 : // open all files at the beginning and read Gaussians if restarting
1156 438 : for(int i=0; i<mw_n_; ++i) {
1157 : string fname;
1158 150 : if(mw_dir_!="") {
1159 9 : if(mw_n_>1) {
1160 18 : stringstream out; out << i;
1161 54 : fname = mw_dir_+"/"+hillsfname+"."+out.str();
1162 0 : } else if(walkers_mpi) {
1163 0 : fname = mw_dir_+"/"+hillsfname;
1164 : } else {
1165 : fname = hillsfname;
1166 : }
1167 : } else {
1168 141 : if(mw_n_>1) {
1169 18 : stringstream out; out << i;
1170 36 : fname = hillsfname+"."+out.str();
1171 : } else {
1172 : fname = hillsfname;
1173 : }
1174 : }
1175 150 : ifiles.emplace_back(new IFile());
1176 : // this is just a shortcut pointer to the last element:
1177 : IFile *ifile = ifiles.back().get();
1178 150 : ifilesnames.push_back(fname);
1179 150 : ifile->link(*this);
1180 150 : if(ifile->FileExist(fname)) {
1181 35 : ifile->open(fname);
1182 35 : if(getRestart()&&!restartedFromGrid) {
1183 36 : log.printf(" Restarting from %s:",ifilesnames[i].c_str());
1184 18 : readGaussians(ifiles[i].get());
1185 : }
1186 70 : ifiles[i]->reset(false);
1187 : // close only the walker own hills file for later writing
1188 64 : if(i==mw_id_) ifiles[i]->close();
1189 : } else {
1190 : // in case a file does not exist and we are restarting, complain that the file was not found
1191 115 : if(getRestart()) log<<" WARNING: restart file "<<fname<<" not found\n";
1192 : }
1193 : }
1194 :
1195 138 : comm.Barrier();
1196 : // this barrier is needed when using walkers_mpi
1197 : // to be sure that all files have been read before
1198 : // backing them up
1199 : // it should not be used when walkers_mpi is false otherwise
1200 : // it would introduce troubles when using replicas without METAD
1201 : // (e.g. in bias exchange with a neutral replica)
1202 : // see issue #168 on github
1203 138 : if(comm.Get_rank()==0 && walkers_mpi) multi_sim_comm.Barrier();
1204 138 : if(targetfilename_.length()>0) {
1205 4 : IFile gridfile; gridfile.open(targetfilename_);
1206 2 : std::string funcl=getLabel() + ".target";
1207 4 : TargetGrid_=Grid::create(funcl,getArguments(),gridfile,false,false,true);
1208 4 : if(TargetGrid_->getDimension()!=getNumberOfArguments()) error("mismatch between dimensionality of input grid and number of arguments");
1209 6 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1210 6 : if( getPntrToArgument(i)->isPeriodic()!=TargetGrid_->getIsPeriodic()[i] ) error("periodicity mismatch between arguments and input bias");
1211 : }
1212 : }
1213 :
1214 : // Calculate the Tiwary-Parrinello reweighting factor if we are restarting from previous hills
1215 138 : if(getRestart() && calc_rct_) computeReweightingFactor();
1216 : // Calculate all special bias quantities desired if restarting with nonzero bias.
1217 138 : if(getRestart() && calc_max_bias_) {
1218 0 : max_bias_ = BiasGrid_->getMaxValue();
1219 0 : getPntrToComponent("maxbias")->set(max_bias_);
1220 : }
1221 138 : if(getRestart() && calc_transition_bias_) {
1222 13 : transition_bias_ = getTransitionBarrierBias();
1223 26 : getPntrToComponent("transbias")->set(transition_bias_);
1224 : }
1225 :
1226 : // open grid file for writing
1227 138 : if(wgridstride_>0) {
1228 16 : gridfile_.link(*this);
1229 16 : if(walkers_mpi) {
1230 0 : int r=0;
1231 0 : if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
1232 0 : comm.Bcast(r,0);
1233 0 : if(r>0) gridfilename_="/dev/null";
1234 0 : gridfile_.enforceSuffix("");
1235 : }
1236 16 : if(mw_n_>1) gridfile_.enforceSuffix("");
1237 16 : gridfile_.open(gridfilename_);
1238 : }
1239 :
1240 : // open hills file for writing
1241 138 : hillsOfile_.link(*this);
1242 138 : if(walkers_mpi) {
1243 36 : int r=0;
1244 36 : if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
1245 36 : comm.Bcast(r,0);
1246 36 : if(r>0) ifilesnames[mw_id_]="/dev/null";
1247 72 : hillsOfile_.enforceSuffix("");
1248 : }
1249 144 : if(mw_n_>1) hillsOfile_.enforceSuffix("");
1250 276 : hillsOfile_.open(ifilesnames[mw_id_]);
1251 138 : if(fmt.length()>0) hillsOfile_.fmtField(fmt);
1252 276 : hillsOfile_.addConstantField("multivariate");
1253 276 : hillsOfile_.addConstantField("kerneltype");
1254 138 : if(doInt_) {
1255 6 : hillsOfile_.addConstantField("lower_int").printField("lower_int",lowI_);
1256 6 : hillsOfile_.addConstantField("upper_int").printField("upper_int",uppI_);
1257 : }
1258 : hillsOfile_.setHeavyFlush();
1259 : // output periodicities of variables
1260 642 : for(unsigned i=0; i<getNumberOfArguments(); ++i) hillsOfile_.setupPrintValue( getPntrToArgument(i) );
1261 :
1262 : bool concurrent=false;
1263 138 : const ActionSet&actionSet(plumed.getActionSet());
1264 948 : for(const auto & p : actionSet) if(dynamic_cast<MetaD*>(p.get())) { concurrent=true; break; }
1265 138 : if(concurrent) log<<" You are using concurrent metadynamics\n";
1266 138 : if(rect_biasf_.size()>0) {
1267 18 : if(walkers_mpi) {
1268 12 : log<<" You are using RECT in its 'altruistic' implementation\n";
1269 : }{
1270 18 : log<<" You are using RECT\n";
1271 : }
1272 : }
1273 :
1274 414 : log<<" Bibliography "<<plumed.cite("Laio and Parrinello, PNAS 99, 12562 (2002)");
1275 234 : if(welltemp_) log<<plumed.cite(
1276 32 : "Barducci, Bussi, and Parrinello, Phys. Rev. Lett. 100, 020603 (2008)");
1277 138 : if(tt_specs_.is_active) {
1278 9 : log << plumed.cite("Dama, Rotskoff, Parrinello, and Voth, J. Chem. Theory Comput. 10, 3626 (2014)");
1279 9 : log << plumed.cite("Dama, Parrinello, and Voth, Phys. Rev. Lett. 112, 240602 (2014)");
1280 : }
1281 264 : if(mw_n_>1||walkers_mpi) log<<plumed.cite(
1282 42 : "Raiteri, Laio, Gervasio, Micheletti, and Parrinello, J. Phys. Chem. B 110, 3533 (2006)");
1283 201 : if(adaptive_!=FlexibleBin::none) log<<plumed.cite(
1284 21 : "Branduardi, Bussi, and Parrinello, J. Chem. Theory Comput. 8, 2247 (2012)");
1285 144 : if(doInt_) log<<plumed.cite(
1286 2 : "Baftizadeh, Cossio, Pietrucci, and Laio, Curr. Phys. Chem. 2, 79 (2012)");
1287 150 : if(acceleration) log<<plumed.cite(
1288 4 : "Pratyush and Parrinello, Phys. Rev. Lett. 111, 230602 (2013)");
1289 153 : if(calc_rct_) log<<plumed.cite(
1290 5 : "Pratyush and Parrinello, J. Phys. Chem. B, 119, 736 (2015)");
1291 438 : if(concurrent || rect_biasf_.size()>0) log<<plumed.cite(
1292 78 : "Gil-Ley and Bussi, J. Chem. Theory Comput. 11, 1077 (2015)");
1293 174 : if(rect_biasf_.size()>0 && walkers_mpi) log<<plumed.cite(
1294 12 : "Hosek, Toulcova, Bortolato, and Spiwok, J. Phys. Chem. B 120, 2209 (2016)");
1295 138 : if(targetfilename_.length()>0) {
1296 6 : log<<plumed.cite("White, Dama, and Voth, J. Chem. Theory Comput. 11, 2451 (2015)");
1297 6 : log<<plumed.cite("Marinelli and Faraldo-Gómez, Biophys. J. 108, 2779 (2015)");
1298 6 : log<<plumed.cite("Gil-Ley, Bottaro, and Bussi, J. Chem. Theory Comput. 12, 2790 (2016)");
1299 : }
1300 138 : if(freq_adaptive_) {
1301 6 : log<<plumed.cite("Wang, Valsson, Tiwary, Parrinello, and Lindorff-Larsen, J. Chem. Phys. 149, 072309 (2018)");
1302 : }
1303 138 : log<<"\n";
1304 138 : }
1305 :
1306 139 : void MetaD::readTemperingSpecs(TemperingSpecs &t_specs) {
1307 : // Set global tempering parameters.
1308 278 : parse(t_specs.name_stem + "BIASFACTOR", t_specs.biasf);
1309 139 : if (t_specs.biasf != -1.0) {
1310 3 : if (kbt_ == 0.0) {
1311 0 : error("Unless the MD engine passes the temperature to plumed, with tempered metad you must specify it using TEMP");
1312 : }
1313 3 : if (t_specs.biasf == 1.0) {
1314 0 : error("A bias factor of 1 corresponds to zero delta T and zero hill size, so it is not allowed.");
1315 : }
1316 3 : t_specs.is_active = true;
1317 6 : parse(t_specs.name_stem + "BIASTHRESHOLD", t_specs.threshold);
1318 3 : if (t_specs.threshold < 0.0) {
1319 0 : error(t_specs.name + " bias threshold is nonsensical");
1320 : }
1321 6 : parse(t_specs.name_stem + "ALPHA", t_specs.alpha);
1322 3 : if (t_specs.alpha <= 0.0 || t_specs.alpha > 1.0) {
1323 0 : error(t_specs.name + " decay shape parameter alpha is nonsensical");
1324 : }
1325 : }
1326 139 : }
1327 :
1328 3 : void MetaD::logTemperingSpecs(const TemperingSpecs &t_specs) {
1329 6 : log.printf(" %s bias factor %f\n", t_specs.name.c_str(), t_specs.biasf);
1330 3 : log.printf(" KbT %f\n", kbt_);
1331 5 : if (t_specs.threshold != 0.0) log.printf(" %s bias threshold %f\n", t_specs.name.c_str(), t_specs.threshold);
1332 4 : if (t_specs.alpha != 1.0) log.printf(" %s decay shape parameter alpha %f\n", t_specs.name.c_str(), t_specs.alpha);
1333 3 : }
1334 :
1335 6036 : void MetaD::readGaussians(IFile *ifile)
1336 : {
1337 6036 : unsigned ncv=getNumberOfArguments();
1338 6036 : vector<double> center(ncv);
1339 6036 : vector<double> sigma(ncv);
1340 : double height;
1341 : int nhills=0;
1342 6036 : bool multivariate=false;
1343 :
1344 6036 : std::vector<Value> tmpvalues;
1345 42291 : for(unsigned j=0; j<getNumberOfArguments(); ++j) tmpvalues.push_back( Value( this, getPntrToArgument(j)->getName(), false ) );
1346 :
1347 11758 : while(scanOneHill(ifile,tmpvalues,center,sigma,height,multivariate)) {
1348 : ;
1349 2861 : nhills++;
1350 : // note that for gamma=1 we store directly -F
1351 2861 : if(welltemp_ && biasf_>1.0) {height*=(biasf_-1.0)/biasf_;}
1352 2861 : addGaussian(Gaussian(center,sigma,height,multivariate));
1353 : }
1354 6036 : log.printf(" %d Gaussians read\n",nhills);
1355 6036 : }
1356 :
1357 2902 : void MetaD::writeGaussian(const Gaussian& hill, OFile&file)
1358 : {
1359 2902 : unsigned ncv=getNumberOfArguments();
1360 5804 : file.printField("time",getTimeStep()*getStep());
1361 13374 : for(unsigned i=0; i<ncv; ++i) {
1362 10472 : file.printField(getPntrToArgument(i),hill.center[i]);
1363 : }
1364 8706 : hillsOfile_.printField("kerneltype","gaussian");
1365 2902 : if(hill.multivariate) {
1366 1338 : hillsOfile_.printField("multivariate","true");
1367 : Matrix<double> mymatrix(ncv,ncv);
1368 : unsigned k=0;
1369 1648 : for(unsigned i=0; i<ncv; i++) {
1370 2113 : for(unsigned j=i; j<ncv; j++) {
1371 : // recompose the full inverse matrix
1372 2268 : mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k];
1373 756 : k++;
1374 : }
1375 : }
1376 : // invert it
1377 : Matrix<double> invmatrix(ncv,ncv);
1378 446 : Invert(mymatrix,invmatrix);
1379 : // enforce symmetry
1380 1648 : for(unsigned i=0; i<ncv; i++) {
1381 2113 : for(unsigned j=i; j<ncv; j++) {
1382 756 : invmatrix(i,j)=invmatrix(j,i);
1383 : }
1384 : }
1385 :
1386 : // do cholesky so to have a "sigma like" number
1387 : Matrix<double> lower(ncv,ncv);
1388 446 : cholesky(invmatrix,lower);
1389 : // loop in band form
1390 1648 : for(unsigned i=0; i<ncv; i++) {
1391 2113 : for(unsigned j=0; j<ncv-i; j++) {
1392 4536 : file.printField("sigma_"+getPntrToArgument(j+i)->getName()+"_"+getPntrToArgument(j)->getName(),lower(j+i,j));
1393 : }
1394 : }
1395 : } else {
1396 7368 : hillsOfile_.printField("multivariate","false");
1397 11726 : for(unsigned i=0; i<ncv; ++i)
1398 18540 : file.printField("sigma_"+getPntrToArgument(i)->getName(),hill.sigma[i]);
1399 : }
1400 2902 : double height=hill.height;
1401 : // note that for gamma=1 we store directly -F
1402 2902 : if(welltemp_ && biasf_>1.0) height*=biasf_/(biasf_-1.0);
1403 8706 : file.printField("height",height).printField("biasf",biasf_);
1404 4411 : if(mw_n_>1) file.printField("clock",int(std::time(0)));
1405 2902 : file.printField();
1406 2902 : }
1407 :
1408 5925 : void MetaD::addGaussian(const Gaussian& hill)
1409 : {
1410 5925 : if(!grid_) hills_.push_back(hill);
1411 : else {
1412 626 : unsigned ncv=getNumberOfArguments();
1413 626 : vector<unsigned> nneighb=getGaussianSupport(hill);
1414 626 : vector<Grid::index_t> neighbors=BiasGrid_->getNeighbors(hill.center,nneighb);
1415 626 : vector<double> der(ncv);
1416 626 : vector<double> xx(ncv);
1417 626 : if(comm.Get_size()==1) {
1418 164534 : for(unsigned i=0; i<neighbors.size(); ++i) {
1419 54486 : Grid::index_t ineigh=neighbors[i];
1420 261594 : for(unsigned j=0; j<ncv; ++j) der[j]=0.0;
1421 54486 : BiasGrid_->getPoint(ineigh,xx);
1422 54486 : double bias=evaluateGaussian(xx,hill,&der[0]);
1423 54486 : BiasGrid_->addValueAndDerivatives(ineigh,bias,der);
1424 : }
1425 : } else {
1426 88 : unsigned stride=comm.Get_size();
1427 88 : unsigned rank=comm.Get_rank();
1428 176 : vector<double> allder(ncv*neighbors.size(),0.0);
1429 176 : vector<double> allbias(neighbors.size(),0.0);
1430 81224 : for(unsigned i=rank; i<neighbors.size(); i+=stride) {
1431 27016 : Grid::index_t ineigh=neighbors[i];
1432 27016 : BiasGrid_->getPoint(ineigh,xx);
1433 54032 : allbias[i]=evaluateGaussian(xx,hill,&allder[ncv*i]);
1434 : }
1435 88 : comm.Sum(allbias);
1436 88 : comm.Sum(allder);
1437 309272 : for(unsigned i=0; i<neighbors.size(); ++i) {
1438 103032 : Grid::index_t ineigh=neighbors[i];
1439 721224 : for(unsigned j=0; j<ncv; ++j) {der[j]=allder[ncv*i+j];}
1440 206064 : BiasGrid_->addValueAndDerivatives(ineigh,allbias[i],der);
1441 : }
1442 : }
1443 : }
1444 5925 : }
1445 :
1446 626 : vector<unsigned> MetaD::getGaussianSupport(const Gaussian& hill)
1447 : {
1448 : vector<unsigned> nneigh;
1449 : vector<double> cutoff;
1450 626 : unsigned ncv=getNumberOfArguments();
1451 :
1452 : // traditional or flexible hill?
1453 626 : if(hill.multivariate) {
1454 : unsigned k=0;
1455 : Matrix<double> mymatrix(ncv,ncv);
1456 0 : for(unsigned i=0; i<ncv; i++) {
1457 0 : for(unsigned j=i; j<ncv; j++) {
1458 : // recompose the full inverse matrix
1459 0 : mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k];
1460 0 : k++;
1461 : }
1462 : }
1463 : // Reinvert so to have the ellipses
1464 : Matrix<double> myinv(ncv,ncv);
1465 0 : Invert(mymatrix,myinv);
1466 : Matrix<double> myautovec(ncv,ncv);
1467 0 : vector<double> myautoval(ncv); //should I take this or their square root?
1468 0 : diagMat(myinv,myautoval,myautovec);
1469 : double maxautoval=0.;
1470 : unsigned ind_maxautoval; ind_maxautoval=ncv;
1471 0 : for(unsigned i=0; i<ncv; i++) {
1472 0 : if(myautoval[i]>maxautoval) {maxautoval=myautoval[i]; ind_maxautoval=i;}
1473 : }
1474 0 : for(unsigned i=0; i<ncv; i++) {
1475 0 : cutoff.push_back(sqrt(2.0*DP2CUTOFF)*abs(sqrt(maxautoval)*myautovec(i,ind_maxautoval)));
1476 : }
1477 : } else {
1478 2526 : for(unsigned i=0; i<ncv; ++i) {
1479 2850 : cutoff.push_back(sqrt(2.0*DP2CUTOFF)*hill.sigma[i]);
1480 : }
1481 : }
1482 :
1483 626 : if(doInt_) {
1484 4 : if(hill.center[0]+cutoff[0] > uppI_ || hill.center[0]-cutoff[0] < lowI_) {
1485 : // in this case, we updated the entire grid to avoid problems
1486 2 : return BiasGrid_->getNbin();
1487 : } else {
1488 0 : nneigh.push_back( static_cast<unsigned>(ceil(cutoff[0]/BiasGrid_->getDx()[0])) );
1489 : return nneigh;
1490 : }
1491 : } else {
1492 2520 : for(unsigned i=0; i<ncv; i++) {
1493 5688 : nneigh.push_back( static_cast<unsigned>(ceil(cutoff[i]/BiasGrid_->getDx()[i])) );
1494 : }
1495 : }
1496 :
1497 : return nneigh;
1498 : }
1499 :
1500 19032 : double MetaD::getBiasAndDerivatives(const vector<double>& cv, double* der)
1501 : {
1502 19032 : double bias=0.0;
1503 19032 : if(!grid_) {
1504 14767 : if(hills_.size()>10000 && (getStep()-last_step_warn_grid)>10000) {
1505 : std::string msg;
1506 0 : Tools::convert(hills_.size(),msg);
1507 0 : msg="You have accumulated "+msg+" hills, you should enable GRIDs to avoid serious performance hits";
1508 0 : warning(msg);
1509 0 : last_step_warn_grid=getStep();
1510 : }
1511 14767 : unsigned stride=comm.Get_size();
1512 14767 : unsigned rank=comm.Get_rank();
1513 20909522 : for(unsigned i=rank; i<hills_.size(); i+=stride) {
1514 6959996 : bias+=evaluateGaussian(cv,hills_[i],der);
1515 : }
1516 14767 : comm.Sum(bias);
1517 14767 : if(der) comm.Sum(der,getNumberOfArguments());
1518 : } else {
1519 4265 : if(der) {
1520 1356 : vector<double> vder(getNumberOfArguments());
1521 1356 : bias=BiasGrid_->getValueAndDerivatives(cv,vder);
1522 4980 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {der[i]=vder[i];}
1523 : } else {
1524 2909 : bias = BiasGrid_->getValue(cv);
1525 : }
1526 : }
1527 :
1528 19032 : return bias;
1529 : }
1530 :
1531 0 : double MetaD::getGaussianNormalization( const Gaussian& hill )
1532 : {
1533 : double norm=1;
1534 0 : unsigned ncv=hill.center.size();
1535 :
1536 0 : if(hill.multivariate) {
1537 : // recompose the full sigma from the upper diag cholesky
1538 : unsigned k=0;
1539 : Matrix<double> mymatrix(ncv,ncv);
1540 0 : for(unsigned i=0; i<ncv; i++) {
1541 0 : for(unsigned j=i; j<ncv; j++) {
1542 0 : mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k]; // recompose the full inverse matrix
1543 0 : k++;
1544 : }
1545 0 : double ldet; logdet( mymatrix, ldet );
1546 0 : norm = exp( ldet ); // Not sure here if mymatrix is sigma or inverse
1547 : }
1548 : } else {
1549 0 : for(unsigned i=0; i<hill.sigma.size(); ++i) norm*=hill.sigma[i];
1550 : }
1551 :
1552 0 : return norm*pow(2*pi,static_cast<double>(ncv)/2.0);
1553 : }
1554 :
1555 7041498 : double MetaD::evaluateGaussian(const vector<double>& cv, const Gaussian& hill, double* der)
1556 : {
1557 : double dp2=0.0;
1558 : double bias=0.0;
1559 : // I use a pointer here because cv is const (and should be const)
1560 : // but when using doInt it is easier to locally replace cv[0] with
1561 : // the upper/lower limit in case it is out of range
1562 : const double *pcv=NULL; // pointer to cv
1563 : double tmpcv[1]; // tmp array with cv (to be used with doInt_)
1564 7041498 : if(cv.size()>0) pcv=&cv[0];
1565 7041498 : if(doInt_) {
1566 1402 : plumed_assert(cv.size()==1);
1567 1402 : tmpcv[0]=cv[0];
1568 1402 : if(cv[0]<lowI_) tmpcv[0]=lowI_;
1569 1402 : if(cv[0]>uppI_) tmpcv[0]=uppI_;
1570 : pcv=&(tmpcv[0]);
1571 : }
1572 7041498 : if(hill.multivariate) {
1573 : unsigned k=0;
1574 230564 : unsigned ncv=cv.size();
1575 : // recompose the full sigma from the upper diag cholesky
1576 : Matrix<double> mymatrix(ncv,ncv);
1577 694540 : for(unsigned i=0; i<ncv; i++) {
1578 698812 : for(unsigned j=i; j<ncv; j++) {
1579 700236 : mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k]; // recompose the full inverse matrix
1580 233412 : k++;
1581 : }
1582 : }
1583 1157092 : for(unsigned i=0; i<cv.size(); ++i) {
1584 463976 : double dp_i=difference(i,hill.center[i],pcv[i]);
1585 231988 : dp_[i]=dp_i;
1586 1164212 : for(unsigned j=i; j<cv.size(); ++j) {
1587 233412 : if(i==j) {
1588 463976 : dp2+=dp_i*dp_i*mymatrix(i,j)*0.5;
1589 : } else {
1590 2848 : double dp_j=difference(j,hill.center[j],pcv[j]);
1591 2848 : dp2+=dp_i*dp_j*mymatrix(i,j);
1592 : }
1593 : }
1594 : }
1595 230564 : if(dp2<DP2CUTOFF) {
1596 221813 : bias=hill.height*exp(-dp2);
1597 221813 : if(der) {
1598 389336 : for(unsigned i=0; i<cv.size(); ++i) {
1599 : double tmp=0.0;
1600 235198 : for(unsigned j=0; j<cv.size(); ++j) {
1601 157208 : tmp += dp_[j]*mymatrix(i,j)*bias;
1602 : }
1603 77990 : der[i]-=tmp;
1604 : }
1605 : }
1606 : }
1607 : } else {
1608 54463568 : for(unsigned i=0; i<cv.size(); ++i) {
1609 40841700 : double dp=difference(i,hill.center[i],pcv[i])*hill.invsigma[i];
1610 13613900 : dp2+=dp*dp;
1611 13613900 : dp_[i]=dp;
1612 : }
1613 6810934 : dp2*=0.5;
1614 6810934 : if(dp2<DP2CUTOFF) {
1615 3947846 : bias=hill.height*exp(-dp2);
1616 3947846 : if(der) {
1617 13510290 : for(unsigned i=0; i<cv.size(); ++i) {der[i]+=-bias*dp_[i]*hill.invsigma[i];}
1618 : }
1619 : }
1620 : }
1621 :
1622 7041498 : if(doInt_ && der) {
1623 1558 : if(cv[0]<lowI_ || cv[0]>uppI_) for(unsigned i=0; i<cv.size(); ++i) der[i]=0;
1624 : }
1625 :
1626 7041498 : return bias;
1627 : }
1628 :
1629 2716 : double MetaD::getHeight(const vector<double>& cv)
1630 : {
1631 2716 : double height=height0_;
1632 2716 : if(welltemp_) {
1633 269 : double vbias = getBiasAndDerivatives(cv);
1634 269 : if(biasf_>1.0) {
1635 253 : height = height0_*exp(-vbias/(kbt_*(biasf_-1.0)));
1636 : } else {
1637 : // notice that if gamma=1 we store directly -F
1638 16 : height = height0_*exp(-vbias/kbt_);
1639 : }
1640 : }
1641 2716 : if(dampfactor_>0.0) {
1642 18 : plumed_assert(BiasGrid_);
1643 18 : double m=BiasGrid_->getMaxValue();
1644 18 : height*=exp(-m/(kbt_*(dampfactor_)));
1645 : }
1646 2716 : if (tt_specs_.is_active) {
1647 60 : double vbarrier = transition_bias_;
1648 60 : temperHeight(height, tt_specs_, vbarrier);
1649 : }
1650 2716 : if(TargetGrid_) {
1651 36 : double f=TargetGrid_->getValue(cv)-TargetGrid_->getMaxValue();
1652 18 : height*=exp(f/kbt_);
1653 : }
1654 2716 : return height;
1655 : }
1656 :
1657 60 : void MetaD::temperHeight(double &height, const TemperingSpecs &t_specs, const double tempering_bias) {
1658 60 : if (t_specs.alpha == 1.0) {
1659 80 : height *= exp(-max(0.0, tempering_bias - t_specs.threshold) / (kbt_ * (t_specs.biasf - 1.0)));
1660 : } else {
1661 40 : height *= pow(1 + (1 - t_specs.alpha) / t_specs.alpha * max(0.0, tempering_bias - t_specs.threshold) / (kbt_ * (t_specs.biasf - 1.0)), - t_specs.alpha / (1 - t_specs.alpha));
1662 : }
1663 60 : }
1664 :
1665 6605 : void MetaD::calculate()
1666 : {
1667 : // this is because presently there is no way to properly pass information
1668 : // on adaptive hills (diff) after exchanges:
1669 6605 : if(adaptive_==FlexibleBin::diffusion && getExchangeStep()) error("ADAPTIVE=DIFF is not compatible with replica exchange");
1670 :
1671 6605 : const unsigned ncv=getNumberOfArguments();
1672 6605 : vector<double> cv(ncv);
1673 6605 : std::unique_ptr<double[]> der(new double[ncv]);
1674 28169 : for(unsigned i=0; i<ncv; ++i) {
1675 21564 : cv[i]=getArgument(i);
1676 10782 : der[i]=0.;
1677 : }
1678 6605 : double ene = getBiasAndDerivatives(cv,der.get());
1679 : // special case for gamma=1.0
1680 6605 : if(biasf_==1.0) {
1681 : ene=0.0;
1682 200 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {der[i]=0.0;}
1683 : }
1684 :
1685 : setBias(ene);
1686 6710 : if(calc_rct_) getPntrToComponent("rbias")->set(ene - reweight_factor_);
1687 : // calculate the acceleration factor
1688 6605 : if(acceleration&&!isFirstStep) {
1689 329 : acc += static_cast<double>(getStride()) * exp(ene/(kbt_));
1690 329 : const double mean_acc = acc/((double) getStep());
1691 658 : getPntrToComponent("acc")->set(mean_acc);
1692 6276 : } else if (acceleration && isFirstStep && acc_restart_mean_ > 0.0) {
1693 2 : acc = acc_restart_mean_ * static_cast<double>(getStep());
1694 2 : if(freq_adaptive_) {
1695 : // has to be done here if restarting, as the acc is not defined before
1696 1 : updateFrequencyAdaptiveStride();
1697 : }
1698 : }
1699 :
1700 13210 : getPntrToComponent("work")->set(work_);
1701 : // set Forces
1702 28169 : for(unsigned i=0; i<ncv; ++i) {
1703 21564 : setOutputForce(i,-der[i]);
1704 : }
1705 6605 : }
1706 :
1707 6079 : void MetaD::update() {
1708 6079 : vector<double> cv(getNumberOfArguments());
1709 : vector<double> thissigma;
1710 : bool multivariate;
1711 :
1712 : // adding hills criteria (could be more complex though)
1713 : bool nowAddAHill;
1714 6079 : if(getStep()%current_stride==0 && !isFirstStep )nowAddAHill=true;
1715 : else {
1716 : nowAddAHill=false;
1717 3363 : isFirstStep=false;
1718 : }
1719 :
1720 42926 : for(unsigned i=0; i<cv.size(); ++i) cv[i] = getArgument(i);
1721 :
1722 6079 : double vbias=getBiasAndDerivatives(cv);
1723 :
1724 : // if you use adaptive, call the FlexibleBin
1725 6079 : if(adaptive_!=FlexibleBin::none) {
1726 778 : flexbin->update(nowAddAHill);
1727 : multivariate=true;
1728 : } else {
1729 : multivariate=false;
1730 : }
1731 :
1732 6079 : if(nowAddAHill) {
1733 : // add a Gaussian
1734 2716 : double height=getHeight(cv);
1735 : // returns upper diagonal inverse
1736 3464 : if(adaptive_!=FlexibleBin::none) thissigma=flexbin->getInverseMatrix();
1737 : // returns normal sigma
1738 2342 : else thissigma=sigma0_;
1739 :
1740 : // In case we use walkers_mpi, it is now necessary to communicate with other replicas.
1741 2716 : if(walkers_mpi) {
1742 : // Allocate arrays to store all walkers hills
1743 348 : std::vector<double> all_cv(mpi_nw_*cv.size(),0.0);
1744 348 : std::vector<double> all_sigma(mpi_nw_*thissigma.size(),0.0);
1745 174 : std::vector<double> all_height(mpi_nw_,0.0);
1746 174 : std::vector<int> all_multivariate(mpi_nw_,0);
1747 174 : if(comm.Get_rank()==0) {
1748 : // Communicate (only root)
1749 99 : multi_sim_comm.Allgather(cv,all_cv);
1750 99 : multi_sim_comm.Allgather(thissigma,all_sigma);
1751 : // notice that if gamma=1 we store directly -F so this scaling is not necessary:
1752 99 : multi_sim_comm.Allgather(height*(biasf_>1.0?biasf_/(biasf_-1.0):1.0),all_height);
1753 99 : multi_sim_comm.Allgather(int(multivariate),all_multivariate);
1754 : }
1755 : // Share info with group members
1756 174 : comm.Bcast(all_cv,0);
1757 174 : comm.Bcast(all_sigma,0);
1758 174 : comm.Bcast(all_height,0);
1759 174 : comm.Bcast(all_multivariate,0);
1760 :
1761 : // Flying Gaussian
1762 174 : if (flying) {
1763 : hills_.clear();
1764 54 : comm.Barrier();
1765 : }
1766 :
1767 1218 : for(unsigned i=0; i<mpi_nw_; i++) {
1768 : // actually add hills one by one
1769 522 : std::vector<double> cv_now(cv.size());
1770 522 : std::vector<double> sigma_now(thissigma.size());
1771 4176 : for(unsigned j=0; j<cv.size(); j++) cv_now[j]=all_cv[i*cv.size()+j];
1772 4500 : for(unsigned j=0; j<thissigma.size(); j++) sigma_now[j]=all_sigma[i*thissigma.size()+j];
1773 : // notice that if gamma=1 we store directly -F so this scaling is not necessary:
1774 2088 : Gaussian newhill=Gaussian(cv_now,sigma_now,all_height[i]*(biasf_>1.0?(biasf_-1.0)/biasf_:1.0),all_multivariate[i]);
1775 522 : addGaussian(newhill);
1776 :
1777 : // Flying Gaussian
1778 522 : if (!flying) {
1779 360 : writeGaussian(newhill,hillsOfile_);
1780 : }
1781 :
1782 : }
1783 : } else {
1784 5084 : Gaussian newhill=Gaussian(cv,thissigma,height,multivariate);
1785 2542 : addGaussian(newhill);
1786 : // print on HILLS file
1787 2542 : writeGaussian(newhill,hillsOfile_);
1788 : }
1789 : }
1790 :
1791 : // this should be outside of the if block in case
1792 : // mw_rstride_ is not a multiple of stride_
1793 6079 : if(mw_n_>1 && getStep()%mw_rstride_==0) {
1794 3012 : hillsOfile_.flush();
1795 : }
1796 :
1797 6079 : double vbias1=getBiasAndDerivatives(cv);
1798 6079 : work_+=vbias1-vbias;
1799 :
1800 : // dump grid on file
1801 6079 : if(wgridstride_>0&&(getStep()%wgridstride_==0||getCPT())) {
1802 : // in case old grids are stored, a sequence of grids should appear
1803 : // this call results in a repetition of the header:
1804 80 : if(storeOldGrids_) gridfile_.clearFields();
1805 : // in case only latest grid is stored, file should be rewound
1806 : // this will overwrite previously written grids
1807 : else {
1808 40 : int r = 0;
1809 40 : if(walkers_mpi) {
1810 0 : if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
1811 0 : comm.Bcast(r,0);
1812 : }
1813 40 : if(r==0) gridfile_.rewind();
1814 : }
1815 80 : BiasGrid_->writeToFile(gridfile_);
1816 : // if a single grid is stored, it is necessary to flush it, otherwise
1817 : // the file might stay empty forever (when a single grid is not large enough to
1818 : // trigger flushing from the operating system).
1819 : // on the other hand, if grids are stored one after the other this is
1820 : // no necessary, and we leave the flushing control to the user as usual
1821 : // (with FLUSH keyword)
1822 80 : if(!storeOldGrids_) gridfile_.flush();
1823 : }
1824 :
1825 : // if multiple walkers and time to read Gaussians
1826 6079 : if(mw_n_>1 && getStep()%mw_rstride_==0) {
1827 21084 : for(int i=0; i<mw_n_; ++i) {
1828 : // don't read your own Gaussians
1829 9036 : if(i==mw_id_) continue;
1830 : // if the file is not open yet
1831 12048 : if(!(ifiles[i]->isOpen())) {
1832 : // check if it exists now and open it!
1833 12 : if(ifiles[i]->FileExist(ifilesnames[i])) {
1834 12 : ifiles[i]->open(ifilesnames[i]);
1835 6 : ifiles[i]->reset(false);
1836 : }
1837 : // otherwise read the new Gaussians
1838 : } else {
1839 12036 : log.printf(" Reading hills from %s:",ifilesnames[i].c_str());
1840 6018 : readGaussians(ifiles[i].get());
1841 6018 : ifiles[i]->reset(false);
1842 : }
1843 : }
1844 : }
1845 : // Recalculate special bias quantities whenever the bias has been changed by the update.
1846 6079 : bool bias_has_changed = (nowAddAHill || (mw_n_ > 1 && getStep() % mw_rstride_ == 0));
1847 6079 : if (calc_rct_ && bias_has_changed && getStep()%(stride_*rct_ustride_)==0) computeReweightingFactor();
1848 6079 : if (calc_max_bias_ && bias_has_changed) {
1849 0 : max_bias_ = BiasGrid_->getMaxValue();
1850 0 : getPntrToComponent("maxbias")->set(max_bias_);
1851 : }
1852 6079 : if (calc_transition_bias_ && bias_has_changed) {
1853 260 : transition_bias_ = getTransitionBarrierBias();
1854 520 : getPntrToComponent("transbias")->set(transition_bias_);
1855 : }
1856 :
1857 : // Frequency adaptive metadynamics - update hill addition frequency
1858 6079 : if(freq_adaptive_ && getStep()%fa_update_frequency_==0) {
1859 151 : updateFrequencyAdaptiveStride();
1860 : }
1861 :
1862 6079 : }
1863 :
1864 : /// takes a pointer to the file and a template string with values v and gives back the next center, sigma and height
1865 8897 : bool MetaD::scanOneHill(IFile *ifile, vector<Value> &tmpvalues, vector<double> ¢er, vector<double> &sigma, double &height, bool &multivariate)
1866 : {
1867 : double dummy;
1868 8897 : multivariate=false;
1869 17794 : if(ifile->scanField("time",dummy)) {
1870 2861 : unsigned ncv; ncv=tmpvalues.size();
1871 14213 : for(unsigned i=0; i<ncv; ++i) {
1872 11352 : ifile->scanField( &tmpvalues[i] );
1873 5676 : if( tmpvalues[i].isPeriodic() && ! getPntrToArgument(i)->isPeriodic() ) {
1874 0 : error("in hills file periodicity for variable " + tmpvalues[i].getName() + " does not match periodicity in input");
1875 5676 : } else if( tmpvalues[i].isPeriodic() ) {
1876 0 : std::string imin, imax; tmpvalues[i].getDomain( imin, imax );
1877 0 : std::string rmin, rmax; getPntrToArgument(i)->getDomain( rmin, rmax );
1878 0 : if( imin!=rmin || imax!=rmax ) {
1879 0 : error("in hills file periodicity for variable " + tmpvalues[i].getName() + " does not match periodicity in input");
1880 : }
1881 : }
1882 5676 : center[i]=tmpvalues[i].get();
1883 : }
1884 : // scan for kerneltype
1885 2861 : std::string ktype="gaussian";
1886 8583 : if( ifile->FieldExist("kerneltype") ) ifile->scanField("kerneltype",ktype);
1887 : // scan for multivariate label: record the actual file position so to eventually rewind
1888 : std::string sss;
1889 5722 : ifile->scanField("multivariate",sss);
1890 2861 : if(sss=="true") multivariate=true;
1891 2861 : else if(sss=="false") multivariate=false;
1892 0 : else plumed_merror("cannot parse multivariate = "+ sss);
1893 2861 : if(multivariate) {
1894 0 : sigma.resize(ncv*(ncv+1)/2);
1895 : Matrix<double> upper(ncv,ncv);
1896 : Matrix<double> lower(ncv,ncv);
1897 0 : for(unsigned i=0; i<ncv; i++) {
1898 0 : for(unsigned j=0; j<ncv-i; j++) {
1899 0 : ifile->scanField("sigma_"+getPntrToArgument(j+i)->getName()+"_"+getPntrToArgument(j)->getName(),lower(j+i,j));
1900 0 : upper(j,j+i)=lower(j+i,j);
1901 : }
1902 : }
1903 : Matrix<double> mymult(ncv,ncv);
1904 : Matrix<double> invmatrix(ncv,ncv);
1905 0 : mult(lower,upper,mymult);
1906 : // now invert and get the sigmas
1907 0 : Invert(mymult,invmatrix);
1908 : // put the sigmas in the usual order: upper diagonal (this time in normal form and not in band form)
1909 : unsigned k=0;
1910 0 : for(unsigned i=0; i<ncv; i++) {
1911 0 : for(unsigned j=i; j<ncv; j++) {
1912 0 : sigma[k]=invmatrix(i,j);
1913 0 : k++;
1914 : }
1915 : }
1916 : } else {
1917 14213 : for(unsigned i=0; i<ncv; ++i) {
1918 17028 : ifile->scanField("sigma_"+getPntrToArgument(i)->getName(),sigma[i]);
1919 : }
1920 : }
1921 :
1922 5722 : ifile->scanField("height",height);
1923 5722 : ifile->scanField("biasf",dummy);
1924 8409 : if(ifile->FieldExist("clock")) ifile->scanField("clock",dummy);
1925 5722 : if(ifile->FieldExist("lower_int")) ifile->scanField("lower_int",dummy);
1926 5722 : if(ifile->FieldExist("upper_int")) ifile->scanField("upper_int",dummy);
1927 2861 : ifile->scanField();
1928 : return true;
1929 : } else {
1930 : return false;
1931 : }
1932 : }
1933 :
1934 100 : void MetaD::computeReweightingFactor()
1935 : {
1936 100 : if(biasf_==1.0) { // in this case we have no bias, so reweight factor is 0.0
1937 0 : getPntrToComponent("rct")->set(0.0);
1938 0 : return;
1939 : }
1940 :
1941 100 : double Z_0=0; //proportional to the integral of exp(-beta*F)
1942 100 : double Z_V=0; //proportional to the integral of exp(-beta*(F+V))
1943 100 : double minusBetaF=biasf_/(biasf_-1.)/kbt_;
1944 100 : double minusBetaFplusV=1./(biasf_-1.)/kbt_;
1945 100 : if (biasf_==-1.0) { //non well-tempered case
1946 0 : minusBetaF=1./kbt_;
1947 : minusBetaFplusV=0;
1948 : }
1949 100 : max_bias_=BiasGrid_->getMaxValue(); //to avoid exp overflow
1950 :
1951 100 : const unsigned rank=comm.Get_rank();
1952 100 : const unsigned stride=comm.Get_size();
1953 1800200 : for (Grid::index_t t=rank; t<BiasGrid_->getSize(); t+=stride) {
1954 900000 : const double val=BiasGrid_->getValue(t);
1955 900000 : Z_0+=std::exp(minusBetaF*(val-max_bias_));
1956 900000 : Z_V+=std::exp(minusBetaFplusV*(val-max_bias_));
1957 : }
1958 100 : if (stride>1) {
1959 80 : comm.Sum(Z_0);
1960 80 : comm.Sum(Z_V);
1961 : }
1962 :
1963 100 : reweight_factor_=kbt_*std::log(Z_0/Z_V)+max_bias_;
1964 200 : getPntrToComponent("rct")->set(reweight_factor_);
1965 : }
1966 :
1967 273 : double MetaD::getTransitionBarrierBias() {
1968 :
1969 : // If there is only one well of interest, return the bias at that well point.
1970 273 : if (transitionwells_.size() == 1) {
1971 0 : double tb_bias = getBiasAndDerivatives(transitionwells_[0], NULL);
1972 0 : return tb_bias;
1973 :
1974 : // Otherwise, check for the least barrier bias between all pairs of wells.
1975 : // Note that because the paths can be considered edges between the wells' nodes
1976 : // to make a graph and the path barriers satisfy certain cycle inequalities, it
1977 : // is sufficient to look at paths corresponding to a minimal spanning tree of the
1978 : // overall graph rather than examining every edge in the graph.
1979 : // For simplicity, I chose the star graph with center well 0 as the spanning tree.
1980 : // It is most efficient to start the path searches from the wells that are
1981 : // expected to be sampled last, so transitionwell_[0] should correspond to the
1982 : // starting well. With this choice the searches will terminate in one step until
1983 : // transitionwell_[1] is sampled.
1984 : } else {
1985 : double least_transition_bias;
1986 273 : vector<double> sink = transitionwells_[0];
1987 273 : vector<double> source = transitionwells_[1];
1988 273 : least_transition_bias = BiasGrid_->findMaximalPathMinimum(source, sink);
1989 546 : for (unsigned i = 2; i < transitionwells_.size(); i++) {
1990 0 : if (least_transition_bias == 0.0) {
1991 : break;
1992 : }
1993 0 : source = transitionwells_[i];
1994 0 : double curr_transition_bias = BiasGrid_->findMaximalPathMinimum(source, sink);
1995 0 : least_transition_bias = fmin(curr_transition_bias, least_transition_bias);
1996 : }
1997 : return least_transition_bias;
1998 : }
1999 : }
2000 :
2001 154 : void MetaD::updateFrequencyAdaptiveStride() {
2002 154 : plumed_massert(freq_adaptive_,"should only be used if frequency adaptive metadynamics is enabled");
2003 154 : plumed_massert(acceleration,"frequency adaptive metadynamics can only be used if the acceleration factor is calculated");
2004 154 : const double mean_acc = acc/((double) getStep());
2005 154 : int tmp_stride= stride_*floor((mean_acc/fa_min_acceleration_)+0.5);
2006 154 : if(mean_acc >= fa_min_acceleration_) {
2007 129 : if(tmp_stride > current_stride) {current_stride = tmp_stride;}
2008 : }
2009 154 : if(fa_max_stride_!=0 && current_stride>fa_max_stride_) {
2010 0 : current_stride=fa_max_stride_;
2011 : }
2012 308 : getPntrToComponent("pace")->set(current_stride);
2013 154 : }
2014 :
2015 6605 : bool MetaD::checkNeedsGradients()const
2016 : {
2017 6605 : if(adaptive_==FlexibleBin::geometry) {
2018 192 : if(getStep()%stride_==0 && !isFirstStep) return true;
2019 : else return false;
2020 : } else return false;
2021 : }
2022 :
2023 : }
2024 5517 : }
|
