1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2011-2021 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 class MetaD : public Bias {
375
376 private:
377 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;
GaussianPLMD::bias::MetaD::Gaussian383 Gaussian(const vector<double> & center,const vector<double> & sigma,double height, bool multivariate ):
384 center(center),sigma(sigma),height(height),multivariate(multivariate),invsigma(sigma) {
385 // to avoid troubles from zero element in flexible hills
386 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 }
388 };
389 struct TemperingSpecs {
390 bool is_active;
391 std::string name_stem;
392 std::string name;
393 double biasf;
394 double threshold;
395 double alpha;
TemperingSpecsPLMD::bias::MetaD::TemperingSpecs396 inline TemperingSpecs(bool is_active, const std::string &name_stem, const std::string &name, double biasf, double threshold, double alpha) :
397 is_active(is_active), name_stem(name_stem), name(name), biasf(biasf), threshold(threshold), alpha(alpha)
398 {}
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<GridBase> 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<GridBase> 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() override;
480 void update() override;
481 static void registerKeywords(Keywords& keys);
482 bool checkNeedsGradients()const override;
483 };
484
485 PLUMED_REGISTER_ACTION(MetaD,"METAD")
486
registerKeywords(Keywords & keys)487 void MetaD::registerKeywords(Keywords& keys) {
488 Bias::registerKeywords(keys);
489 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 keys.addOutputComponent("rct","CALC_RCT","the reweighting factor \\f$c(t)\\f$.");
492 keys.addOutputComponent("work","default","accumulator for work");
493 keys.addOutputComponent("acc","ACCELERATION","the metadynamics acceleration factor");
494 keys.addOutputComponent("maxbias", "CALC_MAX_BIAS", "the maximum of the metadynamics V(s, t)");
495 keys.addOutputComponent("transbias", "CALC_TRANSITION_BIAS", "the metadynamics transition bias V*(t)");
496 keys.addOutputComponent("pace","FREQUENCY_ADAPTIVE","the hill addition frequency when employing frequency adaptive metadynamics");
497 keys.use("ARG");
498 keys.add("compulsory","SIGMA","the widths of the Gaussian hills");
499 keys.add("compulsory","PACE","the frequency for hill addition");
500 keys.add("compulsory","FILE","HILLS","a file in which the list of added hills is stored");
501 keys.add("optional","HEIGHT","the heights of the Gaussian hills. Compulsory unless TAU and either BIASFACTOR or DAMPFACTOR are given");
502 keys.add("optional","FMT","specify format for HILLS files (useful for decrease the number of digits in regtests)");
503 keys.add("optional","BIASFACTOR","use well tempered metadynamics and use this bias factor. Please note you must also specify temp");
504 keys.add("optional","RECT","list of bias factors for all the replicas");
505 keys.add("optional","DAMPFACTOR","damp hills with exp(-max(V)/(\\f$k_B\\f$T*DAMPFACTOR)");
506 for (size_t i = 0; i < n_tempering_options_; i++) {
507 registerTemperingKeywords(tempering_names_[i][0], tempering_names_[i][1], keys);
508 }
509 keys.add("optional","TARGET","target to a predefined distribution");
510 keys.add("optional","TEMP","the system temperature - this is only needed if you are doing well-tempered metadynamics");
511 keys.add("optional","TAU","in well tempered metadynamics, sets height to (\\f$k_B \\Delta T\\f$*pace*timestep)/tau");
512 keys.add("optional","GRID_MIN","the lower bounds for the grid");
513 keys.add("optional","GRID_MAX","the upper bounds for the grid");
514 keys.add("optional","GRID_BIN","the number of bins for the grid");
515 keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
516 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 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 keys.addFlag("GRID_SPARSE",false,"use a sparse grid to store hills");
521 keys.addFlag("GRID_NOSPLINE",false,"don't use spline interpolation with grids");
522 keys.add("optional","GRID_WSTRIDE","write the grid to a file every N steps");
523 keys.add("optional","GRID_WFILE","the file on which to write the grid");
524 keys.add("optional","GRID_RFILE","a grid file from which the bias should be read at the initial step of the simulation");
525 keys.addFlag("STORE_GRIDS",false,"store all the grid files the calculation generates. They will be deleted if this keyword is not present");
526 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 keys.add("optional","WALKERS_ID", "walker id");
528 keys.add("optional","WALKERS_N", "number of walkers");
529 keys.add("optional","WALKERS_DIR", "shared directory with the hills files from all the walkers");
530 keys.add("optional","WALKERS_RSTRIDE","stride for reading hills files");
531 keys.add("optional","INTERVAL","one dimensional lower and upper limits, outside the limits the system will not feel the biasing force.");
532 keys.add("optional","SIGMA_MAX","the upper bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
533 keys.add("optional","SIGMA_MIN","the lower bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
534 keys.addFlag("WALKERS_MPI",false,"Switch on MPI version of multiple walkers - not compatible with WALKERS_* options other than WALKERS_DIR");
535 keys.addFlag("FLYING_GAUSSIAN",false,"Switch on flying Gaussian method, must be used with WALKERS_MPI");
536 keys.addFlag("ACCELERATION",false,"Set to TRUE if you want to compute the metadynamics acceleration factor.");
537 keys.add("optional","ACCELERATION_RFILE","a data file from which the acceleration should be read at the initial step of the simulation");
538 keys.addFlag("CALC_MAX_BIAS", false, "Set to TRUE if you want to compute the maximum of the metadynamics V(s, t)");
539 keys.addFlag("CALC_TRANSITION_BIAS", false, "Set to TRUE if you want to compute a metadynamics transition bias V*(t)");
540 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 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 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 keys.add("optional","FA_MAX_PACE","the maximum hill addition frequency allowed in frequency adaptive metadynamics. By default there is no maximum value.");
544 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 keys.use("RESTART");
546 keys.use("UPDATE_FROM");
547 keys.use("UPDATE_UNTIL");
548 }
549
550 const std::string MetaD::tempering_names_[1][2] = {{"TT", "transition tempered"}};
551
registerTemperingKeywords(const std::string & name_stem,const std::string & name,Keywords & keys)552 void MetaD::registerTemperingKeywords(const std::string &name_stem, const std::string &name, Keywords &keys) {
553 keys.add("optional", name_stem + "BIASFACTOR", "use " + name + " metadynamics with this bias factor. Please note you must also specify temp");
554 keys.add("optional", name_stem + "BIASTHRESHOLD", "use " + name + " metadynamics with this bias threshold. Please note you must also specify " + name_stem + "BIASFACTOR");
555 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 }
557
MetaD(const ActionOptions & ao)558 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 last_step_warn_grid(0)
591 {
592 // parse the flexible hills
593 string adaptiveoption;
594 adaptiveoption="NONE";
595 parse("ADAPTIVE",adaptiveoption);
596 if(adaptiveoption=="GEOM") {
597 log.printf(" Uses Geometry-based hills width: sigma must be in distance units and only one sigma is needed\n");
598 adaptive_=FlexibleBin::geometry;
599 } else if(adaptiveoption=="DIFF") {
600 log.printf(" Uses Diffusion-based hills width: sigma must be in time steps and only one sigma is needed\n");
601 adaptive_=FlexibleBin::diffusion;
602 } else if(adaptiveoption=="NONE") {
603 adaptive_=FlexibleBin::none;
604 } else {
605 error("I do not know this type of adaptive scheme");
606 }
607
608 parse("FMT",fmt);
609
610 // parse the sigma
611 parseVector("SIGMA",sigma0_);
612 if(adaptive_==FlexibleBin::none) {
613 // if you use normal sigma you need one sigma per argument
614 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 if(sigma0_.size()!=1) {
618 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 if(adaptive_==FlexibleBin::diffusion) {
622 if(int(sigma0_[0])-sigma0_[0]>1.e-9 || int(sigma0_[0])-sigma0_[0] <-1.e-9 || int(sigma0_[0])<1 ) {
623 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 parseVector("SIGMA_MIN",sigma0min_);
628 if(sigma0min_.size()>0 && sigma0min_.size()!=getNumberOfArguments()) {
629 error("the number of SIGMA_MIN values be the same of the number of the arguments");
630 } else if(sigma0min_.size()==0) {
631 sigma0min_.resize(getNumberOfArguments());
632 for(unsigned i=0; i<getNumberOfArguments(); i++) {sigma0min_[i]=-1.;}
633 }
634
635 parseVector("SIGMA_MAX",sigma0max_);
636 if(sigma0max_.size()>0 && sigma0max_.size()!=getNumberOfArguments()) {
637 error("the number of SIGMA_MAX values be the same of the number of the arguments");
638 } else if(sigma0max_.size()==0) {
639 sigma0max_.resize(getNumberOfArguments());
640 for(unsigned i=0; i<getNumberOfArguments(); i++) {sigma0max_[i]=-1.;}
641 }
642
643 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 parse("HEIGHT",height0_);
647 parse("PACE",stride_);
648 if(stride_<=0 ) error("frequency for hill addition is nonsensical");
649 current_stride = stride_;
650 string hillsfname="HILLS";
651 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 parseFlag("CALC_MAX_BIAS", calc_max_bias_);
657 parseFlag("CALC_TRANSITION_BIAS", calc_transition_bias_);
658
659 std::vector<double> rect_biasf_;
660 parseVector("RECT",rect_biasf_);
661 if(rect_biasf_.size()>0) {
662 int r=0;
663 if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
664 comm.Bcast(r,0);
665 biasf_=rect_biasf_[r];
666 log<<" You are using RECT\n";
667 } else {
668 parse("BIASFACTOR",biasf_);
669 }
670 if( biasf_<1.0 && biasf_!=-1.0) error("well tempered bias factor is nonsensical");
671 parse("DAMPFACTOR",dampfactor_);
672 double temp=0.0;
673 parse("TEMP",temp);
674 if(temp>0.0) kbt_=plumed.getAtoms().getKBoltzmann()*temp;
675 else kbt_=plumed.getAtoms().getKbT();
676 if(biasf_>=1.0) {
677 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 welltemp_=true;
679 }
680 if(dampfactor_>0.0) {
681 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 readTemperingSpecs(tt_specs_);
687 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 if (calc_transition_bias_) {
692 vector<double> tempcoords(getNumberOfArguments());
693 for (unsigned i = 0; ; i++) {
694 if (!parseNumberedVector("TRANSITIONWELL", i, tempcoords) ) break;
695 if (tempcoords.size() != getNumberOfArguments()) {
696 error("incorrect number of coordinates for transition tempering well");
697 }
698 transitionwells_.push_back(tempcoords);
699 }
700 }
701
702 parse("TARGET",targetfilename_);
703 if(targetfilename_.length()>0 && kbt_==0.0) error("with TARGET temperature must be specified");
704 double tau=0.0;
705 parse("TAU",tau);
706 if(tau==0.0) {
707 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 if(welltemp_) tau=(kbt_*(biasf_-1.0))/height0_*getTimeStep()*stride_;
710 else if(dampfactor_>0.0) tau=(kbt_*dampfactor_)/height0_*getTimeStep()*stride_;
711 } else {
712 if(height0_!=std::numeric_limits<double>::max()) error("At most one between HEIGHT and TAU should be specified");
713 if(welltemp_) {
714 if(biasf_!=1.0) height0_=(kbt_*(biasf_-1.0))/tau*getTimeStep()*stride_;
715 else height0_=kbt_/tau*getTimeStep()*stride_; // special case for gamma=1
716 }
717 else if(dampfactor_>0.0) height0_=(kbt_*dampfactor_)/tau*getTimeStep()*stride_;
718 else error("TAU only makes sense in well-tempered or damped metadynamics");
719 }
720
721 // Grid Stuff
722 vector<std::string> gmin(getNumberOfArguments());
723 parseVector("GRID_MIN",gmin);
724 if(gmin.size()!=getNumberOfArguments() && gmin.size()!=0) error("not enough values for GRID_MIN");
725 vector<std::string> gmax(getNumberOfArguments());
726 parseVector("GRID_MAX",gmax);
727 if(gmax.size()!=getNumberOfArguments() && gmax.size()!=0) error("not enough values for GRID_MAX");
728 vector<unsigned> gbin(getNumberOfArguments());
729 vector<double> gspacing;
730 parseVector("GRID_BIN",gbin);
731 if(gbin.size()!=getNumberOfArguments() && gbin.size()!=0) error("not enough values for GRID_BIN");
732 parseVector("GRID_SPACING",gspacing);
733 if(gspacing.size()!=getNumberOfArguments() && gspacing.size()!=0) error("not enough values for GRID_SPACING");
734 if(gmin.size()!=gmax.size()) error("GRID_MAX and GRID_MIN should be either present or absent");
735 if(gspacing.size()!=0 && gmin.size()==0) error("If GRID_SPACING is present also GRID_MIN should be present");
736 if(gbin.size()!=0 && gmin.size()==0) error("If GRID_SPACING is present also GRID_MIN should be present");
737 if(gmin.size()!=0) {
738 if(gbin.size()==0 && gspacing.size()==0) {
739 if(adaptive_==FlexibleBin::none) {
740 log<<" Binsize not specified, 1/5 of sigma will be be used\n";
741 plumed_assert(sigma0_.size()==getNumberOfArguments());
742 gspacing.resize(getNumberOfArguments());
743 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 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 log<<" Binsize not specified, 1/5 of sigma_min will be be used\n";
748 gspacing.resize(getNumberOfArguments());
749 for(unsigned i=0; i<gspacing.size(); i++) gspacing[i]=0.2*sigma0min_[i];
750 }
751 } else if(gspacing.size()!=0 && gbin.size()==0) {
752 log<<" The number of bins will be estimated from GRID_SPACING\n";
753 } else if(gspacing.size()!=0 && gbin.size()!=0) {
754 log<<" You specified both GRID_BIN and GRID_SPACING\n";
755 log<<" The more conservative (highest) number of bins will be used for each variable\n";
756 }
757 if(gbin.size()==0) gbin.assign(getNumberOfArguments(),1);
758 if(gspacing.size()!=0) for(unsigned i=0; i<getNumberOfArguments(); i++) {
759 double a,b;
760 Tools::convert(gmin[i],a);
761 Tools::convert(gmax[i],b);
762 unsigned n=((b-a)/gspacing[i])+1;
763 if(gbin[i]<n) gbin[i]=n;
764 }
765 }
766 bool sparsegrid=false;
767 parseFlag("GRID_SPARSE",sparsegrid);
768 bool nospline=false;
769 parseFlag("GRID_NOSPLINE",nospline);
770 bool spline=!nospline;
771 if(gbin.size()>0) {grid_=true;}
772 parse("GRID_WSTRIDE",wgridstride_);
773 string gridfilename_;
774 parse("GRID_WFILE",gridfilename_);
775 parseFlag("STORE_GRIDS",storeOldGrids_);
776 if(grid_ && gridfilename_.length()>0) {
777 if(wgridstride_==0 ) error("frequency with which to output grid not specified use GRID_WSTRIDE");
778 }
779
780 if(grid_ && wgridstride_>0) {
781 if(gridfilename_.length()==0) error("grid filename not specified use GRID_WFILE");
782 }
783 string gridreadfilename_;
784 parse("GRID_RFILE",gridreadfilename_);
785
786 if(!grid_&&gridfilename_.length()> 0) error("To write a grid you need first to define it!");
787 if(!grid_&&gridreadfilename_.length()>0) error("To read a grid you need first to define it!");
788
789 // Reweighting factor rct
790 parseFlag("CALC_RCT",calc_rct_);
791 if (calc_rct_)
792 plumed_massert(grid_,"CALC_RCT is supported only if bias is on a grid");
793 parse("RCT_USTRIDE",rct_ustride_);
794
795 if(dampfactor_>0.0) {
796 if(!grid_) error("With DAMPFACTOR you should use grids");
797 }
798
799 // Multiple walkers
800 parse("WALKERS_N",mw_n_);
801 parse("WALKERS_ID",mw_id_);
802 if(mw_n_<=mw_id_) error("walker ID should be a numerical value less than the total number of walkers");
803 parse("WALKERS_DIR",mw_dir_);
804 parse("WALKERS_RSTRIDE",mw_rstride_);
805
806 // MPI version
807 parseFlag("WALKERS_MPI",walkers_mpi);
808
809 // Flying Gaussian
810 parseFlag("FLYING_GAUSSIAN", flying);
811
812 // Inteval keyword
813 vector<double> tmpI(2);
814 parseVector("INTERVAL",tmpI);
815 if(tmpI.size()!=2&&tmpI.size()!=0) error("both a lower and an upper limits must be provided with INTERVAL");
816 else if(tmpI.size()==2) {
817 lowI_=tmpI.at(0);
818 uppI_=tmpI.at(1);
819 if(getNumberOfArguments()!=1) error("INTERVAL limits correction works only for monodimensional metadynamics!");
820 if(uppI_<lowI_) error("The Upper limit must be greater than the Lower limit!");
821 if(getPntrToArgument(0)->isPeriodic()) error("INTERVAL cannot be used with periodic variables!");
822 doInt_=true;
823 }
824
825 acceleration=false;
826 parseFlag("ACCELERATION",acceleration);
827 // Check for a restart acceleration if acceleration is active.
828 string acc_rfilename;
829 if (acceleration) {
830 parse("ACCELERATION_RFILE", acc_rfilename);
831 }
832
833 freq_adaptive_=false;
834 parseFlag("FREQUENCY_ADAPTIVE",freq_adaptive_);
835 //
836 fa_update_frequency_=0;
837 parse("FA_UPDATE_FREQUENCY",fa_update_frequency_);
838 if(fa_update_frequency_!=0 && !freq_adaptive_) {
839 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 if(fa_update_frequency_==0 && freq_adaptive_) {
842 fa_update_frequency_=stride_;
843 }
844 //
845 fa_max_stride_=0;
846 parse("FA_MAX_PACE",fa_max_stride_);
847 if(fa_max_stride_!=0 && !freq_adaptive_) {
848 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 fa_min_acceleration_=1.0;
852 parse("FA_MIN_ACCELERATION",fa_min_acceleration_);
853 if(fa_min_acceleration_!=1.0 && !freq_adaptive_) {
854 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 checkRead();
858
859 log.printf(" Gaussian width ");
860 if (adaptive_==FlexibleBin::diffusion)log.printf(" (Note: The units of sigma are in timesteps) ");
861 if (adaptive_==FlexibleBin::geometry)log.printf(" (Note: The units of sigma are in dist units) ");
862 for(unsigned i=0; i<sigma0_.size(); ++i) log.printf(" %f",sigma0_[i]);
863 log.printf(" Gaussian height %f\n",height0_);
864 log.printf(" Gaussian deposition pace %d\n",stride_);
865 log.printf(" Gaussian file %s\n",hillsfname.c_str());
866 if(welltemp_) {
867 log.printf(" Well-Tempered Bias Factor %f\n",biasf_);
868 log.printf(" Hills relaxation time (tau) %f\n",tau);
869 log.printf(" KbT %f\n",kbt_);
870 }
871 // Transition tempered metadynamics options
872 if (tt_specs_.is_active) {
873 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 if (!calc_transition_bias_) {
877 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 if (welltemp_ || dampfactor_ > 0.0 || tt_specs_.is_active) {
891 // Determine the number of active temperings.
892 int n_active = 0;
893 if (welltemp_) n_active++;
894 if (dampfactor_ > 0.0) n_active++;
895 if (tt_specs_.is_active) n_active++;
896 // Find the greatest alpha.
897 double greatest_alpha = 0.0;
898 if (welltemp_) greatest_alpha = max(greatest_alpha, 1.0);
899 if (dampfactor_ > 0.0) greatest_alpha = max(greatest_alpha, 1.0);
900 if (tt_specs_.is_active) greatest_alpha = max(greatest_alpha, tt_specs_.alpha);
901 // Find the least alpha.
902 double least_alpha = 1.0;
903 if (welltemp_) least_alpha = min(least_alpha, 1.0);
904 if (dampfactor_ > 0.0) least_alpha = min(least_alpha, 1.0);
905 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 if (welltemp_ && 1.0 == greatest_alpha && biasf_ != 1.0) total_governing_deltaT_inv += 1.0 / (biasf_ - 1.0);
910 if (dampfactor_ > 0.0 && 1.0 == greatest_alpha) total_governing_deltaT_inv += 1.0 / (dampfactor_);
911 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 if (n_active == 1 && total_governing_deltaT_inv < 0.0) {
915 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 } else if (total_governing_deltaT_inv < 0.0 && greatest_alpha == least_alpha) {
919 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 } else if (total_governing_deltaT_inv < 0.0 && greatest_alpha != least_alpha) {
923 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 if(doInt_) log.printf(" Upper and Lower limits boundaries for the bias are activated at %f - %f\n", lowI_, uppI_);
928 if(grid_) {
929 log.printf(" Grid min");
930 for(unsigned i=0; i<gmin.size(); ++i) log.printf(" %s",gmin[i].c_str() );
931 log.printf("\n");
932 log.printf(" Grid max");
933 for(unsigned i=0; i<gmax.size(); ++i) log.printf(" %s",gmax[i].c_str() );
934 log.printf("\n");
935 log.printf(" Grid bin");
936 for(unsigned i=0; i<gbin.size(); ++i) log.printf(" %u",gbin[i]);
937 log.printf("\n");
938 if(spline) {log.printf(" Grid uses spline interpolation\n");}
939 if(sparsegrid) {log.printf(" Grid uses sparse grid\n");}
940 if(wgridstride_>0) {log.printf(" Grid is written on file %s with stride %d\n",gridfilename_.c_str(),wgridstride_);}
941 }
942
943 if(mw_n_>1) {
944 if(walkers_mpi) error("MPI version of multiple walkers is not compatible with filesystem version of multiple walkers");
945 log.printf(" %d multiple walkers active\n",mw_n_);
946 log.printf(" walker id %d\n",mw_id_);
947 log.printf(" reading stride %d\n",mw_rstride_);
948 if(mw_dir_!="")log.printf(" directory with hills files %s\n",mw_dir_.c_str());
949 } else {
950 if(walkers_mpi) {
951 log.printf(" Multiple walkers active using MPI communnication\n");
952 if(mw_dir_!="")log.printf(" directory with hills files %s\n",mw_dir_.c_str());
953 if(comm.Get_rank()==0) {
954 // Only root of group can communicate with other walkers
955 mpi_nw_=multi_sim_comm.Get_size();
956 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 comm.Bcast(mpi_nw_,0);
961 comm.Bcast(mpi_mw_,0);
962 }
963 }
964
965 if(flying) {
966 if(!walkers_mpi) error("Flying Gaussian method must be used with MPI version of multiple walkers");
967 log.printf(" Flying Gaussian method with %d walkers active\n",mpi_nw_);
968 }
969
970 if(calc_rct_) {
971 addComponent("rbias"); componentIsNotPeriodic("rbias");
972 addComponent("rct"); componentIsNotPeriodic("rct");
973 log.printf(" The c(t) reweighting factor will be calculated every %u hills\n",rct_ustride_);
974 getPntrToComponent("rct")->set(reweight_factor_);
975 }
976 addComponent("work"); componentIsNotPeriodic("work");
977
978 if(acceleration) {
979 if (kbt_ == 0.0) {
980 error("The calculation of the acceleration works only if simulation temperature has been defined");
981 }
982 log.printf(" calculation on the fly of the acceleration factor\n");
983 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 if (acc_rfilename.length() == 0) {
987 getPntrToComponent("acc")->set(1.0);
988 if(getRestart()) {
989 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 IFile acc_rfile;
999 acc_rfile.link(*this);
1000 if(acc_rfile.FileExist(acc_rfilename)) {
1001 acc_rfile.open(acc_rfilename);
1002 } else {
1003 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 std::string acclabel = getLabel() + ".acc";
1009 acc_rfile.allowIgnoredFields();
1010 while(acc_rfile.scanField("time", acc_rtime)) {
1011 acc_rfile.scanField(acclabel, acc_rmean);
1012 acc_rfile.scanField();
1013 }
1014 acc_restart_mean_ = acc_rmean;
1015 // Set component based on the read values.
1016 getPntrToComponent("acc")->set(acc_rmean);
1017 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 if (calc_max_bias_) {
1021 if (!grid_) error("Calculating the maximum bias on the fly works only with a grid");
1022 log.printf(" calculation on the fly of the maximum bias max(V(s,t)) \n");
1023 addComponent("maxbias");
1024 componentIsNotPeriodic("maxbias");
1025 }
1026 if (calc_transition_bias_) {
1027 if (!grid_) error("Calculating the transition bias on the fly works only with a grid");
1028 log.printf(" calculation on the fly of the transition bias V*(t)\n");
1029 addComponent("transbias");
1030 componentIsNotPeriodic("transbias");
1031 log<<" Number of transition wells "<<transitionwells_.size()<<"\n";
1032 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 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 for (unsigned i = 0; i < transitionwells_.size(); i++) {
1037 // Log the coordinate.
1038 log.printf(" Transition well %d at coordinate ", i);
1039 for (unsigned j = 0; j < getNumberOfArguments(); j++) log.printf("%f ", transitionwells_[i][j]);
1040 log.printf("\n");
1041 // Check that the coordinate is in the grid.
1042 for (unsigned j = 0; j < getNumberOfArguments(); j++) {
1043 double max, min;
1044 Tools::convert(gmin[j], min);
1045 Tools::convert(gmax[j], max);
1046 if (transitionwells_[i][j] < min || transitionwells_[i][j] > max) error(" transition well is not in grid");
1047 }
1048 }
1049 }
1050
1051 if(freq_adaptive_) {
1052 if(!acceleration) {
1053 plumed_merror("Frequency adaptive metadynamics only works if the calculation of the acceleration factor is enabled with the ACCELERATION keyword\n");
1054 }
1055 if(walkers_mpi) {
1056 plumed_merror("Combining frequency adaptive metadynamics with MPI multiple walkers is not allowed");
1057 }
1058
1059 log.printf(" Frequency adaptive metadynamics enabled\n");
1060 if(getRestart() && acc_rfilename.length() == 0) {
1061 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 log.printf(" The frequency for hill addition will change dynamically based on the metadynamics acceleration factor\n");
1064 log.printf(" The hill addition frequency will be updated every %d steps\n",fa_update_frequency_);
1065 if(fa_min_acceleration_>1.0) {
1066 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 if(fa_max_stride_!=0) {
1069 log.printf(" The hill addition frequency will not become larger than %d steps\n",fa_max_stride_);
1070 }
1071 addComponent("pace"); componentIsNotPeriodic("pace");
1072 updateFrequencyAdaptiveStride();
1073 }
1074
1075 // for performance
1076 dp_.reset( new double[getNumberOfArguments()] );
1077
1078 // initializing and checking grid
1079 if(grid_) {
1080 // check for mesh and sigma size
1081 for(unsigned i=0; i<getNumberOfArguments(); i++) {
1082 double a,b;
1083 Tools::convert(gmin[i],a);
1084 Tools::convert(gmax[i],b);
1085 double mesh=(b-a)/((double)gbin[i]);
1086 if(adaptive_==FlexibleBin::none) {
1087 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 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 std::string funcl=getLabel() + ".bias";
1093 if(!sparsegrid) {BiasGrid_.reset(new Grid(funcl,getArguments(),gmin,gmax,gbin,spline,true));}
1094 else {BiasGrid_.reset(new SparseGrid(funcl,getArguments(),gmin,gmax,gbin,spline,true));}
1095 std::vector<std::string> actualmin=BiasGrid_->getMin();
1096 std::vector<std::string> actualmax=BiasGrid_->getMax();
1097 for(unsigned i=0; i<getNumberOfArguments(); i++) {
1098 std::string is;
1099 Tools::convert(i,is);
1100 if(gmin[i]!=actualmin[i]) error("GRID_MIN["+is+"] must be adjusted to "+actualmin[i]+" to fit periodicity");
1101 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 if(gridreadfilename_.length()>0) {
1108 // read the grid in input, find the keys
1109 IFile gridfile;
1110 gridfile.link(*this);
1111 if(gridfile.FileExist(gridreadfilename_)) {
1112 gridfile.open(gridreadfilename_);
1113 } else {
1114 error("The GRID file you want to read: " + gridreadfilename_ + ", cannot be found!");
1115 }
1116 std::string funcl=getLabel() + ".bias";
1117 BiasGrid_=GridBase::create(funcl, getArguments(), gridfile, gmin, gmax, gbin, sparsegrid, spline, true);
1118 if(BiasGrid_->getDimension()!=getNumberOfArguments()) error("mismatch between dimensionality of input grid and number of arguments");
1119 for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1120 if( getPntrToArgument(i)->isPeriodic()!=BiasGrid_->getIsPeriodic()[i] ) error("periodicity mismatch between arguments and input bias");
1121 double a, b;
1122 Tools::convert(gmin[i],a);
1123 Tools::convert(gmax[i],b);
1124 double mesh=(b-a)/((double)gbin[i]);
1125 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 log.printf(" Restarting from %s:",gridreadfilename_.c_str());
1128 if(getRestart()) restartedFromGrid=true;
1129 }
1130
1131 // initializing and checking grid
1132 if(grid_&&!(gridreadfilename_.length()>0)) {
1133 // check for adaptive and sigma_min
1134 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 for(unsigned i=0; i<getNumberOfArguments(); i++) {
1137 double a,b;
1138 Tools::convert(gmin[i],a);
1139 Tools::convert(gmax[i],b);
1140 double mesh=(b-a)/((double)gbin[i]);
1141 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 std::string funcl=getLabel() + ".bias";
1144 if(!sparsegrid) {BiasGrid_.reset(new Grid(funcl,getArguments(),gmin,gmax,gbin,spline,true));}
1145 else {BiasGrid_.reset(new SparseGrid(funcl,getArguments(),gmin,gmax,gbin,spline,true));}
1146 std::vector<std::string> actualmin=BiasGrid_->getMin();
1147 std::vector<std::string> actualmax=BiasGrid_->getMax();
1148 for(unsigned i=0; i<getNumberOfArguments(); i++) {
1149 if(gmin[i]!=actualmin[i]) log<<" WARNING: GRID_MIN["<<i<<"] has been adjusted to "<<actualmin[i]<<" to fit periodicity\n";
1150 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 for(int i=0; i<mw_n_; ++i) {
1157 string fname;
1158 if(mw_dir_!="") {
1159 if(mw_n_>1) {
1160 stringstream out; out << i;
1161 fname = mw_dir_+"/"+hillsfname+"."+out.str();
1162 } else if(walkers_mpi) {
1163 fname = mw_dir_+"/"+hillsfname;
1164 } else {
1165 fname = hillsfname;
1166 }
1167 } else {
1168 if(mw_n_>1) {
1169 stringstream out; out << i;
1170 fname = hillsfname+"."+out.str();
1171 } else {
1172 fname = hillsfname;
1173 }
1174 }
1175 ifiles.emplace_back(new IFile());
1176 // this is just a shortcut pointer to the last element:
1177 IFile *ifile = ifiles.back().get();
1178 ifilesnames.push_back(fname);
1179 ifile->link(*this);
1180 if(ifile->FileExist(fname)) {
1181 ifile->open(fname);
1182 if(getRestart()&&!restartedFromGrid) {
1183 log.printf(" Restarting from %s:",ifilesnames[i].c_str());
1184 readGaussians(ifiles[i].get());
1185 }
1186 ifiles[i]->reset(false);
1187 // close only the walker own hills file for later writing
1188 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 if(getRestart()) log<<" WARNING: restart file "<<fname<<" not found\n";
1192 }
1193 }
1194
1195 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 if(comm.Get_rank()==0 && walkers_mpi) multi_sim_comm.Barrier();
1204 if(targetfilename_.length()>0) {
1205 IFile gridfile; gridfile.open(targetfilename_);
1206 std::string funcl=getLabel() + ".target";
1207 TargetGrid_=GridBase::create(funcl,getArguments(),gridfile,false,false,true);
1208 if(TargetGrid_->getDimension()!=getNumberOfArguments()) error("mismatch between dimensionality of input grid and number of arguments");
1209 for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1210 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 if(getRestart() && calc_rct_) computeReweightingFactor();
1216 // Calculate all special bias quantities desired if restarting with nonzero bias.
1217 if(getRestart() && calc_max_bias_) {
1218 max_bias_ = BiasGrid_->getMaxValue();
1219 getPntrToComponent("maxbias")->set(max_bias_);
1220 }
1221 if(getRestart() && calc_transition_bias_) {
1222 transition_bias_ = getTransitionBarrierBias();
1223 getPntrToComponent("transbias")->set(transition_bias_);
1224 }
1225
1226 // open grid file for writing
1227 if(wgridstride_>0) {
1228 gridfile_.link(*this);
1229 if(walkers_mpi) {
1230 int r=0;
1231 if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
1232 comm.Bcast(r,0);
1233 if(r>0) gridfilename_="/dev/null";
1234 gridfile_.enforceSuffix("");
1235 }
1236 if(mw_n_>1) gridfile_.enforceSuffix("");
1237 gridfile_.open(gridfilename_);
1238 }
1239
1240 // open hills file for writing
1241 hillsOfile_.link(*this);
1242 if(walkers_mpi) {
1243 int r=0;
1244 if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
1245 comm.Bcast(r,0);
1246 if(r>0) ifilesnames[mw_id_]="/dev/null";
1247 hillsOfile_.enforceSuffix("");
1248 }
1249 if(mw_n_>1) hillsOfile_.enforceSuffix("");
1250 hillsOfile_.open(ifilesnames[mw_id_]);
1251 if(fmt.length()>0) hillsOfile_.fmtField(fmt);
1252 hillsOfile_.addConstantField("multivariate");
1253 hillsOfile_.addConstantField("kerneltype");
1254 if(doInt_) {
1255 hillsOfile_.addConstantField("lower_int").printField("lower_int",lowI_);
1256 hillsOfile_.addConstantField("upper_int").printField("upper_int",uppI_);
1257 }
1258 hillsOfile_.setHeavyFlush();
1259 // output periodicities of variables
1260 for(unsigned i=0; i<getNumberOfArguments(); ++i) hillsOfile_.setupPrintValue( getPntrToArgument(i) );
1261
1262 bool concurrent=false;
1263 const ActionSet&actionSet(plumed.getActionSet());
1264 for(const auto & p : actionSet) if(dynamic_cast<MetaD*>(p.get())) { concurrent=true; break; }
1265 if(concurrent) log<<" You are using concurrent metadynamics\n";
1266 if(rect_biasf_.size()>0) {
1267 if(walkers_mpi) {
1268 log<<" You are using RECT in its 'altruistic' implementation\n";
1269 }{
1270 log<<" You are using RECT\n";
1271 }
1272 }
1273
1274 log<<" Bibliography "<<plumed.cite("Laio and Parrinello, PNAS 99, 12562 (2002)");
1275 if(welltemp_) log<<plumed.cite(
1276 "Barducci, Bussi, and Parrinello, Phys. Rev. Lett. 100, 020603 (2008)");
1277 if(tt_specs_.is_active) {
1278 log << plumed.cite("Dama, Rotskoff, Parrinello, and Voth, J. Chem. Theory Comput. 10, 3626 (2014)");
1279 log << plumed.cite("Dama, Parrinello, and Voth, Phys. Rev. Lett. 112, 240602 (2014)");
1280 }
1281 if(mw_n_>1||walkers_mpi) log<<plumed.cite(
1282 "Raiteri, Laio, Gervasio, Micheletti, and Parrinello, J. Phys. Chem. B 110, 3533 (2006)");
1283 if(adaptive_!=FlexibleBin::none) log<<plumed.cite(
1284 "Branduardi, Bussi, and Parrinello, J. Chem. Theory Comput. 8, 2247 (2012)");
1285 if(doInt_) log<<plumed.cite(
1286 "Baftizadeh, Cossio, Pietrucci, and Laio, Curr. Phys. Chem. 2, 79 (2012)");
1287 if(acceleration) log<<plumed.cite(
1288 "Pratyush and Parrinello, Phys. Rev. Lett. 111, 230602 (2013)");
1289 if(calc_rct_) log<<plumed.cite(
1290 "Pratyush and Parrinello, J. Phys. Chem. B, 119, 736 (2015)");
1291 if(concurrent || rect_biasf_.size()>0) log<<plumed.cite(
1292 "Gil-Ley and Bussi, J. Chem. Theory Comput. 11, 1077 (2015)");
1293 if(rect_biasf_.size()>0 && walkers_mpi) log<<plumed.cite(
1294 "Hosek, Toulcova, Bortolato, and Spiwok, J. Phys. Chem. B 120, 2209 (2016)");
1295 if(targetfilename_.length()>0) {
1296 log<<plumed.cite("White, Dama, and Voth, J. Chem. Theory Comput. 11, 2451 (2015)");
1297 log<<plumed.cite("Marinelli and Faraldo-Gómez, Biophys. J. 108, 2779 (2015)");
1298 log<<plumed.cite("Gil-Ley, Bottaro, and Bussi, J. Chem. Theory Comput. 12, 2790 (2016)");
1299 }
1300 if(freq_adaptive_) {
1301 log<<plumed.cite("Wang, Valsson, Tiwary, Parrinello, and Lindorff-Larsen, J. Chem. Phys. 149, 072309 (2018)");
1302 }
1303 log<<"\n";
1304 }
1305
readTemperingSpecs(TemperingSpecs & t_specs)1306 void MetaD::readTemperingSpecs(TemperingSpecs &t_specs) {
1307 // Set global tempering parameters.
1308 parse(t_specs.name_stem + "BIASFACTOR", t_specs.biasf);
1309 if (t_specs.biasf != -1.0) {
1310 if (kbt_ == 0.0) {
1311 error("Unless the MD engine passes the temperature to plumed, with tempered metad you must specify it using TEMP");
1312 }
1313 if (t_specs.biasf == 1.0) {
1314 error("A bias factor of 1 corresponds to zero delta T and zero hill size, so it is not allowed.");
1315 }
1316 t_specs.is_active = true;
1317 parse(t_specs.name_stem + "BIASTHRESHOLD", t_specs.threshold);
1318 if (t_specs.threshold < 0.0) {
1319 error(t_specs.name + " bias threshold is nonsensical");
1320 }
1321 parse(t_specs.name_stem + "ALPHA", t_specs.alpha);
1322 if (t_specs.alpha <= 0.0 || t_specs.alpha > 1.0) {
1323 error(t_specs.name + " decay shape parameter alpha is nonsensical");
1324 }
1325 }
1326 }
1327
logTemperingSpecs(const TemperingSpecs & t_specs)1328 void MetaD::logTemperingSpecs(const TemperingSpecs &t_specs) {
1329 log.printf(" %s bias factor %f\n", t_specs.name.c_str(), t_specs.biasf);
1330 log.printf(" KbT %f\n", kbt_);
1331 if (t_specs.threshold != 0.0) log.printf(" %s bias threshold %f\n", t_specs.name.c_str(), t_specs.threshold);
1332 if (t_specs.alpha != 1.0) log.printf(" %s decay shape parameter alpha %f\n", t_specs.name.c_str(), t_specs.alpha);
1333 }
1334
readGaussians(IFile * ifile)1335 void MetaD::readGaussians(IFile *ifile)
1336 {
1337 unsigned ncv=getNumberOfArguments();
1338 vector<double> center(ncv);
1339 vector<double> sigma(ncv);
1340 double height;
1341 int nhills=0;
1342 bool multivariate=false;
1343
1344 std::vector<Value> tmpvalues;
1345 for(unsigned j=0; j<getNumberOfArguments(); ++j) tmpvalues.push_back( Value( this, getPntrToArgument(j)->getName(), false ) );
1346
1347 while(scanOneHill(ifile,tmpvalues,center,sigma,height,multivariate)) {
1348 ;
1349 nhills++;
1350 // note that for gamma=1 we store directly -F
1351 if(welltemp_ && biasf_>1.0) {height*=(biasf_-1.0)/biasf_;}
1352 addGaussian(Gaussian(center,sigma,height,multivariate));
1353 }
1354 log.printf(" %d Gaussians read\n",nhills);
1355 }
1356
writeGaussian(const Gaussian & hill,OFile & file)1357 void MetaD::writeGaussian(const Gaussian& hill, OFile&file)
1358 {
1359 unsigned ncv=getNumberOfArguments();
1360 file.printField("time",getTimeStep()*getStep());
1361 for(unsigned i=0; i<ncv; ++i) {
1362 file.printField(getPntrToArgument(i),hill.center[i]);
1363 }
1364 hillsOfile_.printField("kerneltype","gaussian");
1365 if(hill.multivariate) {
1366 hillsOfile_.printField("multivariate","true");
1367 Matrix<double> mymatrix(ncv,ncv);
1368 unsigned k=0;
1369 for(unsigned i=0; i<ncv; i++) {
1370 for(unsigned j=i; j<ncv; j++) {
1371 // recompose the full inverse matrix
1372 mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k];
1373 k++;
1374 }
1375 }
1376 // invert it
1377 Matrix<double> invmatrix(ncv,ncv);
1378 Invert(mymatrix,invmatrix);
1379 // enforce symmetry
1380 for(unsigned i=0; i<ncv; i++) {
1381 for(unsigned j=i; j<ncv; j++) {
1382 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 cholesky(invmatrix,lower);
1389 // loop in band form
1390 for(unsigned i=0; i<ncv; i++) {
1391 for(unsigned j=0; j<ncv-i; j++) {
1392 file.printField("sigma_"+getPntrToArgument(j+i)->getName()+"_"+getPntrToArgument(j)->getName(),lower(j+i,j));
1393 }
1394 }
1395 } else {
1396 hillsOfile_.printField("multivariate","false");
1397 for(unsigned i=0; i<ncv; ++i)
1398 file.printField("sigma_"+getPntrToArgument(i)->getName(),hill.sigma[i]);
1399 }
1400 double height=hill.height;
1401 // note that for gamma=1 we store directly -F
1402 if(welltemp_ && biasf_>1.0) height*=biasf_/(biasf_-1.0);
1403 file.printField("height",height).printField("biasf",biasf_);
1404 if(mw_n_>1) file.printField("clock",int(std::time(0)));
1405 file.printField();
1406 }
1407
addGaussian(const Gaussian & hill)1408 void MetaD::addGaussian(const Gaussian& hill)
1409 {
1410 if(!grid_) hills_.push_back(hill);
1411 else {
1412 unsigned ncv=getNumberOfArguments();
1413 vector<unsigned> nneighb=getGaussianSupport(hill);
1414 vector<Grid::index_t> neighbors=BiasGrid_->getNeighbors(hill.center,nneighb);
1415 vector<double> der(ncv);
1416 vector<double> xx(ncv);
1417 if(comm.Get_size()==1) {
1418 for(unsigned i=0; i<neighbors.size(); ++i) {
1419 Grid::index_t ineigh=neighbors[i];
1420 for(unsigned j=0; j<ncv; ++j) der[j]=0.0;
1421 BiasGrid_->getPoint(ineigh,xx);
1422 double bias=evaluateGaussian(xx,hill,&der[0]);
1423 BiasGrid_->addValueAndDerivatives(ineigh,bias,der);
1424 }
1425 } else {
1426 unsigned stride=comm.Get_size();
1427 unsigned rank=comm.Get_rank();
1428 vector<double> allder(ncv*neighbors.size(),0.0);
1429 vector<double> allbias(neighbors.size(),0.0);
1430 for(unsigned i=rank; i<neighbors.size(); i+=stride) {
1431 Grid::index_t ineigh=neighbors[i];
1432 BiasGrid_->getPoint(ineigh,xx);
1433 allbias[i]=evaluateGaussian(xx,hill,&allder[ncv*i]);
1434 }
1435 comm.Sum(allbias);
1436 comm.Sum(allder);
1437 for(unsigned i=0; i<neighbors.size(); ++i) {
1438 Grid::index_t ineigh=neighbors[i];
1439 for(unsigned j=0; j<ncv; ++j) {der[j]=allder[ncv*i+j];}
1440 BiasGrid_->addValueAndDerivatives(ineigh,allbias[i],der);
1441 }
1442 }
1443 }
1444 }
1445
getGaussianSupport(const Gaussian & hill)1446 vector<unsigned> MetaD::getGaussianSupport(const Gaussian& hill)
1447 {
1448 vector<unsigned> nneigh;
1449 vector<double> cutoff;
1450 unsigned ncv=getNumberOfArguments();
1451
1452 // traditional or flexible hill?
1453 if(hill.multivariate) {
1454 unsigned k=0;
1455 Matrix<double> mymatrix(ncv,ncv);
1456 for(unsigned i=0; i<ncv; i++) {
1457 for(unsigned j=i; j<ncv; j++) {
1458 // recompose the full inverse matrix
1459 mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k];
1460 k++;
1461 }
1462 }
1463 // Reinvert so to have the ellipses
1464 Matrix<double> myinv(ncv,ncv);
1465 Invert(mymatrix,myinv);
1466 Matrix<double> myautovec(ncv,ncv);
1467 vector<double> myautoval(ncv); //should I take this or their square root?
1468 diagMat(myinv,myautoval,myautovec);
1469 double maxautoval=0.;
1470 unsigned ind_maxautoval; ind_maxautoval=ncv;
1471 for(unsigned i=0; i<ncv; i++) {
1472 if(myautoval[i]>maxautoval) {maxautoval=myautoval[i]; ind_maxautoval=i;}
1473 }
1474 for(unsigned i=0; i<ncv; i++) {
1475 cutoff.push_back(sqrt(2.0*DP2CUTOFF)*abs(sqrt(maxautoval)*myautovec(i,ind_maxautoval)));
1476 }
1477 } else {
1478 for(unsigned i=0; i<ncv; ++i) {
1479 cutoff.push_back(sqrt(2.0*DP2CUTOFF)*hill.sigma[i]);
1480 }
1481 }
1482
1483 if(doInt_) {
1484 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 return BiasGrid_->getNbin();
1487 } else {
1488 nneigh.push_back( static_cast<unsigned>(ceil(cutoff[0]/BiasGrid_->getDx()[0])) );
1489 return nneigh;
1490 }
1491 } else {
1492 for(unsigned i=0; i<ncv; i++) {
1493 nneigh.push_back( static_cast<unsigned>(ceil(cutoff[i]/BiasGrid_->getDx()[i])) );
1494 }
1495 }
1496
1497 return nneigh;
1498 }
1499
getBiasAndDerivatives(const vector<double> & cv,double * der)1500 double MetaD::getBiasAndDerivatives(const vector<double>& cv, double* der)
1501 {
1502 double bias=0.0;
1503 if(!grid_) {
1504 if(hills_.size()>10000 && (getStep()-last_step_warn_grid)>10000) {
1505 std::string msg;
1506 Tools::convert(hills_.size(),msg);
1507 msg="You have accumulated "+msg+" hills, you should enable GRIDs to avoid serious performance hits";
1508 warning(msg);
1509 last_step_warn_grid=getStep();
1510 }
1511 unsigned stride=comm.Get_size();
1512 unsigned rank=comm.Get_rank();
1513 for(unsigned i=rank; i<hills_.size(); i+=stride) {
1514 bias+=evaluateGaussian(cv,hills_[i],der);
1515 }
1516 comm.Sum(bias);
1517 if(der) comm.Sum(der,getNumberOfArguments());
1518 } else {
1519 if(der) {
1520 vector<double> vder(getNumberOfArguments());
1521 bias=BiasGrid_->getValueAndDerivatives(cv,vder);
1522 for(unsigned i=0; i<getNumberOfArguments(); ++i) {der[i]=vder[i];}
1523 } else {
1524 bias = BiasGrid_->getValue(cv);
1525 }
1526 }
1527
1528 return bias;
1529 }
1530
getGaussianNormalization(const Gaussian & hill)1531 double MetaD::getGaussianNormalization( const Gaussian& hill )
1532 {
1533 double norm=1;
1534 unsigned ncv=hill.center.size();
1535
1536 if(hill.multivariate) {
1537 // recompose the full sigma from the upper diag cholesky
1538 unsigned k=0;
1539 Matrix<double> mymatrix(ncv,ncv);
1540 for(unsigned i=0; i<ncv; i++) {
1541 for(unsigned j=i; j<ncv; j++) {
1542 mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k]; // recompose the full inverse matrix
1543 k++;
1544 }
1545 double ldet; logdet( mymatrix, ldet );
1546 norm = exp( ldet ); // Not sure here if mymatrix is sigma or inverse
1547 }
1548 } else {
1549 for(unsigned i=0; i<hill.sigma.size(); ++i) norm*=hill.sigma[i];
1550 }
1551
1552 return norm*pow(2*pi,static_cast<double>(ncv)/2.0);
1553 }
1554
evaluateGaussian(const vector<double> & cv,const Gaussian & hill,double * der)1555 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 if(cv.size()>0) pcv=&cv[0];
1565 if(doInt_) {
1566 plumed_assert(cv.size()==1);
1567 tmpcv[0]=cv[0];
1568 if(cv[0]<lowI_) tmpcv[0]=lowI_;
1569 if(cv[0]>uppI_) tmpcv[0]=uppI_;
1570 pcv=&(tmpcv[0]);
1571 }
1572 if(hill.multivariate) {
1573 unsigned k=0;
1574 unsigned ncv=cv.size();
1575 // recompose the full sigma from the upper diag cholesky
1576 Matrix<double> mymatrix(ncv,ncv);
1577 for(unsigned i=0; i<ncv; i++) {
1578 for(unsigned j=i; j<ncv; j++) {
1579 mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k]; // recompose the full inverse matrix
1580 k++;
1581 }
1582 }
1583 for(unsigned i=0; i<cv.size(); ++i) {
1584 double dp_i=difference(i,hill.center[i],pcv[i]);
1585 dp_[i]=dp_i;
1586 for(unsigned j=i; j<cv.size(); ++j) {
1587 if(i==j) {
1588 dp2+=dp_i*dp_i*mymatrix(i,j)*0.5;
1589 } else {
1590 double dp_j=difference(j,hill.center[j],pcv[j]);
1591 dp2+=dp_i*dp_j*mymatrix(i,j);
1592 }
1593 }
1594 }
1595 if(dp2<DP2CUTOFF) {
1596 bias=hill.height*exp(-dp2);
1597 if(der) {
1598 for(unsigned i=0; i<cv.size(); ++i) {
1599 double tmp=0.0;
1600 for(unsigned j=0; j<cv.size(); ++j) {
1601 tmp += dp_[j]*mymatrix(i,j)*bias;
1602 }
1603 der[i]-=tmp;
1604 }
1605 }
1606 }
1607 } else {
1608 for(unsigned i=0; i<cv.size(); ++i) {
1609 double dp=difference(i,hill.center[i],pcv[i])*hill.invsigma[i];
1610 dp2+=dp*dp;
1611 dp_[i]=dp;
1612 }
1613 dp2*=0.5;
1614 if(dp2<DP2CUTOFF) {
1615 bias=hill.height*exp(-dp2);
1616 if(der) {
1617 for(unsigned i=0; i<cv.size(); ++i) {der[i]+=-bias*dp_[i]*hill.invsigma[i];}
1618 }
1619 }
1620 }
1621
1622 if(doInt_ && der) {
1623 if(cv[0]<lowI_ || cv[0]>uppI_) for(unsigned i=0; i<cv.size(); ++i) der[i]=0;
1624 }
1625
1626 return bias;
1627 }
1628
getHeight(const vector<double> & cv)1629 double MetaD::getHeight(const vector<double>& cv)
1630 {
1631 double height=height0_;
1632 if(welltemp_) {
1633 double vbias = getBiasAndDerivatives(cv);
1634 if(biasf_>1.0) {
1635 height = height0_*exp(-vbias/(kbt_*(biasf_-1.0)));
1636 } else {
1637 // notice that if gamma=1 we store directly -F
1638 height = height0_*exp(-vbias/kbt_);
1639 }
1640 }
1641 if(dampfactor_>0.0) {
1642 plumed_assert(BiasGrid_);
1643 double m=BiasGrid_->getMaxValue();
1644 height*=exp(-m/(kbt_*(dampfactor_)));
1645 }
1646 if (tt_specs_.is_active) {
1647 double vbarrier = transition_bias_;
1648 temperHeight(height, tt_specs_, vbarrier);
1649 }
1650 if(TargetGrid_) {
1651 double f=TargetGrid_->getValue(cv)-TargetGrid_->getMaxValue();
1652 height*=exp(f/kbt_);
1653 }
1654 return height;
1655 }
1656
temperHeight(double & height,const TemperingSpecs & t_specs,const double tempering_bias)1657 void MetaD::temperHeight(double &height, const TemperingSpecs &t_specs, const double tempering_bias) {
1658 if (t_specs.alpha == 1.0) {
1659 height *= exp(-max(0.0, tempering_bias - t_specs.threshold) / (kbt_ * (t_specs.biasf - 1.0)));
1660 } else {
1661 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 }
1664
calculate()1665 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 if(adaptive_==FlexibleBin::diffusion && getExchangeStep()) error("ADAPTIVE=DIFF is not compatible with replica exchange");
1670
1671 const unsigned ncv=getNumberOfArguments();
1672 vector<double> cv(ncv);
1673 std::unique_ptr<double[]> der(new double[ncv]);
1674 for(unsigned i=0; i<ncv; ++i) {
1675 cv[i]=getArgument(i);
1676 der[i]=0.;
1677 }
1678 double ene = getBiasAndDerivatives(cv,der.get());
1679 // special case for gamma=1.0
1680 if(biasf_==1.0) {
1681 ene=0.0;
1682 for(unsigned i=0; i<getNumberOfArguments(); ++i) {der[i]=0.0;}
1683 }
1684
1685 setBias(ene);
1686 if(calc_rct_) getPntrToComponent("rbias")->set(ene - reweight_factor_);
1687 // calculate the acceleration factor
1688 if(acceleration&&!isFirstStep) {
1689 acc += static_cast<double>(getStride()) * exp(ene/(kbt_));
1690 const double mean_acc = acc/((double) getStep());
1691 getPntrToComponent("acc")->set(mean_acc);
1692 } else if (acceleration && isFirstStep && acc_restart_mean_ > 0.0) {
1693 acc = acc_restart_mean_ * static_cast<double>(getStep());
1694 if(freq_adaptive_) {
1695 // has to be done here if restarting, as the acc is not defined before
1696 updateFrequencyAdaptiveStride();
1697 }
1698 }
1699
1700 getPntrToComponent("work")->set(work_);
1701 // set Forces
1702 for(unsigned i=0; i<ncv; ++i) {
1703 setOutputForce(i,-der[i]);
1704 }
1705 }
1706
update()1707 void MetaD::update() {
1708 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 if(getStep()%current_stride==0 && !isFirstStep )nowAddAHill=true;
1715 else {
1716 nowAddAHill=false;
1717 isFirstStep=false;
1718 }
1719
1720 for(unsigned i=0; i<cv.size(); ++i) cv[i] = getArgument(i);
1721
1722 double vbias=getBiasAndDerivatives(cv);
1723
1724 // if you use adaptive, call the FlexibleBin
1725 if(adaptive_!=FlexibleBin::none) {
1726 flexbin->update(nowAddAHill);
1727 multivariate=true;
1728 } else {
1729 multivariate=false;
1730 }
1731
1732 if(nowAddAHill) {
1733 // add a Gaussian
1734 double height=getHeight(cv);
1735 // returns upper diagonal inverse
1736 if(adaptive_!=FlexibleBin::none) thissigma=flexbin->getInverseMatrix();
1737 // returns normal sigma
1738 else thissigma=sigma0_;
1739
1740 // In case we use walkers_mpi, it is now necessary to communicate with other replicas.
1741 if(walkers_mpi) {
1742 // Allocate arrays to store all walkers hills
1743 std::vector<double> all_cv(mpi_nw_*cv.size(),0.0);
1744 std::vector<double> all_sigma(mpi_nw_*thissigma.size(),0.0);
1745 std::vector<double> all_height(mpi_nw_,0.0);
1746 std::vector<int> all_multivariate(mpi_nw_,0);
1747 if(comm.Get_rank()==0) {
1748 // Communicate (only root)
1749 multi_sim_comm.Allgather(cv,all_cv);
1750 multi_sim_comm.Allgather(thissigma,all_sigma);
1751 // notice that if gamma=1 we store directly -F so this scaling is not necessary:
1752 multi_sim_comm.Allgather(height*(biasf_>1.0?biasf_/(biasf_-1.0):1.0),all_height);
1753 multi_sim_comm.Allgather(int(multivariate),all_multivariate);
1754 }
1755 // Share info with group members
1756 comm.Bcast(all_cv,0);
1757 comm.Bcast(all_sigma,0);
1758 comm.Bcast(all_height,0);
1759 comm.Bcast(all_multivariate,0);
1760
1761 // Flying Gaussian
1762 if (flying) {
1763 hills_.clear();
1764 comm.Barrier();
1765 }
1766
1767 for(unsigned i=0; i<mpi_nw_; i++) {
1768 // actually add hills one by one
1769 std::vector<double> cv_now(cv.size());
1770 std::vector<double> sigma_now(thissigma.size());
1771 for(unsigned j=0; j<cv.size(); j++) cv_now[j]=all_cv[i*cv.size()+j];
1772 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 Gaussian newhill=Gaussian(cv_now,sigma_now,all_height[i]*(biasf_>1.0?(biasf_-1.0)/biasf_:1.0),all_multivariate[i]);
1775 addGaussian(newhill);
1776
1777 // Flying Gaussian
1778 if (!flying) {
1779 writeGaussian(newhill,hillsOfile_);
1780 }
1781
1782 }
1783 } else {
1784 Gaussian newhill=Gaussian(cv,thissigma,height,multivariate);
1785 addGaussian(newhill);
1786 // print on HILLS file
1787 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 if(mw_n_>1 && getStep()%mw_rstride_==0) {
1794 hillsOfile_.flush();
1795 }
1796
1797 double vbias1=getBiasAndDerivatives(cv);
1798 work_+=vbias1-vbias;
1799
1800 // dump grid on file
1801 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 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 int r = 0;
1809 if(walkers_mpi) {
1810 if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
1811 comm.Bcast(r,0);
1812 }
1813 if(r==0) gridfile_.rewind();
1814 }
1815 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 if(!storeOldGrids_) gridfile_.flush();
1823 }
1824
1825 // if multiple walkers and time to read Gaussians
1826 if(mw_n_>1 && getStep()%mw_rstride_==0) {
1827 for(int i=0; i<mw_n_; ++i) {
1828 // don't read your own Gaussians
1829 if(i==mw_id_) continue;
1830 // if the file is not open yet
1831 if(!(ifiles[i]->isOpen())) {
1832 // check if it exists now and open it!
1833 if(ifiles[i]->FileExist(ifilesnames[i])) {
1834 ifiles[i]->open(ifilesnames[i]);
1835 ifiles[i]->reset(false);
1836 }
1837 // otherwise read the new Gaussians
1838 } else {
1839 log.printf(" Reading hills from %s:",ifilesnames[i].c_str());
1840 readGaussians(ifiles[i].get());
1841 ifiles[i]->reset(false);
1842 }
1843 }
1844 }
1845 // Recalculate special bias quantities whenever the bias has been changed by the update.
1846 bool bias_has_changed = (nowAddAHill || (mw_n_ > 1 && getStep() % mw_rstride_ == 0));
1847 if (calc_rct_ && bias_has_changed && getStep()%(stride_*rct_ustride_)==0) computeReweightingFactor();
1848 if (calc_max_bias_ && bias_has_changed) {
1849 max_bias_ = BiasGrid_->getMaxValue();
1850 getPntrToComponent("maxbias")->set(max_bias_);
1851 }
1852 if (calc_transition_bias_ && bias_has_changed) {
1853 transition_bias_ = getTransitionBarrierBias();
1854 getPntrToComponent("transbias")->set(transition_bias_);
1855 }
1856
1857 // Frequency adaptive metadynamics - update hill addition frequency
1858 if(freq_adaptive_ && getStep()%fa_update_frequency_==0) {
1859 updateFrequencyAdaptiveStride();
1860 }
1861
1862 }
1863
1864 /// takes a pointer to the file and a template string with values v and gives back the next center, sigma and height
scanOneHill(IFile * ifile,vector<Value> & tmpvalues,vector<double> & center,vector<double> & sigma,double & height,bool & multivariate)1865 bool MetaD::scanOneHill(IFile *ifile, vector<Value> &tmpvalues, vector<double> ¢er, vector<double> &sigma, double &height, bool &multivariate)
1866 {
1867 double dummy;
1868 multivariate=false;
1869 if(ifile->scanField("time",dummy)) {
1870 unsigned ncv; ncv=tmpvalues.size();
1871 for(unsigned i=0; i<ncv; ++i) {
1872 ifile->scanField( &tmpvalues[i] );
1873 if( tmpvalues[i].isPeriodic() && ! getPntrToArgument(i)->isPeriodic() ) {
1874 error("in hills file periodicity for variable " + tmpvalues[i].getName() + " does not match periodicity in input");
1875 } else if( tmpvalues[i].isPeriodic() ) {
1876 std::string imin, imax; tmpvalues[i].getDomain( imin, imax );
1877 std::string rmin, rmax; getPntrToArgument(i)->getDomain( rmin, rmax );
1878 if( imin!=rmin || imax!=rmax ) {
1879 error("in hills file periodicity for variable " + tmpvalues[i].getName() + " does not match periodicity in input");
1880 }
1881 }
1882 center[i]=tmpvalues[i].get();
1883 }
1884 // scan for kerneltype
1885 std::string ktype="gaussian";
1886 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 ifile->scanField("multivariate",sss);
1890 if(sss=="true") multivariate=true;
1891 else if(sss=="false") multivariate=false;
1892 else plumed_merror("cannot parse multivariate = "+ sss);
1893 if(multivariate) {
1894 sigma.resize(ncv*(ncv+1)/2);
1895 Matrix<double> upper(ncv,ncv);
1896 Matrix<double> lower(ncv,ncv);
1897 for(unsigned i=0; i<ncv; i++) {
1898 for(unsigned j=0; j<ncv-i; j++) {
1899 ifile->scanField("sigma_"+getPntrToArgument(j+i)->getName()+"_"+getPntrToArgument(j)->getName(),lower(j+i,j));
1900 upper(j,j+i)=lower(j+i,j);
1901 }
1902 }
1903 Matrix<double> mymult(ncv,ncv);
1904 Matrix<double> invmatrix(ncv,ncv);
1905 mult(lower,upper,mymult);
1906 // now invert and get the sigmas
1907 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 for(unsigned i=0; i<ncv; i++) {
1911 for(unsigned j=i; j<ncv; j++) {
1912 sigma[k]=invmatrix(i,j);
1913 k++;
1914 }
1915 }
1916 } else {
1917 for(unsigned i=0; i<ncv; ++i) {
1918 ifile->scanField("sigma_"+getPntrToArgument(i)->getName(),sigma[i]);
1919 }
1920 }
1921
1922 ifile->scanField("height",height);
1923 ifile->scanField("biasf",dummy);
1924 if(ifile->FieldExist("clock")) ifile->scanField("clock",dummy);
1925 if(ifile->FieldExist("lower_int")) ifile->scanField("lower_int",dummy);
1926 if(ifile->FieldExist("upper_int")) ifile->scanField("upper_int",dummy);
1927 ifile->scanField();
1928 return true;
1929 } else {
1930 return false;
1931 }
1932 }
1933
computeReweightingFactor()1934 void MetaD::computeReweightingFactor()
1935 {
1936 if(biasf_==1.0) { // in this case we have no bias, so reweight factor is 0.0
1937 getPntrToComponent("rct")->set(0.0);
1938 return;
1939 }
1940
1941 double Z_0=0; //proportional to the integral of exp(-beta*F)
1942 double Z_V=0; //proportional to the integral of exp(-beta*(F+V))
1943 double minusBetaF=biasf_/(biasf_-1.)/kbt_;
1944 double minusBetaFplusV=1./(biasf_-1.)/kbt_;
1945 if (biasf_==-1.0) { //non well-tempered case
1946 minusBetaF=1./kbt_;
1947 minusBetaFplusV=0;
1948 }
1949 max_bias_=BiasGrid_->getMaxValue(); //to avoid exp overflow
1950
1951 const unsigned rank=comm.Get_rank();
1952 const unsigned stride=comm.Get_size();
1953 for (Grid::index_t t=rank; t<BiasGrid_->getSize(); t+=stride) {
1954 const double val=BiasGrid_->getValue(t);
1955 Z_0+=std::exp(minusBetaF*(val-max_bias_));
1956 Z_V+=std::exp(minusBetaFplusV*(val-max_bias_));
1957 }
1958 if (stride>1) {
1959 comm.Sum(Z_0);
1960 comm.Sum(Z_V);
1961 }
1962
1963 reweight_factor_=kbt_*std::log(Z_0/Z_V)+max_bias_;
1964 getPntrToComponent("rct")->set(reweight_factor_);
1965 }
1966
getTransitionBarrierBias()1967 double MetaD::getTransitionBarrierBias() {
1968
1969 // If there is only one well of interest, return the bias at that well point.
1970 if (transitionwells_.size() == 1) {
1971 double tb_bias = getBiasAndDerivatives(transitionwells_[0], NULL);
1972 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 vector<double> sink = transitionwells_[0];
1987 vector<double> source = transitionwells_[1];
1988 least_transition_bias = BiasGrid_->findMaximalPathMinimum(source, sink);
1989 for (unsigned i = 2; i < transitionwells_.size(); i++) {
1990 if (least_transition_bias == 0.0) {
1991 break;
1992 }
1993 source = transitionwells_[i];
1994 double curr_transition_bias = BiasGrid_->findMaximalPathMinimum(source, sink);
1995 least_transition_bias = fmin(curr_transition_bias, least_transition_bias);
1996 }
1997 return least_transition_bias;
1998 }
1999 }
2000
updateFrequencyAdaptiveStride()2001 void MetaD::updateFrequencyAdaptiveStride() {
2002 plumed_massert(freq_adaptive_,"should only be used if frequency adaptive metadynamics is enabled");
2003 plumed_massert(acceleration,"frequency adaptive metadynamics can only be used if the acceleration factor is calculated");
2004 const double mean_acc = acc/((double) getStep());
2005 int tmp_stride= stride_*floor((mean_acc/fa_min_acceleration_)+0.5);
2006 if(mean_acc >= fa_min_acceleration_) {
2007 if(tmp_stride > current_stride) {current_stride = tmp_stride;}
2008 }
2009 if(fa_max_stride_!=0 && current_stride>fa_max_stride_) {
2010 current_stride=fa_max_stride_;
2011 }
2012 getPntrToComponent("pace")->set(current_stride);
2013 }
2014
checkNeedsGradients() const2015 bool MetaD::checkNeedsGradients()const
2016 {
2017 if(adaptive_==FlexibleBin::geometry) {
2018 if(getStep()%stride_==0 && !isFirstStep) return true;
2019 else return false;
2020 } else return false;
2021 }
2022
2023 }
2024 }
2025