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