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