Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "Bias.h"
23 : #include "ActionRegister.h"
24 : #include "core/ActionSet.h"
25 : #include "core/PlumedMain.h"
26 : #include "core/Atoms.h"
27 : #include "core/FlexibleBin.h"
28 : #include "tools/Exception.h"
29 : #include "tools/Grid.h"
30 : #include "tools/Matrix.h"
31 : #include "tools/OpenMP.h"
32 : #include "tools/Random.h"
33 : #include "tools/File.h"
34 : #include <ctime>
35 : #include <numeric>
36 : #if defined(__PLUMED_HAS_GETCWD)
37 : #include <unistd.h>
38 : #endif
39 :
40 : namespace PLMD {
41 : namespace bias {
42 :
43 : //+PLUMEDOC BIAS METAD
44 : /*
45 : Used to performed metadynamics on one or more collective variables.
46 :
47 : In a metadynamics simulations a history dependent bias composed of
48 : intermittently added Gaussian functions is added to the potential \cite metad.
49 :
50 : \f[
51 : V(\vec{s},t) = \sum_{ k \tau < t} W(k \tau)
52 : \exp\left(
53 : -\sum_{i=1}^{d} \frac{(s_i-s_i^{(0)}(k \tau))^2}{2\sigma_i^2}
54 : \right).
55 : \f]
56 :
57 : This potential forces the system away from the kinetic traps in the potential energy surface
58 : and out into the unexplored parts of the energy landscape. Information on the Gaussian
59 : functions from which this potential is composed is output to a file called HILLS, which
60 : is used both the restart the calculation and to reconstruct the free energy as a function of the CVs.
61 : The free energy can be reconstructed from a metadynamics calculation because the final bias is given
62 : by:
63 :
64 : \f[
65 : V(\vec{s}) = -F(\vec{s})
66 : \f]
67 :
68 : During post processing the free energy can be calculated in this way using the \ref sum_hills
69 : utility.
70 :
71 : In the simplest possible implementation of a metadynamics calculation the expense of a metadynamics
72 : calculation increases with the length of the simulation as one has to, at every step, evaluate
73 : the values of a larger and larger number of Gaussian kernels. To avoid this issue you can
74 : store the bias on a grid. This approach is similar to that proposed in \cite babi08jcp but has the
75 : advantage that the grid spacing is independent on the Gaussian width.
76 : Notice that you should provide the grid boundaries (GRID_MIN and GRID_MAX) and either the number of bins
77 : for every collective variable (GRID_BIN) or the desired grid spacing (GRID_SPACING).
78 : In case you provide both PLUMED will use the most conservative choice (highest number of bins) for each dimension.
79 : In case you do not provide any information about bin size (neither GRID_BIN nor GRID_SPACING)
80 : PLUMED will use 1/5 of the Gaussian width (SIGMA) as grid spacing if the width is fixed or 1/5 of the minimum
81 : Gaussian width (SIGMA_MIN) if the width is variable. This default choice should be reasonable for most applications.
82 :
83 : Alternatively to the use of grids, it is possible to use a neighbor list to decrease the cost of evaluating the bias,
84 : this can be enabled using NLIST. NLIST can be beneficial with more than 2 collective variables, where GRID becomes
85 : expensive and memory consuming. The neighbor list will be updated everytime the CVs go farther than a cut-off value
86 : from the position they were at last neighbor list update. Gaussians are added to the neigbhor list if their center
87 : is within 6.*DP2CUTOFF*sigma*sigma. While the list is updated if the CVs are farther from the center than 0.5 of the
88 : standard deviation of the Gaussian center distribution of the list. These parameters (6 and 0.5) can be modified using
89 : NLIST_PARAMETERS. Note that the use of neighbor list does not provide the exact bias.
90 :
91 : Metadynamics can be restarted either from a HILLS file as well as from a GRID, in this second
92 : case one can first save a GRID using GRID_WFILE (and GRID_WSTRIDE) and at a later stage read
93 : it using GRID_RFILE.
94 :
95 : The work performed by the METAD bias can be calculated using CALC_WORK, note that this is expensive when not using grids.
96 :
97 : Another option that is available in plumed is well-tempered metadynamics \cite Barducci:2008. In this
98 : variant of metadynamics the heights of the Gaussian hills are scaled at each step so the bias is now
99 : given by:
100 :
101 : \f[
102 : 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(
103 : -\sum_{i=1}^{d} \frac{(s_i({q})-s_i({q}(t'))^2}{2\sigma_i^2}
104 : \right),
105 : \f]
106 :
107 : This method ensures that the bias converges more smoothly. It should be noted that, in the case of well-tempered metadynamics, in
108 : the output printed the Gaussian height is re-scaled using the bias factor.
109 : Also notice that with well-tempered metadynamics the HILLS file does not contain the bias,
110 : but the negative of the free-energy estimate. This choice has the advantage that
111 : one can restart a simulation using a different value for the \f$\Delta T\f$. The applied bias will be scaled accordingly.
112 :
113 : Note that you can use here also the flexible Gaussian approach \cite Branduardi:2012dl
114 : in which you can adapt the Gaussian to the extent of Cartesian space covered by a variable or
115 : to the space in collective variable covered in a given time. In this case the width of the deposited
116 : Gaussian potential is denoted by one value only that is a Cartesian space (ADAPTIVE=GEOM) or a time
117 : (ADAPTIVE=DIFF). Note that a specific integration technique for the deposited Gaussian kernels
118 : should be used in this case. Check the documentation for utility sum_hills.
119 :
120 : With the keyword INTERVAL one changes the metadynamics algorithm setting the bias force equal to zero
121 : outside boundary \cite baftizadeh2012protein. If, for example, metadynamics is performed on a CV s and one is interested only
122 : to the free energy for s > boundary, the history dependent potential is still updated according to the above
123 : equations but the metadynamics force is set to zero for s < boundary. Notice that Gaussian kernels are added also
124 : if s < boundary, as the tails of these Gaussian kernels influence VG in the relevant region s > boundary. In this way, the
125 : force on the system in the region s > boundary comes from both metadynamics and the force field, in the region
126 : s < boundary only from the latter. This approach allows obtaining a history-dependent bias potential VG that
127 : fluctuates around a stable estimator, equal to the negative of the free energy far enough from the
128 : boundaries. Note that:
129 : - It works only for one-dimensional biases;
130 : - It works both with and without GRID;
131 : - The interval limit boundary in a region where the free energy derivative is not large;
132 : - If in the region outside the limit boundary the system has a free energy minimum, the INTERVAL keyword should
133 : be used together with a \ref UPPER_WALLS or \ref LOWER_WALLS at boundary.
134 :
135 : As a final note, since version 2.0.2 when the system is outside of the selected interval the force
136 : is set to zero and the bias value to the value at the corresponding boundary. This allows acceptances
137 : for replica exchange methods to be computed correctly.
138 :
139 : Multiple walkers \cite multiplewalkers can also be used. See below the examples.
140 :
141 :
142 : The \f$c(t)\f$ reweighting factor can also be calculated on the fly using the equations
143 : presented in \cite Tiwary_jp504920s.
144 : The expression used to calculate \f$c(t)\f$ follows directly from Eq. 3 in \cite Tiwary_jp504920s,
145 : where \f$F(\vec{s})=-\gamma/(\gamma-1) V(\vec{s})\f$.
146 : This gives smoother results than equivalent Eqs. 13 and Eqs. 14 in that paper.
147 : The \f$c(t)\f$ is given by the rct component while the bias
148 : normalized by \f$c(t)\f$ is given by the rbias component (rbias=bias-rct) which can be used
149 : to obtain a reweighted histogram.
150 : The calculation of \f$c(t)\f$ is enabled by using the keyword CALC_RCT.
151 : By default \f$c(t)\f$ is updated every time the bias changes, but if this slows down the simulation
152 : the keyword RCT_USTRIDE can be set to a value higher than 1.
153 : This option requires that a grid is used.
154 :
155 : Additional material and examples can be also found in the tutorials:
156 :
157 : - \ref lugano-3
158 :
159 : Concurrent metadynamics
160 : as done e.g. in Ref. \cite gil2015enhanced . This indeed can be obtained by using the METAD
161 : action multiple times in the same input file.
162 :
163 : \par Examples
164 :
165 : The following input is for a standard metadynamics calculation using as
166 : collective variables the distance between atoms 3 and 5
167 : and the distance between atoms 2 and 4. The value of the CVs and
168 : the metadynamics bias potential are written to the COLVAR file every 100 steps.
169 : \plumedfile
170 : DISTANCE ATOMS=3,5 LABEL=d1
171 : DISTANCE ATOMS=2,4 LABEL=d2
172 : METAD ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3 PACE=500 LABEL=restraint
173 : PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
174 : \endplumedfile
175 : (See also \ref DISTANCE \ref PRINT).
176 :
177 : \par
178 : If you use adaptive Gaussian kernels, with diffusion scheme where you use
179 : a Gaussian that should cover the space of 20 time steps in collective variables.
180 : Note that in this case the histogram correction is needed when summing up hills.
181 : \plumedfile
182 : DISTANCE ATOMS=3,5 LABEL=d1
183 : DISTANCE ATOMS=2,4 LABEL=d2
184 : METAD ARG=d1,d2 SIGMA=20 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=DIFF
185 : PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
186 : \endplumedfile
187 :
188 : \par
189 : If you use adaptive Gaussian kernels, with geometrical scheme where you use
190 : a Gaussian that should cover the space of 0.05 nm in Cartesian space.
191 : Note that in this case the histogram correction is needed when summing up hills.
192 : \plumedfile
193 : DISTANCE ATOMS=3,5 LABEL=d1
194 : DISTANCE ATOMS=2,4 LABEL=d2
195 : METAD ARG=d1,d2 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=GEOM
196 : PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
197 : \endplumedfile
198 :
199 : \par
200 : When using adaptive Gaussian kernels you might want to limit how the hills width can change.
201 : You can use SIGMA_MIN and SIGMA_MAX keywords.
202 : The sigmas should specified in terms of CV so you should use the CV units.
203 : Note that if you use a negative number, this means that the limit is not set.
204 : Note also that in this case the histogram correction is needed when summing up hills.
205 : \plumedfile
206 : DISTANCE ATOMS=3,5 LABEL=d1
207 : DISTANCE ATOMS=2,4 LABEL=d2
208 : METAD ...
209 : ARG=d1,d2 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=GEOM
210 : SIGMA_MIN=0.2,0.1 SIGMA_MAX=0.5,1.0
211 : ... METAD
212 : PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
213 : \endplumedfile
214 :
215 : \par
216 : Multiple walkers can be also use as in \cite multiplewalkers
217 : These are enabled by setting the number of walker used, the id of the
218 : current walker which interprets the input file, the directory where the
219 : hills containing files resides, and the frequency to read the other walkers.
220 : Here is an example
221 : \plumedfile
222 : DISTANCE ATOMS=3,5 LABEL=d1
223 : METAD ...
224 : ARG=d1 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint
225 : WALKERS_N=10
226 : WALKERS_ID=3
227 : WALKERS_DIR=../
228 : WALKERS_RSTRIDE=100
229 : ... METAD
230 : \endplumedfile
231 : where WALKERS_N is the total number of walkers, WALKERS_ID is the
232 : id of the present walker (starting from 0 ) and the WALKERS_DIR is the directory
233 : where all the walkers are located. WALKERS_RSTRIDE is the number of step between
234 : one update and the other. Since version 2.2.5, hills files are automatically
235 : flushed every WALKERS_RSTRIDE steps.
236 :
237 : \par
238 : The \f$c(t)\f$ reweighting factor can be calculated on the fly using the equations
239 : presented in \cite Tiwary_jp504920s as described above.
240 : This is enabled by using the keyword CALC_RCT,
241 : and can be done only if the bias is defined on a grid.
242 : \plumedfile
243 : phi: TORSION ATOMS=1,2,3,4
244 : psi: TORSION ATOMS=5,6,7,8
245 :
246 : METAD ...
247 : LABEL=metad
248 : ARG=phi,psi SIGMA=0.20,0.20 HEIGHT=1.20 BIASFACTOR=5 TEMP=300.0 PACE=500
249 : GRID_MIN=-pi,-pi GRID_MAX=pi,pi GRID_BIN=150,150
250 : CALC_RCT
251 : RCT_USTRIDE=10
252 : ... METAD
253 : \endplumedfile
254 : Here we have asked that the calculation is performed every 10 hills deposition by using
255 : RCT_USTRIDE keyword. If this keyword is not given, the calculation will
256 : by default be performed every time the bias changes. The \f$c(t)\f$ reweighting factor will be given
257 : in the rct component while the instantaneous value of the bias potential
258 : normalized using the \f$c(t)\f$ reweighting factor is given in the rbias component
259 : [rbias=bias-rct] which can be used to obtain a reweighted histogram or
260 : free energy surface using the \ref HISTOGRAM analysis.
261 :
262 : \par
263 : The kinetics of the transitions between basins can also be analyzed on the fly as
264 : in \cite PRL230602. The flag ACCELERATION turn on accumulation of the acceleration
265 : factor that can then be used to determine the rate. This method can be used together
266 : with \ref COMMITTOR analysis to stop the simulation when the system get to the target basin.
267 : It must be used together with Well-Tempered Metadynamics. If restarting from a previous
268 : metadynamics you need to use the ACCELERATION_RFILE keyword to give the name of the
269 : data file from which the previous value of the acceleration factor should be read, otherwise the
270 : calculation of the acceleration factor will be wrong.
271 :
272 : \par
273 : By using the flag FREQUENCY_ADAPTIVE the frequency adaptive scheme introduced in \cite Wang-JCP-2018
274 : is turned on. The frequency for hill addition then changes dynamically based on the acceleration factor
275 : according to the following equation
276 : \f[
277 : \tau_{\mathrm{dep}}(t) =
278 : \min\left[
279 : \tau_0 \cdot
280 : \max\left[\frac{\alpha(t)}{\theta},1\right]
281 : ,\tau_{c}
282 : \right]
283 : \f]
284 : where \f$\tau_0\f$ is the initial hill addition frequency given by the PACE keyword,
285 : \f$\tau_{c}\f$ is the maximum allowed frequency given by the FA_MAX_PACE keyword,
286 : \f$\alpha(t)\f$ is the instantaneous acceleration factor at time \f$t\f$,
287 : and \f$\theta\f$ is a threshold value that acceleration factor has to reach before
288 : triggering a change in the hill addition frequency given by the FA_MIN_ACCELERATION keyword.
289 : The frequency for updating the hill addition frequency according to this equation is
290 : given by the FA_UPDATE_FREQUENCY keyword, by default it is the same as the value given
291 : in PACE. The hill hill addition frequency increase monotonously such that if the
292 : instantaneous acceleration factor is lower than in the previous updating step the
293 : previous \f$\tau_{\mathrm{dep}}\f$ is kept rather than updating it to a lower value.
294 : The instantaneous hill addition frequency \f$\tau_{\mathrm{dep}}(t)\f$ is outputted
295 : to pace component. Note that if restarting from a previous metadynamics run you need to
296 : use the ACCELERATION_RFILE keyword to read in the acceleration factors from the
297 : previous run, otherwise the hill addition frequency will start from the initial
298 : frequency.
299 :
300 :
301 : \par
302 : You can also provide a target distribution using the keyword TARGET
303 : \cite white2015designing
304 : \cite marinelli2015ensemble
305 : \cite gil2016empirical
306 : The TARGET should be a grid containing a free-energy (i.e. the -\f$k_B\f$T*log of the desired target distribution).
307 : Gaussian kernels will then be scaled by a factor
308 : \f[
309 : e^{\beta(\tilde{F}(s)-\tilde{F}_{max})}
310 : \f]
311 : Here \f$\tilde{F}(s)\f$ is the free energy defined on the grid and \f$\tilde{F}_{max}\f$ its maximum value.
312 : Notice that we here used the maximum value as in ref \cite gil2016empirical
313 : This choice allows to avoid exceedingly large Gaussian kernels to be added. However,
314 : it could make the Gaussian too small. You should always choose carefully the HEIGHT parameter
315 : in this case.
316 : The grid file should be similar to other PLUMED grid files in that it should contain
317 : both the target free-energy and its derivatives.
318 :
319 : Notice that if you wish your simulation to converge to the target free energy you should use
320 : the DAMPFACTOR command to provide a global tempering \cite dama2014well
321 : Alternatively, if you use a BIASFACTOR your simulation will converge to a free
322 : energy that is a linear combination of the target free energy and of the intrinsic free energy
323 : determined by the original force field.
324 :
325 : \plumedfile
326 : DISTANCE ATOMS=3,5 LABEL=d1
327 : METAD ...
328 : LABEL=t1
329 : ARG=d1 SIGMA=0.05 TAU=200 DAMPFACTOR=100 PACE=250
330 : GRID_MIN=1.14 GRID_MAX=1.32 GRID_BIN=6
331 : TARGET=dist.grid
332 : ... METAD
333 :
334 : PRINT ARG=d1,t1.bias STRIDE=100 FILE=COLVAR
335 : \endplumedfile
336 :
337 : The file dist.dat for this calculation would read:
338 :
339 : \auxfile{dist.grid}
340 : #! FIELDS d1 t1.target der_d1
341 : #! SET min_d1 1.14
342 : #! SET max_d1 1.32
343 : #! SET nbins_d1 6
344 : #! SET periodic_d1 false
345 : 1.1400 0.0031 0.1101
346 : 1.1700 0.0086 0.2842
347 : 1.2000 0.0222 0.6648
348 : 1.2300 0.0521 1.4068
349 : 1.2600 0.1120 2.6873
350 : 1.2900 0.2199 4.6183
351 : 1.3200 0.3948 7.1055
352 : \endauxfile
353 :
354 : Notice that BIASFACTOR can also be chosen as equal to 1. In this case one will perform
355 : unbiased sampling. Instead of using HEIGHT, one should provide the TAU parameter.
356 : \plumedfile
357 : d: DISTANCE ATOMS=3,5
358 : METAD ARG=d SIGMA=0.1 TAU=4.0 TEMP=300 PACE=100 BIASFACTOR=1.0
359 : \endplumedfile
360 : The HILLS file obtained will still work with `plumed sum_hills` so as to plot a free-energy.
361 : The case where this makes sense is probably that of RECT simulations.
362 :
363 : Regarding RECT simulations, you can also use the RECT keyword so as to avoid using multiple input files.
364 : For instance, a single input file will be
365 : \plumedfile
366 : d: DISTANCE ATOMS=3,5
367 : METAD ARG=d SIGMA=0.1 TAU=4.0 TEMP=300 PACE=100 RECT=1.0,1.5,2.0,3.0
368 : \endplumedfile
369 : The number of elements in the RECT array should be equal to the number of replicas.
370 :
371 : */
372 : //+ENDPLUMEDOC
373 :
374 : class MetaD : public Bias {
375 :
376 : private:
377 : struct Gaussian {
378 : bool multivariate; // this is required to discriminate the one dimensional case
379 : double height;
380 : std::vector<double> center;
381 : std::vector<double> sigma;
382 : std::vector<double> invsigma;
383 5617 : Gaussian(const bool m, const double h, const std::vector<double>& c, const std::vector<double>& s):
384 5617 : multivariate(m),height(h),center(c),sigma(s),invsigma(s) {
385 : // to avoid troubles from zero element in flexible hills
386 16388 : for(unsigned i=0; i<invsigma.size(); ++i) {
387 10771 : if(std::abs(invsigma[i])>1.e-20) {
388 10771 : invsigma[i]=1.0/invsigma[i] ;
389 : } else {
390 0 : invsigma[i]=0.0;
391 : }
392 : }
393 5617 : }
394 : };
395 8 : struct TemperingSpecs {
396 : bool is_active;
397 : std::string name_stem;
398 : std::string name;
399 : double biasf;
400 : double threshold;
401 : double alpha;
402 156 : inline TemperingSpecs(bool is_active, const std::string &name_stem, const std::string &name, double biasf, double threshold, double alpha) :
403 156 : is_active(is_active), name_stem(name_stem), name(name), biasf(biasf), threshold(threshold), alpha(alpha)
404 156 : {}
405 : };
406 : // general setup
407 : double kbt_;
408 : int stride_;
409 : bool calc_work_;
410 : // well-tempered MetaD
411 : bool welltemp_;
412 : double biasf_;
413 : // output files format
414 : std::string fmt_;
415 : // first step?
416 : bool isFirstStep_;
417 : // Gaussian starting parameters
418 : double height0_;
419 : std::vector<double> sigma0_;
420 : std::vector<double> sigma0min_;
421 : std::vector<double> sigma0max_;
422 : // Gaussians
423 : std::vector<Gaussian> hills_;
424 : std::unique_ptr<FlexibleBin> flexbin_;
425 : int adaptive_;
426 : OFile hillsOfile_;
427 : std::vector<std::unique_ptr<IFile>> ifiles_;
428 : std::vector<std::string> ifilesnames_;
429 : // Grids
430 : bool grid_;
431 : std::unique_ptr<GridBase> BiasGrid_;
432 : OFile gridfile_;
433 : bool storeOldGrids_;
434 : int wgridstride_;
435 : // multiple walkers
436 : int mw_n_;
437 : std::string mw_dir_;
438 : int mw_id_;
439 : int mw_rstride_;
440 : bool walkers_mpi_;
441 : unsigned mpi_nw_;
442 : // flying gaussians
443 : bool flying_;
444 : // kinetics from metadynamics
445 : bool acceleration_;
446 : double acc_;
447 : double acc_restart_mean_;
448 : // transition-tempering metadynamics
449 : bool calc_max_bias_;
450 : double max_bias_;
451 : bool calc_transition_bias_;
452 : double transition_bias_;
453 : std::vector<std::vector<double> > transitionwells_;
454 : static const size_t n_tempering_options_ = 1;
455 : static const std::string tempering_names_[1][2];
456 : double dampfactor_;
457 : struct TemperingSpecs tt_specs_;
458 : std::string targetfilename_;
459 : std::unique_ptr<GridBase> TargetGrid_;
460 : // frequency adaptive metadynamics
461 : int current_stride_;
462 : bool freq_adaptive_;
463 : int fa_update_frequency_;
464 : int fa_max_stride_;
465 : double fa_min_acceleration_;
466 : // intervals
467 : double uppI_;
468 : double lowI_;
469 : bool doInt_;
470 : // reweighting
471 : bool calc_rct_;
472 : double reweight_factor_;
473 : unsigned rct_ustride_;
474 : // work
475 : double work_;
476 : // neighbour list stuff
477 : bool nlist_;
478 : bool nlist_update_;
479 : unsigned nlist_steps_;
480 : std::array<double,2> nlist_param_;
481 : std::vector<Gaussian> nlist_hills_;
482 : std::vector<double> nlist_center_;
483 : std::vector<double> nlist_dev2_;
484 :
485 : double stretchA=1.0;
486 : double stretchB=0.0;
487 :
488 : bool noStretchWarningDone=false;
489 :
490 12 : void noStretchWarning() {
491 12 : if(!noStretchWarningDone) {
492 3 : log<<"\nWARNING: you are using a HILLS file with Gaussian kernels, PLUMED 2.8 uses stretched Gaussians by default\n";
493 : }
494 12 : noStretchWarningDone=true;
495 12 : }
496 :
497 : static void registerTemperingKeywords(const std::string &name_stem, const std::string &name, Keywords &keys);
498 : void readTemperingSpecs(TemperingSpecs &t_specs);
499 : void logTemperingSpecs(const TemperingSpecs &t_specs);
500 : void readGaussians(IFile*);
501 : void writeGaussian(const Gaussian&,OFile&);
502 : void addGaussian(const Gaussian&);
503 : double getHeight(const std::vector<double>&);
504 : void temperHeight(double &height, const TemperingSpecs &t_specs, const double tempering_bias);
505 : double getBias(const std::vector<double>&);
506 : double getBiasAndDerivatives(const std::vector<double>&, std::vector<double>&);
507 : double evaluateGaussian(const std::vector<double>&, const Gaussian&);
508 : double evaluateGaussianAndDerivatives(const std::vector<double>&, const Gaussian&,std::vector<double>&,std::vector<double>&);
509 : double getGaussianNormalization(const Gaussian&);
510 : std::vector<unsigned> getGaussianSupport(const Gaussian&);
511 : bool scanOneHill(IFile* ifile, std::vector<Value>& v, std::vector<double>& center, std::vector<double>& sigma, double& height, bool& multivariate);
512 : void computeReweightingFactor();
513 : double getTransitionBarrierBias();
514 : void updateFrequencyAdaptiveStride();
515 : void updateNlist();
516 :
517 : public:
518 : explicit MetaD(const ActionOptions&);
519 : void calculate() override;
520 : void update() override;
521 : static void registerKeywords(Keywords& keys);
522 : bool checkNeedsGradients()const override;
523 : };
524 :
525 14090 : PLUMED_REGISTER_ACTION(MetaD,"METAD")
526 :
527 161 : void MetaD::registerKeywords(Keywords& keys) {
528 161 : Bias::registerKeywords(keys);
529 322 : keys.addOutputComponent("rbias","CALC_RCT","the instantaneous value of the bias normalized using the c(t) reweighting factor [rbias=bias-rct]."
530 : "This component can be used to obtain a reweighted histogram.");
531 322 : keys.addOutputComponent("rct","CALC_RCT","the reweighting factor \\f$c(t)\\f$.");
532 322 : keys.addOutputComponent("work","CALC_WORK","accumulator for work");
533 322 : keys.addOutputComponent("acc","ACCELERATION","the metadynamics acceleration factor");
534 322 : keys.addOutputComponent("maxbias", "CALC_MAX_BIAS", "the maximum of the metadynamics V(s, t)");
535 322 : keys.addOutputComponent("transbias", "CALC_TRANSITION_BIAS", "the metadynamics transition bias V*(t)");
536 322 : keys.addOutputComponent("pace","FREQUENCY_ADAPTIVE","the hill addition frequency when employing frequency adaptive metadynamics");
537 322 : keys.addOutputComponent("nlker","NLIST","number of hills in the neighbor list");
538 322 : keys.addOutputComponent("nlsteps","NLIST","number of steps from last neighbor list update");
539 161 : keys.use("ARG");
540 322 : keys.add("compulsory","SIGMA","the widths of the Gaussian hills");
541 322 : keys.add("compulsory","PACE","the frequency for hill addition");
542 322 : keys.add("compulsory","FILE","HILLS","a file in which the list of added hills is stored");
543 322 : keys.add("optional","HEIGHT","the heights of the Gaussian hills. Compulsory unless TAU and either BIASFACTOR or DAMPFACTOR are given");
544 322 : keys.add("optional","FMT","specify format for HILLS files (useful for decrease the number of digits in regtests)");
545 322 : keys.add("optional","BIASFACTOR","use well tempered metadynamics and use this bias factor. Please note you must also specify temp");
546 322 : keys.addFlag("CALC_WORK",false,"calculate the total accumulated work done by the bias since last restart");
547 322 : keys.add("optional","RECT","list of bias factors for all the replicas");
548 322 : keys.add("optional","DAMPFACTOR","damp hills with exp(-max(V)/(kT*DAMPFACTOR)");
549 322 : for (size_t i = 0; i < n_tempering_options_; i++) {
550 161 : registerTemperingKeywords(tempering_names_[i][0], tempering_names_[i][1], keys);
551 : }
552 322 : keys.add("optional","TARGET","target to a predefined distribution");
553 322 : keys.add("optional","TEMP","the system temperature - this is only needed if you are doing well-tempered metadynamics");
554 322 : keys.add("optional","TAU","in well tempered metadynamics, sets height to (k_B Delta T*pace*timestep)/tau");
555 322 : keys.addFlag("CALC_RCT",false,"calculate the c(t) reweighting factor and use that to obtain the normalized bias [rbias=bias-rct]."
556 : "This method is not compatible with metadynamics not on a grid.");
557 322 : keys.add("optional","RCT_USTRIDE","the update stride for calculating the \\f$c(t)\\f$ reweighting factor."
558 : "The default 1, so \\f$c(t)\\f$ is updated every time the bias is updated.");
559 322 : keys.add("optional","GRID_MIN","the lower bounds for the grid");
560 322 : keys.add("optional","GRID_MAX","the upper bounds for the grid");
561 322 : keys.add("optional","GRID_BIN","the number of bins for the grid");
562 322 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
563 322 : keys.addFlag("GRID_SPARSE",false,"use a sparse grid to store hills");
564 322 : keys.addFlag("GRID_NOSPLINE",false,"don't use spline interpolation with grids");
565 322 : keys.add("optional","GRID_WSTRIDE","write the grid to a file every N steps");
566 322 : keys.add("optional","GRID_WFILE","the file on which to write the grid");
567 322 : keys.add("optional","GRID_RFILE","a grid file from which the bias should be read at the initial step of the simulation");
568 322 : keys.addFlag("STORE_GRIDS",false,"store all the grid files the calculation generates. They will be deleted if this keyword is not present");
569 322 : keys.addFlag("NLIST",false,"Use neighbor list for kernels summation, faster but experimental");
570 322 : keys.add("optional", "NLIST_PARAMETERS","(default=6.,0.5) the two cutoff parameters for the Gaussians neighbor list");
571 322 : 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");
572 322 : keys.add("optional","SIGMA_MAX","the upper bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
573 322 : keys.add("optional","SIGMA_MIN","the lower bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
574 322 : keys.add("optional","WALKERS_ID", "walker id");
575 322 : keys.add("optional","WALKERS_N", "number of walkers");
576 322 : keys.add("optional","WALKERS_DIR", "shared directory with the hills files from all the walkers");
577 322 : keys.add("optional","WALKERS_RSTRIDE","stride for reading hills files");
578 322 : keys.addFlag("WALKERS_MPI",false,"Switch on MPI version of multiple walkers - not compatible with WALKERS_* options other than WALKERS_DIR");
579 322 : keys.add("optional","INTERVAL","one dimensional lower and upper limits, outside the limits the system will not feel the biasing force.");
580 322 : keys.addFlag("FLYING_GAUSSIAN",false,"Switch on flying Gaussian method, must be used with WALKERS_MPI");
581 322 : keys.addFlag("ACCELERATION",false,"Set to TRUE if you want to compute the metadynamics acceleration factor.");
582 322 : keys.add("optional","ACCELERATION_RFILE","a data file from which the acceleration should be read at the initial step of the simulation");
583 322 : keys.addFlag("CALC_MAX_BIAS", false, "Set to TRUE if you want to compute the maximum of the metadynamics V(s, t)");
584 322 : keys.addFlag("CALC_TRANSITION_BIAS", false, "Set to TRUE if you want to compute a metadynamics transition bias V*(t)");
585 322 : 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.");
586 322 : 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.");
587 322 : 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");
588 322 : keys.add("optional","FA_MAX_PACE","the maximum hill addition frequency allowed in frequency adaptive metadynamics. By default there is no maximum value.");
589 322 : 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.");
590 161 : keys.use("RESTART");
591 161 : keys.use("UPDATE_FROM");
592 161 : keys.use("UPDATE_UNTIL");
593 161 : }
594 :
595 : const std::string MetaD::tempering_names_[1][2] = {{"TT", "transition tempered"}};
596 :
597 161 : void MetaD::registerTemperingKeywords(const std::string &name_stem, const std::string &name, Keywords &keys) {
598 322 : keys.add("optional", name_stem + "BIASFACTOR", "use " + name + " metadynamics with this bias factor. Please note you must also specify temp");
599 322 : keys.add("optional", name_stem + "BIASTHRESHOLD", "use " + name + " metadynamics with this bias threshold. Please note you must also specify " + name_stem + "BIASFACTOR");
600 322 : 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");
601 161 : }
602 :
603 157 : MetaD::MetaD(const ActionOptions& ao):
604 : PLUMED_BIAS_INIT(ao),
605 156 : kbt_(0.0),
606 156 : stride_(0),
607 156 : calc_work_(false),
608 156 : welltemp_(false),
609 156 : biasf_(-1.0),
610 156 : isFirstStep_(true),
611 156 : height0_(std::numeric_limits<double>::max()),
612 156 : adaptive_(FlexibleBin::none),
613 156 : grid_(false),
614 156 : wgridstride_(0),
615 156 : mw_n_(1), mw_dir_(""), mw_id_(0), mw_rstride_(1),
616 156 : walkers_mpi_(false), mpi_nw_(0),
617 156 : flying_(false),
618 156 : acceleration_(false), acc_(0.0), acc_restart_mean_(0.0),
619 156 : calc_max_bias_(false), max_bias_(0.0),
620 156 : calc_transition_bias_(false), transition_bias_(0.0),
621 156 : dampfactor_(0.0),
622 312 : tt_specs_(false, "TT", "Transition Tempered", -1.0, 0.0, 1.0),
623 156 : current_stride_(0),
624 156 : freq_adaptive_(false),
625 156 : fa_update_frequency_(0),
626 156 : fa_max_stride_(0),
627 156 : fa_min_acceleration_(1.0),
628 156 : uppI_(-1), lowI_(-1), doInt_(false),
629 156 : calc_rct_(false),
630 156 : reweight_factor_(0.0),
631 156 : rct_ustride_(1),
632 156 : work_(0),
633 156 : nlist_(false),
634 156 : nlist_update_(false),
635 469 : nlist_steps_(0) {
636 156 : if(!dp2cutoffNoStretch()) {
637 156 : stretchA=dp2cutoffA;
638 156 : stretchB=dp2cutoffB;
639 : }
640 : // parse the flexible hills
641 : std::string adaptiveoption;
642 : adaptiveoption="NONE";
643 312 : parse("ADAPTIVE",adaptiveoption);
644 156 : if(adaptiveoption=="GEOM") {
645 22 : log.printf(" Uses Geometry-based hills width: sigma must be in distance units and only one sigma is needed\n");
646 22 : adaptive_=FlexibleBin::geometry;
647 134 : } else if(adaptiveoption=="DIFF") {
648 3 : log.printf(" Uses Diffusion-based hills width: sigma must be in time steps and only one sigma is needed\n");
649 3 : adaptive_=FlexibleBin::diffusion;
650 131 : } else if(adaptiveoption=="NONE") {
651 130 : adaptive_=FlexibleBin::none;
652 : } else {
653 1 : error("I do not know this type of adaptive scheme");
654 : }
655 :
656 155 : parse("FMT",fmt_);
657 :
658 : // parse the sigma
659 155 : parseVector("SIGMA",sigma0_);
660 155 : if(adaptive_==FlexibleBin::none) {
661 : // if you use normal sigma you need one sigma per argument
662 130 : if( sigma0_.size()!=getNumberOfArguments() ) {
663 0 : error("number of arguments does not match number of SIGMA parameters");
664 : }
665 : } else {
666 : // if you use flexible hills you need one sigma
667 25 : if(sigma0_.size()!=1) {
668 1 : error("If you choose ADAPTIVE you need only one sigma according to your choice of type (GEOM/DIFF)");
669 : }
670 : // if adaptive then the number must be an integer
671 24 : if(adaptive_==FlexibleBin::diffusion) {
672 3 : if(int(sigma0_[0])-sigma0_[0]>1.e-9 || int(sigma0_[0])-sigma0_[0] <-1.e-9 || int(sigma0_[0])<1 ) {
673 0 : error("In case of adaptive hills with diffusion, the sigma must be an integer which is the number of time steps\n");
674 : }
675 : }
676 : // here evtl parse the sigma min and max values
677 48 : parseVector("SIGMA_MIN",sigma0min_);
678 24 : if(sigma0min_.size()>0 && sigma0min_.size()!=getNumberOfArguments()) {
679 1 : error("the number of SIGMA_MIN values be the same of the number of the arguments");
680 23 : } else if(sigma0min_.size()==0) {
681 23 : sigma0min_.resize(getNumberOfArguments());
682 67 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
683 44 : sigma0min_[i]=-1.;
684 : }
685 : }
686 :
687 46 : parseVector("SIGMA_MAX",sigma0max_);
688 23 : if(sigma0max_.size()>0 && sigma0max_.size()!=getNumberOfArguments()) {
689 1 : error("the number of SIGMA_MAX values be the same of the number of the arguments");
690 22 : } else if(sigma0max_.size()==0) {
691 22 : sigma0max_.resize(getNumberOfArguments());
692 64 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
693 42 : sigma0max_[i]=-1.;
694 : }
695 : }
696 :
697 44 : flexbin_=Tools::make_unique<FlexibleBin>(adaptive_,this,sigma0_[0],sigma0min_,sigma0max_);
698 : }
699 :
700 : // note: HEIGHT is not compulsory, since one could use the TAU keyword, see below
701 152 : parse("HEIGHT",height0_);
702 152 : parse("PACE",stride_);
703 151 : if(stride_<=0) {
704 0 : error("frequency for hill addition is nonsensical");
705 : }
706 151 : current_stride_ = stride_;
707 159 : std::string hillsfname="HILLS";
708 151 : parse("FILE",hillsfname);
709 :
710 : // Manually set to calculate special bias quantities
711 : // throughout the course of simulation. (These are chosen due to
712 : // relevance for tempering and event-driven logic as well.)
713 151 : parseFlag("CALC_MAX_BIAS", calc_max_bias_);
714 305 : parseFlag("CALC_TRANSITION_BIAS", calc_transition_bias_);
715 :
716 : std::vector<double> rect_biasf_;
717 302 : parseVector("RECT",rect_biasf_);
718 151 : if(rect_biasf_.size()>0) {
719 18 : int r=0;
720 18 : if(comm.Get_rank()==0) {
721 9 : r=multi_sim_comm.Get_rank();
722 : }
723 18 : comm.Bcast(r,0);
724 18 : biasf_=rect_biasf_[r];
725 18 : log<<" You are using RECT\n";
726 : } else {
727 266 : parse("BIASFACTOR",biasf_);
728 : }
729 151 : if( biasf_<1.0 && biasf_!=-1.0) {
730 0 : error("well tempered bias factor is nonsensical");
731 : }
732 151 : parse("DAMPFACTOR",dampfactor_);
733 151 : double temp=0.0;
734 151 : parse("TEMP",temp);
735 151 : if(temp>0.0) {
736 56 : kbt_=plumed.getAtoms().getKBoltzmann()*temp;
737 : } else {
738 95 : kbt_=plumed.getAtoms().getKbT();
739 : }
740 151 : if(biasf_>=1.0) {
741 38 : if(kbt_==0.0) {
742 0 : error("Unless the MD engine passes the temperature to plumed, with well-tempered metad you must specify it using TEMP");
743 : }
744 38 : welltemp_=true;
745 : }
746 151 : if(dampfactor_>0.0) {
747 2 : if(kbt_==0.0) {
748 0 : error("Unless the MD engine passes the temperature to plumed, with damped metad you must specify it using TEMP");
749 : }
750 : }
751 :
752 151 : parseFlag("CALC_WORK",calc_work_);
753 :
754 : // Set transition tempering parameters.
755 : // Transition wells are read later via calc_transition_bias_.
756 151 : readTemperingSpecs(tt_specs_);
757 151 : if (tt_specs_.is_active) {
758 3 : calc_transition_bias_ = true;
759 : }
760 :
761 : // If any previous option specified to calculate a transition bias,
762 : // now read the transition wells for that quantity.
763 151 : if (calc_transition_bias_) {
764 13 : std::vector<double> tempcoords(getNumberOfArguments());
765 26 : for (unsigned i = 0; ; i++) {
766 78 : if (!parseNumberedVector("TRANSITIONWELL", i, tempcoords) ) {
767 : break;
768 : }
769 26 : if (tempcoords.size() != getNumberOfArguments()) {
770 0 : error("incorrect number of coordinates for transition tempering well");
771 : }
772 26 : transitionwells_.push_back(tempcoords);
773 : }
774 : }
775 :
776 302 : parse("TARGET",targetfilename_);
777 151 : if(targetfilename_.length()>0 && kbt_==0.0) {
778 0 : error("with TARGET temperature must be specified");
779 : }
780 151 : double tau=0.0;
781 151 : parse("TAU",tau);
782 151 : if(tau==0.0) {
783 129 : if(height0_==std::numeric_limits<double>::max()) {
784 0 : error("At least one between HEIGHT and TAU should be specified");
785 : }
786 : // if tau is not set, we compute it here from the other input parameters
787 129 : if(welltemp_) {
788 19 : tau=(kbt_*(biasf_-1.0))/height0_*getTimeStep()*stride_;
789 110 : } else if(dampfactor_>0.0) {
790 0 : tau=(kbt_*dampfactor_)/height0_*getTimeStep()*stride_;
791 : }
792 : } else {
793 22 : if(height0_!=std::numeric_limits<double>::max()) {
794 0 : error("At most one between HEIGHT and TAU should be specified");
795 : }
796 22 : if(welltemp_) {
797 19 : if(biasf_!=1.0) {
798 15 : height0_=(kbt_*(biasf_-1.0))/tau*getTimeStep()*stride_;
799 : } else {
800 4 : height0_=kbt_/tau*getTimeStep()*stride_; // special case for gamma=1
801 : }
802 3 : } else if(dampfactor_>0.0) {
803 2 : height0_=(kbt_*dampfactor_)/tau*getTimeStep()*stride_;
804 : } else {
805 1 : error("TAU only makes sense in well-tempered or damped metadynamics");
806 : }
807 : }
808 :
809 : // Grid Stuff
810 153 : std::vector<std::string> gmin(getNumberOfArguments());
811 300 : parseVector("GRID_MIN",gmin);
812 150 : if(gmin.size()!=getNumberOfArguments() && gmin.size()!=0) {
813 0 : error("not enough values for GRID_MIN");
814 : }
815 150 : std::vector<std::string> gmax(getNumberOfArguments());
816 300 : parseVector("GRID_MAX",gmax);
817 150 : if(gmax.size()!=getNumberOfArguments() && gmax.size()!=0) {
818 0 : error("not enough values for GRID_MAX");
819 : }
820 150 : std::vector<unsigned> gbin(getNumberOfArguments());
821 : std::vector<double> gspacing;
822 300 : parseVector("GRID_BIN",gbin);
823 150 : if(gbin.size()!=getNumberOfArguments() && gbin.size()!=0) {
824 0 : error("not enough values for GRID_BIN");
825 : }
826 300 : parseVector("GRID_SPACING",gspacing);
827 150 : if(gspacing.size()!=getNumberOfArguments() && gspacing.size()!=0) {
828 0 : error("not enough values for GRID_SPACING");
829 : }
830 150 : if(gmin.size()!=gmax.size()) {
831 0 : error("GRID_MAX and GRID_MIN should be either present or absent");
832 : }
833 150 : if(gspacing.size()!=0 && gmin.size()==0) {
834 0 : error("If GRID_SPACING is present also GRID_MIN and GRID_MAX should be present");
835 : }
836 150 : if(gbin.size()!=0 && gmin.size()==0) {
837 0 : error("If GRID_BIN is present also GRID_MIN and GRID_MAX should be present");
838 : }
839 150 : if(gmin.size()!=0) {
840 61 : if(gbin.size()==0 && gspacing.size()==0) {
841 6 : if(adaptive_==FlexibleBin::none) {
842 6 : log<<" Binsize not specified, 1/5 of sigma will be be used\n";
843 6 : plumed_assert(sigma0_.size()==getNumberOfArguments());
844 6 : gspacing.resize(getNumberOfArguments());
845 13 : for(unsigned i=0; i<gspacing.size(); i++) {
846 7 : gspacing[i]=0.2*sigma0_[i];
847 : }
848 : } else {
849 : // with adaptive hills and grid a sigma min must be specified
850 0 : for(unsigned i=0; i<sigma0min_.size(); i++)
851 0 : if(sigma0min_[i]<=0) {
852 0 : error("When using ADAPTIVE Gaussians on a grid SIGMA_MIN must be specified");
853 : }
854 0 : log<<" Binsize not specified, 1/5 of sigma_min will be be used\n";
855 0 : gspacing.resize(getNumberOfArguments());
856 0 : for(unsigned i=0; i<gspacing.size(); i++) {
857 0 : gspacing[i]=0.2*sigma0min_[i];
858 : }
859 : }
860 55 : } else if(gspacing.size()!=0 && gbin.size()==0) {
861 2 : log<<" The number of bins will be estimated from GRID_SPACING\n";
862 53 : } else if(gspacing.size()!=0 && gbin.size()!=0) {
863 1 : log<<" You specified both GRID_BIN and GRID_SPACING\n";
864 1 : log<<" The more conservative (highest) number of bins will be used for each variable\n";
865 : }
866 61 : if(gbin.size()==0) {
867 8 : gbin.assign(getNumberOfArguments(),1);
868 : }
869 61 : if(gspacing.size()!=0)
870 21 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
871 : double a,b;
872 13 : Tools::convert(gmin[i],a);
873 12 : Tools::convert(gmax[i],b);
874 12 : unsigned n=std::ceil(((b-a)/gspacing[i]));
875 12 : if(gbin[i]<n) {
876 11 : gbin[i]=n;
877 : }
878 : }
879 : }
880 149 : if(gbin.size()>0) {
881 60 : grid_=true;
882 : }
883 :
884 149 : bool sparsegrid=false;
885 149 : parseFlag("GRID_SPARSE",sparsegrid);
886 149 : bool nospline=false;
887 149 : parseFlag("GRID_NOSPLINE",nospline);
888 149 : bool spline=!nospline;
889 300 : parse("GRID_WSTRIDE",wgridstride_);
890 : std::string gridfilename_;
891 149 : parse("GRID_WFILE",gridfilename_);
892 149 : parseFlag("STORE_GRIDS",storeOldGrids_);
893 149 : if(grid_ && gridfilename_.length()>0) {
894 19 : if(wgridstride_==0 ) {
895 0 : error("frequency with which to output grid not specified use GRID_WSTRIDE");
896 : }
897 : }
898 149 : if(grid_ && wgridstride_>0) {
899 19 : if(gridfilename_.length()==0) {
900 1 : error("grid filename not specified use GRID_WFILE");
901 : }
902 : }
903 :
904 : std::string gridreadfilename_;
905 149 : parse("GRID_RFILE",gridreadfilename_);
906 :
907 149 : if(!grid_&&gridfilename_.length()> 0) {
908 0 : error("To write a grid you need first to define it!");
909 : }
910 149 : if(!grid_&&gridreadfilename_.length()>0) {
911 0 : error("To read a grid you need first to define it!");
912 : }
913 :
914 : /*setup neighbor list stuff*/
915 298 : parseFlag("NLIST", nlist_);
916 149 : nlist_center_.resize(getNumberOfArguments());
917 149 : nlist_dev2_.resize(getNumberOfArguments());
918 149 : if(nlist_&&grid_) {
919 1 : error("NLIST and GRID cannot be combined!");
920 : }
921 : std::vector<double> nlist_param;
922 298 : parseVector("NLIST_PARAMETERS",nlist_param);
923 149 : if(nlist_param.size()==0) {
924 149 : nlist_param_[0]=6.0;//*DP2CUTOFF -> max distance of neighbors
925 149 : nlist_param_[1]=0.5;//*nlist_dev2_[i] -> condition for rebuilding
926 : } else {
927 0 : plumed_massert(nlist_param.size()==2,"two cutoff parameters are needed for the neighbor list");
928 0 : plumed_massert(nlist_param[0]>1.0,"the first of NLIST_PARAMETERS must be greater than 1. The smaller the first, the smaller should be the second as well");
929 0 : const double min_PARAM_1=(1.-1./std::sqrt(nlist_param[0]/2))+0.16;
930 0 : plumed_massert(nlist_param[1]>0,"the second of NLIST_PARAMETERS must be greater than 0");
931 0 : plumed_massert(nlist_param[1]<=min_PARAM_1,"the second of NLIST_PARAMETERS must be smaller to avoid systematic errors. Largest suggested value is: 1.16-1/sqrt(PARAM_0/2) = "+std::to_string(min_PARAM_1));
932 0 : nlist_param_[0]=nlist_param[0];
933 0 : nlist_param_[1]=nlist_param[1];
934 : }
935 :
936 : // Reweighting factor rct
937 149 : parseFlag("CALC_RCT",calc_rct_);
938 149 : if (calc_rct_) {
939 6 : plumed_massert(grid_,"CALC_RCT is supported only if bias is on a grid");
940 : }
941 149 : parse("RCT_USTRIDE",rct_ustride_);
942 :
943 149 : if(dampfactor_>0.0) {
944 2 : if(!grid_) {
945 0 : error("With DAMPFACTOR you should use grids");
946 : }
947 : }
948 :
949 : // Multiple walkers
950 149 : parse("WALKERS_N",mw_n_);
951 149 : parse("WALKERS_ID",mw_id_);
952 149 : if(mw_n_<=mw_id_) {
953 0 : error("walker ID should be a numerical value less than the total number of walkers");
954 : }
955 149 : parse("WALKERS_DIR",mw_dir_);
956 149 : parse("WALKERS_RSTRIDE",mw_rstride_);
957 :
958 : // MPI version
959 149 : parseFlag("WALKERS_MPI",walkers_mpi_);
960 :
961 : //If this Action is not compiled with MPI the user is informed and we exit gracefully
962 149 : if(walkers_mpi_) {
963 39 : plumed_assert(Communicator::plumedHasMPI()) << "Invalid walkers configuration: WALKERS_MPI flag requires MPI compilation";
964 40 : plumed_assert(Communicator::initialized()) << "Invalid walkers configuration: WALKERS_MPI needs the communicator correctly initialized.";
965 : }
966 :
967 : // Flying Gaussian
968 148 : parseFlag("FLYING_GAUSSIAN", flying_);
969 :
970 : // Inteval keyword
971 149 : std::vector<double> tmpI(2);
972 296 : parseVector("INTERVAL",tmpI);
973 148 : if(tmpI.size()!=2&&tmpI.size()!=0) {
974 0 : error("both a lower and an upper limits must be provided with INTERVAL");
975 148 : } else if(tmpI.size()==2) {
976 2 : lowI_=tmpI.at(0);
977 2 : uppI_=tmpI.at(1);
978 2 : if(getNumberOfArguments()!=1) {
979 0 : error("INTERVAL limits correction works only for monodimensional metadynamics!");
980 : }
981 2 : if(uppI_<lowI_) {
982 0 : error("The Upper limit must be greater than the Lower limit!");
983 : }
984 2 : if(getPntrToArgument(0)->isPeriodic()) {
985 0 : error("INTERVAL cannot be used with periodic variables!");
986 : }
987 2 : doInt_=true;
988 : }
989 :
990 296 : parseFlag("ACCELERATION",acceleration_);
991 : // Check for a restart acceleration if acceleration is active.
992 : std::string acc_rfilename;
993 148 : if(acceleration_) {
994 8 : parse("ACCELERATION_RFILE", acc_rfilename);
995 : }
996 :
997 148 : freq_adaptive_=false;
998 148 : parseFlag("FREQUENCY_ADAPTIVE",freq_adaptive_);
999 : //
1000 148 : fa_update_frequency_=0;
1001 148 : parse("FA_UPDATE_FREQUENCY",fa_update_frequency_);
1002 148 : if(fa_update_frequency_!=0 && !freq_adaptive_) {
1003 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");
1004 : }
1005 148 : if(fa_update_frequency_==0 && freq_adaptive_) {
1006 0 : fa_update_frequency_=stride_;
1007 : }
1008 : //
1009 148 : fa_max_stride_=0;
1010 148 : parse("FA_MAX_PACE",fa_max_stride_);
1011 148 : if(fa_max_stride_!=0 && !freq_adaptive_) {
1012 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");
1013 : }
1014 : //
1015 148 : fa_min_acceleration_=1.0;
1016 148 : parse("FA_MIN_ACCELERATION",fa_min_acceleration_);
1017 148 : if(fa_min_acceleration_!=1.0 && !freq_adaptive_) {
1018 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");
1019 : }
1020 :
1021 148 : checkRead();
1022 :
1023 148 : log.printf(" Gaussian width ");
1024 148 : if (adaptive_==FlexibleBin::diffusion) {
1025 2 : log.printf(" (Note: The units of sigma are in timesteps) ");
1026 : }
1027 148 : if (adaptive_==FlexibleBin::geometry) {
1028 19 : log.printf(" (Note: The units of sigma are in dist units) ");
1029 : }
1030 396 : for(unsigned i=0; i<sigma0_.size(); ++i) {
1031 248 : log.printf(" %f",sigma0_[i]);
1032 : }
1033 148 : log.printf(" Gaussian height %f\n",height0_);
1034 148 : log.printf(" Gaussian deposition pace %d\n",stride_);
1035 148 : log.printf(" Gaussian file %s\n",hillsfname.c_str());
1036 148 : if(welltemp_) {
1037 38 : log.printf(" Well-Tempered Bias Factor %f\n",biasf_);
1038 38 : log.printf(" Hills relaxation time (tau) %f\n",tau);
1039 38 : log.printf(" KbT %f\n",kbt_);
1040 : }
1041 :
1042 : // Transition tempered metadynamics options
1043 148 : if (tt_specs_.is_active) {
1044 3 : logTemperingSpecs(tt_specs_);
1045 : // Check that the appropriate transition bias quantity is calculated.
1046 : // (Should never trip, given that the flag is automatically set.)
1047 3 : if (!calc_transition_bias_) {
1048 0 : error(" transition tempering requires calculation of a transition bias");
1049 : }
1050 : }
1051 :
1052 : // Overall tempering sanity check (this gets tricky when multiple are active).
1053 : // When multiple temperings are active, it's fine to have one tempering attempt
1054 : // to increase hill size with increasing bias, so long as the others can shrink
1055 : // the hills faster than it increases their size in the long-time limit.
1056 : // This set of checks ensures that the hill sizes eventually decay to zero as c(t)
1057 : // diverges to infinity.
1058 : // The alpha parameter allows hills to decay as 1/t^alpha instead of 1/t,
1059 : // a slower decay, so as t -> infinity, only the temperings with the largest
1060 : // alphas govern the final asymptotic decay. (Alpha helps prevent false convergence.)
1061 148 : if (welltemp_ || dampfactor_ > 0.0 || tt_specs_.is_active) {
1062 : // Determine the number of active temperings.
1063 : int n_active = 0;
1064 43 : if (welltemp_) {
1065 : n_active++;
1066 : }
1067 43 : if (dampfactor_ > 0.0) {
1068 2 : n_active++;
1069 : }
1070 43 : if (tt_specs_.is_active) {
1071 3 : n_active++;
1072 : }
1073 : // Find the greatest alpha.
1074 43 : double greatest_alpha = 0.0;
1075 43 : if (welltemp_) {
1076 38 : greatest_alpha = std::max(greatest_alpha, 1.0);
1077 : }
1078 43 : if (dampfactor_ > 0.0) {
1079 4 : greatest_alpha = std::max(greatest_alpha, 1.0);
1080 : }
1081 43 : if (tt_specs_.is_active) {
1082 6 : greatest_alpha = std::max(greatest_alpha, tt_specs_.alpha);
1083 : }
1084 : // Find the least alpha.
1085 43 : double least_alpha = 1.0;
1086 : if (welltemp_) {
1087 : least_alpha = std::min(least_alpha, 1.0);
1088 : }
1089 43 : if (dampfactor_ > 0.0) {
1090 2 : least_alpha = std::min(least_alpha, 1.0);
1091 : }
1092 43 : if (tt_specs_.is_active) {
1093 4 : least_alpha = std::min(least_alpha, tt_specs_.alpha);
1094 : }
1095 : // Find the inverse harmonic average of the delta T parameters for all
1096 : // of the temperings with the greatest alpha values.
1097 : double total_governing_deltaT_inv = 0.0;
1098 43 : if (welltemp_ && 1.0 == greatest_alpha && biasf_ != 1.0) {
1099 34 : total_governing_deltaT_inv += 1.0 / (biasf_ - 1.0);
1100 : }
1101 43 : if (dampfactor_ > 0.0 && 1.0 == greatest_alpha) {
1102 2 : total_governing_deltaT_inv += 1.0 / (dampfactor_);
1103 : }
1104 43 : if (tt_specs_.is_active && tt_specs_.alpha == greatest_alpha) {
1105 3 : total_governing_deltaT_inv += 1.0 / (tt_specs_.biasf - 1.0);
1106 : }
1107 : // Give a newbie-friendly error message for people using one tempering if
1108 : // only one is active.
1109 43 : if (n_active == 1 && total_governing_deltaT_inv < 0.0) {
1110 0 : error("for stable tempering, the bias factor must be greater than one");
1111 : // Give a slightly more complex error message to users stacking multiple
1112 : // tempering options at a time, but all with uniform alpha values.
1113 43 : } else if (total_governing_deltaT_inv < 0.0 && greatest_alpha == least_alpha) {
1114 0 : error("for stable tempering, the sum of the inverse Delta T parameters must be greater than zero!");
1115 : // Give the most technical error message to users stacking multiple tempering
1116 : // options with different alpha parameters.
1117 43 : } else if (total_governing_deltaT_inv < 0.0 && greatest_alpha != least_alpha) {
1118 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!");
1119 : }
1120 : }
1121 :
1122 148 : if(doInt_) {
1123 2 : log.printf(" Upper and Lower limits boundaries for the bias are activated at %f - %f\n", lowI_, uppI_);
1124 : }
1125 :
1126 148 : if(grid_) {
1127 60 : log.printf(" Grid min");
1128 161 : for(unsigned i=0; i<gmin.size(); ++i) {
1129 101 : log.printf(" %s",gmin[i].c_str() );
1130 : }
1131 60 : log.printf("\n");
1132 60 : log.printf(" Grid max");
1133 161 : for(unsigned i=0; i<gmax.size(); ++i) {
1134 101 : log.printf(" %s",gmax[i].c_str() );
1135 : }
1136 60 : log.printf("\n");
1137 60 : log.printf(" Grid bin");
1138 161 : for(unsigned i=0; i<gbin.size(); ++i) {
1139 101 : log.printf(" %u",gbin[i]);
1140 : }
1141 60 : log.printf("\n");
1142 60 : if(spline) {
1143 60 : log.printf(" Grid uses spline interpolation\n");
1144 : }
1145 60 : if(sparsegrid) {
1146 6 : log.printf(" Grid uses sparse grid\n");
1147 : }
1148 60 : if(wgridstride_>0) {
1149 19 : log.printf(" Grid is written on file %s with stride %d\n",gridfilename_.c_str(),wgridstride_);
1150 : }
1151 : }
1152 :
1153 148 : if(mw_n_>1) {
1154 6 : if(walkers_mpi_) {
1155 0 : error("MPI version of multiple walkers is not compatible with filesystem version of multiple walkers");
1156 : }
1157 6 : log.printf(" %d multiple walkers active\n",mw_n_);
1158 6 : log.printf(" walker id %d\n",mw_id_);
1159 6 : log.printf(" reading stride %d\n",mw_rstride_);
1160 6 : if(mw_dir_!="") {
1161 3 : log.printf(" directory with hills files %s\n",mw_dir_.c_str());
1162 : }
1163 : } else {
1164 142 : if(walkers_mpi_) {
1165 38 : log.printf(" Multiple walkers active using MPI communnication\n");
1166 38 : if(mw_dir_!="") {
1167 0 : log.printf(" directory with hills files %s\n",mw_dir_.c_str());
1168 : }
1169 38 : if(comm.Get_rank()==0) {
1170 : // Only root of group can communicate with other walkers
1171 23 : mpi_nw_=multi_sim_comm.Get_size();
1172 : }
1173 : // Communicate to the other members of the same group
1174 : // info abount number of walkers and walker index
1175 38 : comm.Bcast(mpi_nw_,0);
1176 : }
1177 : }
1178 :
1179 148 : if(flying_) {
1180 6 : if(!walkers_mpi_) {
1181 0 : error("Flying Gaussian method must be used with MPI version of multiple walkers");
1182 : }
1183 6 : log.printf(" Flying Gaussian method with %d walkers active\n",mpi_nw_);
1184 : }
1185 :
1186 148 : if(nlist_) {
1187 1 : addComponent("nlker");
1188 1 : componentIsNotPeriodic("nlker");
1189 1 : addComponent("nlsteps");
1190 2 : componentIsNotPeriodic("nlsteps");
1191 : }
1192 :
1193 148 : if(calc_rct_) {
1194 6 : addComponent("rbias");
1195 6 : componentIsNotPeriodic("rbias");
1196 6 : addComponent("rct");
1197 6 : componentIsNotPeriodic("rct");
1198 6 : log.printf(" The c(t) reweighting factor will be calculated every %u hills\n",rct_ustride_);
1199 12 : getPntrToComponent("rct")->set(reweight_factor_);
1200 : }
1201 :
1202 148 : if(calc_work_) {
1203 1 : addComponent("work");
1204 2 : componentIsNotPeriodic("work");
1205 : }
1206 :
1207 148 : if(acceleration_) {
1208 4 : if (kbt_ == 0.0) {
1209 0 : error("The calculation of the acceleration works only if simulation temperature has been defined");
1210 : }
1211 4 : log.printf(" calculation on the fly of the acceleration factor\n");
1212 4 : addComponent("acc");
1213 8 : componentIsNotPeriodic("acc");
1214 : // Set the initial value of the the acceleration.
1215 : // If this is not a restart, set to 1.0.
1216 4 : if (acc_rfilename.length() == 0) {
1217 4 : getPntrToComponent("acc")->set(1.0);
1218 2 : if(getRestart()) {
1219 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.\n");
1220 1 : log.printf(" You should use the ACCELERATION_RFILE keyword.\n");
1221 : }
1222 : // Otherwise, read and set the restart value.
1223 : } else {
1224 : // Restart of acceleration does not make sense if the restart timestep is zero.
1225 : //if (getStep() == 0) {
1226 : // error("Restarting calculation of acceleration factors works only if simulation timestep is restarted correctly");
1227 : //}
1228 : // Open the ACCELERATION_RFILE.
1229 2 : IFile acc_rfile;
1230 2 : acc_rfile.link(*this);
1231 2 : if(acc_rfile.FileExist(acc_rfilename)) {
1232 2 : acc_rfile.open(acc_rfilename);
1233 : } else {
1234 0 : error("The ACCELERATION_RFILE file you want to read: " + acc_rfilename + ", cannot be found!");
1235 : }
1236 : // Read the file to find the restart acceleration.
1237 2 : double acc_rmean=0.0;
1238 2 : double acc_rtime=0.0;
1239 : bool found=false;
1240 2 : std::string acclabel = getLabel() + ".acc";
1241 2 : acc_rfile.allowIgnoredFields();
1242 248 : while(acc_rfile.scanField("time", acc_rtime)) {
1243 122 : acc_rfile.scanField(acclabel, acc_rmean);
1244 122 : acc_rfile.scanField();
1245 : found=true;
1246 : }
1247 2 : if(!found) {
1248 0 : error("The ACCELERATION_RFILE file you want to read: " + acc_rfilename + ", does not contain a time field!");
1249 : }
1250 2 : acc_restart_mean_ = acc_rmean;
1251 : // Set component based on the read values.
1252 2 : getPntrToComponent("acc")->set(acc_rmean);
1253 2 : log.printf(" initial acceleration factor read from file %s: value of %f at time %f\n",acc_rfilename.c_str(),acc_rmean,acc_rtime);
1254 2 : }
1255 : }
1256 :
1257 148 : if (calc_max_bias_) {
1258 0 : if (!grid_) {
1259 0 : error("Calculating the maximum bias on the fly works only with a grid");
1260 : }
1261 0 : log.printf(" calculation on the fly of the maximum bias max(V(s,t)) \n");
1262 0 : addComponent("maxbias");
1263 0 : componentIsNotPeriodic("maxbias");
1264 : }
1265 :
1266 148 : if (calc_transition_bias_) {
1267 13 : if (!grid_) {
1268 0 : error("Calculating the transition bias on the fly works only with a grid");
1269 : }
1270 13 : log.printf(" calculation on the fly of the transition bias V*(t)\n");
1271 13 : addComponent("transbias");
1272 13 : componentIsNotPeriodic("transbias");
1273 13 : log<<" Number of transition wells "<<transitionwells_.size()<<"\n";
1274 13 : if (transitionwells_.size() == 0) {
1275 0 : error("Calculating the transition bias on the fly requires definition of at least one transition well");
1276 : }
1277 : // Check that a grid is in use.
1278 13 : if (!grid_) {
1279 0 : error(" transition barrier finding requires a grid for the bias");
1280 : }
1281 : // Log the wells and check that they are in the grid.
1282 39 : for (unsigned i = 0; i < transitionwells_.size(); i++) {
1283 : // Log the coordinate.
1284 26 : log.printf(" Transition well %d at coordinate ", i);
1285 64 : for (unsigned j = 0; j < getNumberOfArguments(); j++) {
1286 38 : log.printf("%f ", transitionwells_[i][j]);
1287 : }
1288 26 : log.printf("\n");
1289 : // Check that the coordinate is in the grid.
1290 64 : for (unsigned j = 0; j < getNumberOfArguments(); j++) {
1291 : double max, min;
1292 38 : Tools::convert(gmin[j], min);
1293 38 : Tools::convert(gmax[j], max);
1294 38 : if (transitionwells_[i][j] < min || transitionwells_[i][j] > max) {
1295 0 : error(" transition well is not in grid");
1296 : }
1297 : }
1298 : }
1299 : }
1300 :
1301 148 : if(freq_adaptive_) {
1302 2 : if(!acceleration_) {
1303 0 : plumed_merror("Frequency adaptive metadynamics only works if the calculation of the acceleration factor is enabled with the ACCELERATION keyword\n");
1304 : }
1305 2 : if(walkers_mpi_) {
1306 0 : plumed_merror("Combining frequency adaptive metadynamics with MPI multiple walkers is not allowed");
1307 : }
1308 :
1309 2 : log.printf(" Frequency adaptive metadynamics enabled\n");
1310 2 : if(getRestart() && acc_rfilename.length() == 0) {
1311 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.\n");
1312 0 : log.printf(" You should use the ACCELERATION_RFILE keyword.\n");
1313 : }
1314 2 : log.printf(" The frequency for hill addition will change dynamically based on the metadynamics acceleration factor\n");
1315 2 : log.printf(" The hill addition frequency will be updated every %d steps\n",fa_update_frequency_);
1316 2 : if(fa_min_acceleration_>1.0) {
1317 2 : log.printf(" The hill addition frequency will only be updated once the metadynamics acceleration factor becomes larger than %.1f \n",fa_min_acceleration_);
1318 : }
1319 2 : if(fa_max_stride_!=0) {
1320 2 : log.printf(" The hill addition frequency will not become larger than %d steps\n",fa_max_stride_);
1321 : }
1322 2 : addComponent("pace");
1323 2 : componentIsNotPeriodic("pace");
1324 2 : updateFrequencyAdaptiveStride();
1325 : }
1326 :
1327 : // initializing and checking grid
1328 : bool restartedFromGrid=false; // restart from grid file
1329 148 : if(grid_) {
1330 60 : if(!(gridreadfilename_.length()>0)) {
1331 : // check for mesh and sigma size
1332 116 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
1333 : double a,b;
1334 74 : Tools::convert(gmin[i],a);
1335 74 : Tools::convert(gmax[i],b);
1336 74 : double mesh=(b-a)/((double)gbin[i]);
1337 74 : if(adaptive_==FlexibleBin::none) {
1338 74 : if(mesh>0.5*sigma0_[i]) {
1339 38 : log<<" WARNING: Using a METAD with a Grid Spacing larger than half of the Gaussians width (SIGMA) can produce artifacts\n";
1340 : }
1341 : } else {
1342 0 : if(sigma0min_[i]<0.) {
1343 0 : error("When using ADAPTIVE Gaussians on a grid SIGMA_MIN must be specified");
1344 : }
1345 0 : if(mesh>0.5*sigma0min_[i]) {
1346 0 : log<<" WARNING: to use a METAD with a GRID and ADAPTIVE you need to set a Grid Spacing lower than half of the Gaussians (SIGMA_MIN) \n";
1347 : }
1348 : }
1349 : }
1350 42 : std::string funcl=getLabel() + ".bias";
1351 42 : if(!sparsegrid) {
1352 36 : BiasGrid_=Tools::make_unique<Grid>(funcl,getArguments(),gmin,gmax,gbin,spline,true);
1353 : } else {
1354 6 : BiasGrid_=Tools::make_unique<SparseGrid>(funcl,getArguments(),gmin,gmax,gbin,spline,true);
1355 : }
1356 42 : std::vector<std::string> actualmin=BiasGrid_->getMin();
1357 42 : std::vector<std::string> actualmax=BiasGrid_->getMax();
1358 116 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
1359 : std::string is;
1360 74 : Tools::convert(i,is);
1361 74 : if(gmin[i]!=actualmin[i]) {
1362 0 : error("GRID_MIN["+is+"] must be adjusted to "+actualmin[i]+" to fit periodicity");
1363 : }
1364 74 : if(gmax[i]!=actualmax[i]) {
1365 0 : error("GRID_MAX["+is+"] must be adjusted to "+actualmax[i]+" to fit periodicity");
1366 : }
1367 : }
1368 42 : } else {
1369 : // read the grid in input, find the keys
1370 : #ifdef __PLUMED_HAS_GETCWD
1371 18 : if(walkers_mpi_&&gridreadfilename_.at(0)!='/') {
1372 : //if possible the root replica will share its current folder so that all walkers will read the same file
1373 0 : char cwd[4096]= {0};
1374 : const char* ret=getcwd(cwd,4096);
1375 0 : plumed_assert(ret)<<"Name of current directory too long, increase buffer size";
1376 0 : gridreadfilename_ = "/" + gridreadfilename_;
1377 0 : gridreadfilename_ = ret + gridreadfilename_;
1378 0 : if(comm.Get_rank()==0) {
1379 0 : multi_sim_comm.Bcast(gridreadfilename_,0);
1380 : }
1381 0 : comm.Bcast(gridreadfilename_,0);
1382 : }
1383 : #endif
1384 18 : IFile gridfile;
1385 18 : gridfile.link(*this);
1386 18 : if(gridfile.FileExist(gridreadfilename_)) {
1387 18 : gridfile.open(gridreadfilename_);
1388 : } else {
1389 0 : error("The GRID file you want to read: " + gridreadfilename_ + ", cannot be found!");
1390 : }
1391 18 : std::string funcl=getLabel() + ".bias";
1392 36 : BiasGrid_=GridBase::create(funcl, getArguments(), gridfile, gmin, gmax, gbin, sparsegrid, spline, true);
1393 18 : if(BiasGrid_->getDimension()!=getNumberOfArguments()) {
1394 0 : error("mismatch between dimensionality of input grid and number of arguments");
1395 : }
1396 45 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1397 54 : if( getPntrToArgument(i)->isPeriodic()!=BiasGrid_->getIsPeriodic()[i] ) {
1398 0 : error("periodicity mismatch between arguments and input bias");
1399 : }
1400 : double a, b;
1401 27 : Tools::convert(gmin[i],a);
1402 27 : Tools::convert(gmax[i],b);
1403 27 : double mesh=(b-a)/((double)gbin[i]);
1404 27 : if(mesh>0.5*sigma0_[i]) {
1405 27 : log<<" WARNING: Using a METAD with a Grid Spacing larger than half of the Gaussians width can produce artifacts\n";
1406 : }
1407 : }
1408 18 : log.printf(" Restarting from %s\n",gridreadfilename_.c_str());
1409 18 : if(getRestart()) {
1410 : restartedFromGrid=true;
1411 : }
1412 18 : }
1413 : }
1414 :
1415 : // if we are restarting from GRID and using WALKERS_MPI we can check that all walkers have actually read the grid
1416 148 : if(getRestart()&&walkers_mpi_) {
1417 9 : std::vector<int> restarted(mpi_nw_,0);
1418 9 : if(comm.Get_rank()==0) {
1419 6 : multi_sim_comm.Allgather(int(restartedFromGrid), restarted);
1420 : }
1421 9 : comm.Bcast(restarted,0);
1422 : int result = std::accumulate(restarted.begin(),restarted.end(),0);
1423 9 : if(result!=0&&result!=mpi_nw_) {
1424 0 : error("in this WALKERS_MPI run some replica have restarted from GRID while other do not!");
1425 : }
1426 : }
1427 :
1428 : #ifdef __PLUMED_HAS_GETCWD
1429 186 : if(walkers_mpi_&&mw_dir_==""&&hillsfname.at(0)!='/') {
1430 : //if possible the root replica will share its current folder so that all walkers will read the same file
1431 38 : char cwd[4096]= {0};
1432 : const char* ret=getcwd(cwd,4096);
1433 38 : plumed_assert(ret)<<"Name of current directory too long, increase buffer size";
1434 : mw_dir_ = ret;
1435 38 : mw_dir_ = mw_dir_ + "/";
1436 38 : if(comm.Get_rank()==0) {
1437 23 : multi_sim_comm.Bcast(mw_dir_,0);
1438 : }
1439 38 : comm.Bcast(mw_dir_,0);
1440 : }
1441 : #endif
1442 :
1443 : // creating std::vector of ifile* for hills reading
1444 : // open all files at the beginning and read Gaussians if restarting
1445 : bool restartedFromHills=false; // restart from hills files
1446 308 : for(int i=0; i<mw_n_; ++i) {
1447 : std::string fname;
1448 160 : if(mw_dir_!="") {
1449 47 : if(mw_n_>1) {
1450 9 : std::stringstream out;
1451 9 : out << i;
1452 18 : fname = mw_dir_+"/"+hillsfname+"."+out.str();
1453 47 : } else if(walkers_mpi_) {
1454 76 : fname = mw_dir_+"/"+hillsfname;
1455 : } else {
1456 : fname = hillsfname;
1457 : }
1458 : } else {
1459 113 : if(mw_n_>1) {
1460 9 : std::stringstream out;
1461 9 : out << i;
1462 18 : fname = hillsfname+"."+out.str();
1463 9 : } else {
1464 : fname = hillsfname;
1465 : }
1466 : }
1467 160 : ifiles_.emplace_back(Tools::make_unique<IFile>());
1468 : // this is just a shortcut pointer to the last element:
1469 : IFile *ifile = ifiles_.back().get();
1470 160 : ifilesnames_.push_back(fname);
1471 160 : ifile->link(*this);
1472 160 : if(ifile->FileExist(fname)) {
1473 35 : ifile->open(fname);
1474 35 : if(getRestart()&&!restartedFromGrid) {
1475 19 : log.printf(" Restarting from %s:",ifilesnames_[i].c_str());
1476 19 : readGaussians(ifiles_[i].get());
1477 : restartedFromHills=true;
1478 : }
1479 35 : ifiles_[i]->reset(false);
1480 : // close only the walker own hills file for later writing
1481 35 : if(i==mw_id_) {
1482 30 : ifiles_[i]->close();
1483 : }
1484 : } else {
1485 : // in case a file does not exist and we are restarting, complain that the file was not found
1486 125 : if(getRestart()&&!restartedFromGrid) {
1487 0 : error("restart file "+fname+" not found");
1488 : }
1489 : }
1490 : }
1491 :
1492 : // if we are restarting from FILE and using WALKERS_MPI we can check that all walkers have actually read the FILE
1493 148 : if(getRestart()&&walkers_mpi_) {
1494 9 : std::vector<int> restarted(mpi_nw_,0);
1495 9 : if(comm.Get_rank()==0) {
1496 6 : multi_sim_comm.Allgather(int(restartedFromHills), restarted);
1497 : }
1498 9 : comm.Bcast(restarted,0);
1499 : int result = std::accumulate(restarted.begin(),restarted.end(),0);
1500 9 : if(result!=0&&result!=mpi_nw_) {
1501 0 : error("in this WALKERS_MPI run some replica have restarted from FILE while other do not!");
1502 : }
1503 : }
1504 :
1505 148 : comm.Barrier();
1506 : // this barrier is needed when using walkers_mpi
1507 : // to be sure that all files have been read before
1508 : // backing them up
1509 : // it should not be used when walkers_mpi is false otherwise
1510 : // it would introduce troubles when using replicas without METAD
1511 : // (e.g. in bias exchange with a neutral replica)
1512 : // see issue #168 on github
1513 148 : if(comm.Get_rank()==0 && walkers_mpi_) {
1514 23 : multi_sim_comm.Barrier();
1515 : }
1516 :
1517 148 : if(targetfilename_.length()>0) {
1518 2 : IFile gridfile;
1519 2 : gridfile.open(targetfilename_);
1520 2 : std::string funcl=getLabel() + ".target";
1521 4 : TargetGrid_=GridBase::create(funcl,getArguments(),gridfile,false,false,true);
1522 2 : if(TargetGrid_->getDimension()!=getNumberOfArguments()) {
1523 0 : error("mismatch between dimensionality of input grid and number of arguments");
1524 : }
1525 4 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1526 4 : if( getPntrToArgument(i)->isPeriodic()!=TargetGrid_->getIsPeriodic()[i] ) {
1527 0 : error("periodicity mismatch between arguments and input bias");
1528 : }
1529 : }
1530 2 : }
1531 :
1532 148 : if(getRestart()) {
1533 : // if this is a restart the neighbor list should be immediately updated
1534 37 : if(nlist_) {
1535 1 : nlist_update_=true;
1536 : }
1537 : // Calculate the Tiwary-Parrinello reweighting factor if we are restarting from previous hills
1538 37 : if(calc_rct_) {
1539 0 : computeReweightingFactor();
1540 : }
1541 : // Calculate all special bias quantities desired if restarting with nonzero bias.
1542 37 : if(calc_max_bias_) {
1543 0 : max_bias_ = BiasGrid_->getMaxValue();
1544 0 : getPntrToComponent("maxbias")->set(max_bias_);
1545 : }
1546 37 : if(calc_transition_bias_) {
1547 13 : transition_bias_ = getTransitionBarrierBias();
1548 26 : getPntrToComponent("transbias")->set(transition_bias_);
1549 : }
1550 : }
1551 :
1552 : // open grid file for writing
1553 148 : if(wgridstride_>0) {
1554 19 : gridfile_.link(*this);
1555 19 : if(walkers_mpi_) {
1556 0 : int r=0;
1557 0 : if(comm.Get_rank()==0) {
1558 0 : r=multi_sim_comm.Get_rank();
1559 : }
1560 0 : comm.Bcast(r,0);
1561 0 : if(r>0) {
1562 : gridfilename_="/dev/null";
1563 : }
1564 0 : gridfile_.enforceSuffix("");
1565 : }
1566 19 : if(mw_n_>1) {
1567 0 : gridfile_.enforceSuffix("");
1568 : }
1569 19 : gridfile_.open(gridfilename_);
1570 : }
1571 :
1572 : // open hills file for writing
1573 148 : hillsOfile_.link(*this);
1574 148 : if(walkers_mpi_) {
1575 38 : int r=0;
1576 38 : if(comm.Get_rank()==0) {
1577 23 : r=multi_sim_comm.Get_rank();
1578 : }
1579 38 : comm.Bcast(r,0);
1580 38 : if(r>0) {
1581 25 : ifilesnames_[mw_id_]="/dev/null";
1582 : }
1583 76 : hillsOfile_.enforceSuffix("");
1584 : }
1585 148 : if(mw_n_>1) {
1586 12 : hillsOfile_.enforceSuffix("");
1587 : }
1588 148 : hillsOfile_.open(ifilesnames_[mw_id_]);
1589 148 : if(fmt_.length()>0) {
1590 117 : hillsOfile_.fmtField(fmt_);
1591 : }
1592 148 : hillsOfile_.addConstantField("multivariate");
1593 148 : hillsOfile_.addConstantField("kerneltype");
1594 148 : if(doInt_) {
1595 4 : hillsOfile_.addConstantField("lower_int").printField("lower_int",lowI_);
1596 4 : hillsOfile_.addConstantField("upper_int").printField("upper_int",uppI_);
1597 : }
1598 : hillsOfile_.setHeavyFlush();
1599 : // output periodicities of variables
1600 415 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1601 267 : hillsOfile_.setupPrintValue( getPntrToArgument(i) );
1602 : }
1603 :
1604 : bool concurrent=false;
1605 148 : const ActionSet&actionSet(plumed.getActionSet());
1606 541 : for(const auto & p : actionSet)
1607 467 : if(dynamic_cast<MetaD*>(p.get())) {
1608 : concurrent=true;
1609 : break;
1610 : }
1611 148 : if(concurrent) {
1612 74 : log<<" You are using concurrent metadynamics\n";
1613 : }
1614 148 : if(rect_biasf_.size()>0) {
1615 18 : if(walkers_mpi_) {
1616 12 : log<<" You are using RECT in its 'altruistic' implementation\n";
1617 : }{
1618 18 : log<<" You are using RECT\n";
1619 : }
1620 : }
1621 :
1622 296 : log<<" Bibliography "<<plumed.cite("Laio and Parrinello, PNAS 99, 12562 (2002)");
1623 148 : if(welltemp_) {
1624 76 : log<<plumed.cite("Barducci, Bussi, and Parrinello, Phys. Rev. Lett. 100, 020603 (2008)");
1625 : }
1626 148 : if(tt_specs_.is_active) {
1627 6 : log << plumed.cite("Dama, Rotskoff, Parrinello, and Voth, J. Chem. Theory Comput. 10, 3626 (2014)");
1628 6 : log << plumed.cite("Dama, Parrinello, and Voth, Phys. Rev. Lett. 112, 240602 (2014)");
1629 : }
1630 148 : if(mw_n_>1||walkers_mpi_) {
1631 88 : log<<plumed.cite("Raiteri, Laio, Gervasio, Micheletti, and Parrinello, J. Phys. Chem. B 110, 3533 (2006)");
1632 : }
1633 148 : if(adaptive_!=FlexibleBin::none) {
1634 42 : log<<plumed.cite("Branduardi, Bussi, and Parrinello, J. Chem. Theory Comput. 8, 2247 (2012)");
1635 : }
1636 148 : if(doInt_) {
1637 4 : log<<plumed.cite("Baftizadeh, Cossio, Pietrucci, and Laio, Curr. Phys. Chem. 2, 79 (2012)");
1638 : }
1639 148 : if(acceleration_) {
1640 8 : log<<plumed.cite("Pratyush and Parrinello, Phys. Rev. Lett. 111, 230602 (2013)");
1641 : }
1642 148 : if(calc_rct_) {
1643 12 : log<<plumed.cite("Pratyush and Parrinello, J. Phys. Chem. B, 119, 736 (2015)");
1644 : }
1645 148 : if(concurrent || rect_biasf_.size()>0) {
1646 160 : log<<plumed.cite("Gil-Ley and Bussi, J. Chem. Theory Comput. 11, 1077 (2015)");
1647 : }
1648 148 : if(rect_biasf_.size()>0 && walkers_mpi_) {
1649 24 : log<<plumed.cite("Hosek, Toulcova, Bortolato, and Spiwok, J. Phys. Chem. B 120, 2209 (2016)");
1650 : }
1651 148 : if(targetfilename_.length()>0) {
1652 4 : log<<plumed.cite("White, Dama, and Voth, J. Chem. Theory Comput. 11, 2451 (2015)");
1653 4 : log<<plumed.cite("Marinelli and Faraldo-Gómez, Biophys. J. 108, 2779 (2015)");
1654 4 : log<<plumed.cite("Gil-Ley, Bottaro, and Bussi, J. Chem. Theory Comput. 12, 2790 (2016)");
1655 : }
1656 148 : if(freq_adaptive_) {
1657 4 : log<<plumed.cite("Wang, Valsson, Tiwary, Parrinello, and Lindorff-Larsen, J. Chem. Phys. 149, 072309 (2018)");
1658 : }
1659 148 : log<<"\n";
1660 396 : }
1661 :
1662 151 : void MetaD::readTemperingSpecs(TemperingSpecs &t_specs) {
1663 : // Set global tempering parameters.
1664 151 : parse(t_specs.name_stem + "BIASFACTOR", t_specs.biasf);
1665 151 : if (t_specs.biasf != -1.0) {
1666 3 : if (kbt_ == 0.0) {
1667 0 : error("Unless the MD engine passes the temperature to plumed, with tempered metad you must specify it using TEMP");
1668 : }
1669 3 : if (t_specs.biasf == 1.0) {
1670 0 : error("A bias factor of 1 corresponds to zero delta T and zero hill size, so it is not allowed.");
1671 : }
1672 3 : t_specs.is_active = true;
1673 3 : parse(t_specs.name_stem + "BIASTHRESHOLD", t_specs.threshold);
1674 3 : if (t_specs.threshold < 0.0) {
1675 0 : error(t_specs.name + " bias threshold is nonsensical");
1676 : }
1677 3 : parse(t_specs.name_stem + "ALPHA", t_specs.alpha);
1678 3 : if (t_specs.alpha <= 0.0 || t_specs.alpha > 1.0) {
1679 0 : error(t_specs.name + " decay shape parameter alpha is nonsensical");
1680 : }
1681 : }
1682 151 : }
1683 :
1684 3 : void MetaD::logTemperingSpecs(const TemperingSpecs &t_specs) {
1685 3 : log.printf(" %s bias factor %f\n", t_specs.name.c_str(), t_specs.biasf);
1686 3 : log.printf(" KbT %f\n", kbt_);
1687 3 : if (t_specs.threshold != 0.0) {
1688 2 : log.printf(" %s bias threshold %f\n", t_specs.name.c_str(), t_specs.threshold);
1689 : }
1690 3 : if (t_specs.alpha != 1.0) {
1691 1 : log.printf(" %s decay shape parameter alpha %f\n", t_specs.name.c_str(), t_specs.alpha);
1692 : }
1693 3 : }
1694 :
1695 6036 : void MetaD::readGaussians(IFile *ifile) {
1696 6036 : unsigned ncv=getNumberOfArguments();
1697 6036 : std::vector<double> center(ncv);
1698 6036 : std::vector<double> sigma(ncv);
1699 : double height;
1700 : int nhills=0;
1701 6036 : bool multivariate=false;
1702 :
1703 : std::vector<Value> tmpvalues;
1704 18121 : for(unsigned j=0; j<getNumberOfArguments(); ++j) {
1705 24170 : tmpvalues.push_back( Value( this, getPntrToArgument(j)->getName(), false ) );
1706 : }
1707 :
1708 8569 : while(scanOneHill(ifile,tmpvalues,center,sigma,height,multivariate)) {
1709 2533 : nhills++;
1710 : // note that for gamma=1 we store directly -F
1711 2533 : if(welltemp_ && biasf_>1.0) {
1712 41 : height*=(biasf_-1.0)/biasf_;
1713 : }
1714 2533 : addGaussian(Gaussian(multivariate,height,center,sigma));
1715 : }
1716 6036 : log.printf(" %d Gaussians read\n",nhills);
1717 12072 : }
1718 :
1719 2922 : void MetaD::writeGaussian(const Gaussian& hill, OFile&file) {
1720 2922 : unsigned ncv=getNumberOfArguments();
1721 2922 : file.printField("time",getTimeStep()*getStep());
1722 8194 : for(unsigned i=0; i<ncv; ++i) {
1723 5272 : file.printField(getPntrToArgument(i),hill.center[i]);
1724 : }
1725 5844 : hillsOfile_.printField("kerneltype","stretched-gaussian");
1726 2922 : if(hill.multivariate) {
1727 892 : hillsOfile_.printField("multivariate","true");
1728 : Matrix<double> mymatrix(ncv,ncv);
1729 : unsigned k=0;
1730 1047 : for(unsigned i=0; i<ncv; i++) {
1731 1357 : for(unsigned j=i; j<ncv; j++) {
1732 : // recompose the full inverse matrix
1733 756 : mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k];
1734 756 : k++;
1735 : }
1736 : }
1737 : // invert it
1738 : Matrix<double> invmatrix(ncv,ncv);
1739 446 : Invert(mymatrix,invmatrix);
1740 : // enforce symmetry
1741 1047 : for(unsigned i=0; i<ncv; i++) {
1742 1357 : for(unsigned j=i; j<ncv; j++) {
1743 756 : invmatrix(i,j)=invmatrix(j,i);
1744 : }
1745 : }
1746 :
1747 : // do cholesky so to have a "sigma like" number
1748 : Matrix<double> lower(ncv,ncv);
1749 446 : cholesky(invmatrix,lower);
1750 : // loop in band form
1751 1047 : for(unsigned i=0; i<ncv; i++) {
1752 1357 : for(unsigned j=0; j<ncv-i; j++) {
1753 1512 : file.printField("sigma_"+getPntrToArgument(j+i)->getName()+"_"+getPntrToArgument(j)->getName(),lower(j+i,j));
1754 : }
1755 : }
1756 : } else {
1757 4952 : hillsOfile_.printField("multivariate","false");
1758 7147 : for(unsigned i=0; i<ncv; ++i) {
1759 9342 : file.printField("sigma_"+getPntrToArgument(i)->getName(),hill.sigma[i]);
1760 : }
1761 : }
1762 2922 : double height=hill.height;
1763 : // note that for gamma=1 we store directly -F
1764 2922 : if(welltemp_ && biasf_>1.0) {
1765 339 : height*=biasf_/(biasf_-1.0);
1766 : }
1767 5844 : file.printField("height",height).printField("biasf",biasf_);
1768 2922 : if(mw_n_>1) {
1769 3018 : file.printField("clock",int(std::time(0)));
1770 : }
1771 2922 : file.printField();
1772 2922 : }
1773 :
1774 5617 : void MetaD::addGaussian(const Gaussian& hill) {
1775 5617 : if(grid_) {
1776 640 : size_t ncv=getNumberOfArguments();
1777 640 : std::vector<unsigned> nneighb=getGaussianSupport(hill);
1778 640 : std::vector<Grid::index_t> neighbors=BiasGrid_->getNeighbors(hill.center,nneighb);
1779 640 : std::vector<double> der(ncv);
1780 640 : std::vector<double> xx(ncv);
1781 640 : if(comm.Get_size()==1) {
1782 : // for performance reasons and thread safety
1783 544 : std::vector<double> dp(ncv);
1784 55324 : for(size_t i=0; i<neighbors.size(); ++i) {
1785 54780 : Grid::index_t ineigh=neighbors[i];
1786 158922 : for(unsigned j=0; j<ncv; ++j) {
1787 104142 : der[j]=0.0;
1788 : }
1789 54780 : BiasGrid_->getPoint(ineigh,xx);
1790 54780 : double bias=evaluateGaussianAndDerivatives(xx,hill,der,dp);
1791 54780 : BiasGrid_->addValueAndDerivatives(ineigh,bias,der);
1792 : }
1793 : } else {
1794 96 : unsigned stride=comm.Get_size();
1795 96 : unsigned rank=comm.Get_rank();
1796 96 : std::vector<double> allder(ncv*neighbors.size(),0.0);
1797 96 : std::vector<double> n_der(ncv,0.0);
1798 96 : std::vector<double> allbias(neighbors.size(),0.0);
1799 : // for performance reasons and thread safety
1800 96 : std::vector<double> dp(ncv);
1801 27148 : for(unsigned i=rank; i<neighbors.size(); i+=stride) {
1802 27052 : Grid::index_t ineigh=neighbors[i];
1803 81156 : for(unsigned j=0; j<ncv; ++j) {
1804 54104 : n_der[j]=0.0;
1805 : }
1806 27052 : BiasGrid_->getPoint(ineigh,xx);
1807 27052 : allbias[i]=evaluateGaussianAndDerivatives(xx,hill,n_der,dp);
1808 81156 : for(unsigned j=0; j<ncv; j++) {
1809 54104 : allder[ncv*i+j]=n_der[j];
1810 : }
1811 : }
1812 96 : comm.Sum(allbias);
1813 96 : comm.Sum(allder);
1814 103200 : for(unsigned i=0; i<neighbors.size(); ++i) {
1815 103104 : Grid::index_t ineigh=neighbors[i];
1816 309312 : for(unsigned j=0; j<ncv; ++j) {
1817 206208 : der[j]=allder[ncv*i+j];
1818 : }
1819 103104 : BiasGrid_->addValueAndDerivatives(ineigh,allbias[i],der);
1820 : }
1821 : }
1822 : } else {
1823 4977 : hills_.push_back(hill);
1824 : }
1825 5617 : }
1826 :
1827 640 : std::vector<unsigned> MetaD::getGaussianSupport(const Gaussian& hill) {
1828 : std::vector<unsigned> nneigh;
1829 : std::vector<double> cutoff;
1830 640 : unsigned ncv=getNumberOfArguments();
1831 :
1832 : // traditional or flexible hill?
1833 640 : if(hill.multivariate) {
1834 : unsigned k=0;
1835 : Matrix<double> mymatrix(ncv,ncv);
1836 0 : for(unsigned i=0; i<ncv; i++) {
1837 0 : for(unsigned j=i; j<ncv; j++) {
1838 : // recompose the full inverse matrix
1839 0 : mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k];
1840 0 : k++;
1841 : }
1842 : }
1843 : // Reinvert so to have the ellipses
1844 : Matrix<double> myinv(ncv,ncv);
1845 0 : Invert(mymatrix,myinv);
1846 : Matrix<double> myautovec(ncv,ncv);
1847 0 : std::vector<double> myautoval(ncv); //should I take this or their square root?
1848 0 : diagMat(myinv,myautoval,myautovec);
1849 : double maxautoval=0.;
1850 : unsigned ind_maxautoval;
1851 : ind_maxautoval=ncv;
1852 0 : for(unsigned i=0; i<ncv; i++) {
1853 0 : if(myautoval[i]>maxautoval) {
1854 : maxautoval=myautoval[i];
1855 : ind_maxautoval=i;
1856 : }
1857 : }
1858 0 : for(unsigned i=0; i<ncv; i++) {
1859 0 : cutoff.push_back(std::sqrt(2.0*dp2cutoff)*std::abs(std::sqrt(maxautoval)*myautovec(i,ind_maxautoval)));
1860 : }
1861 : } else {
1862 1618 : for(unsigned i=0; i<ncv; ++i) {
1863 978 : cutoff.push_back(std::sqrt(2.0*dp2cutoff)*hill.sigma[i]);
1864 : }
1865 : }
1866 :
1867 640 : if(doInt_) {
1868 2 : if(hill.center[0]+cutoff[0] > uppI_ || hill.center[0]-cutoff[0] < lowI_) {
1869 : // in this case, we updated the entire grid to avoid problems
1870 2 : return BiasGrid_->getNbin();
1871 : } else {
1872 0 : nneigh.push_back( static_cast<unsigned>(ceil(cutoff[0]/BiasGrid_->getDx()[0])) );
1873 : return nneigh;
1874 : }
1875 : } else {
1876 1614 : for(unsigned i=0; i<ncv; i++) {
1877 976 : nneigh.push_back( static_cast<unsigned>(ceil(cutoff[i]/BiasGrid_->getDx()[i])) );
1878 : }
1879 : }
1880 :
1881 : return nneigh;
1882 : }
1883 :
1884 285 : double MetaD::getBias(const std::vector<double>& cv) {
1885 285 : double bias=0.0;
1886 285 : if(grid_) {
1887 203 : bias = BiasGrid_->getValue(cv);
1888 : } else {
1889 82 : unsigned nt=OpenMP::getNumThreads();
1890 82 : unsigned stride=comm.Get_size();
1891 82 : unsigned rank=comm.Get_rank();
1892 :
1893 82 : if(!nlist_) {
1894 82 : #pragma omp parallel num_threads(nt)
1895 : {
1896 : #pragma omp for reduction(+:bias) nowait
1897 : for(unsigned i=rank; i<hills_.size(); i+=stride) {
1898 : bias+=evaluateGaussian(cv,hills_[i]);
1899 : }
1900 : }
1901 : } else {
1902 0 : #pragma omp parallel num_threads(nt)
1903 : {
1904 : #pragma omp for reduction(+:bias) nowait
1905 : for(unsigned i=rank; i<nlist_hills_.size(); i+=stride) {
1906 : bias+=evaluateGaussian(cv,nlist_hills_[i]);
1907 : }
1908 : }
1909 : }
1910 82 : comm.Sum(bias);
1911 : }
1912 :
1913 285 : return bias;
1914 : }
1915 :
1916 8395 : double MetaD::getBiasAndDerivatives(const std::vector<double>& cv, std::vector<double>& der) {
1917 8395 : unsigned ncv=getNumberOfArguments();
1918 8395 : double bias=0.0;
1919 8395 : if(grid_) {
1920 1506 : std::vector<double> vder(ncv);
1921 1506 : bias=BiasGrid_->getValueAndDerivatives(cv,vder);
1922 3498 : for(unsigned i=0; i<ncv; i++) {
1923 1992 : der[i]=vder[i];
1924 : }
1925 : } else {
1926 6889 : unsigned nt=OpenMP::getNumThreads();
1927 6889 : unsigned stride=comm.Get_size();
1928 6889 : unsigned rank=comm.Get_rank();
1929 :
1930 6889 : if(!nlist_) {
1931 6884 : if(hills_.size()<2*nt*stride||nt==1) {
1932 : // for performance reasons and thread safety
1933 2587 : std::vector<double> dp(ncv);
1934 6703 : for(unsigned i=rank; i<hills_.size(); i+=stride) {
1935 4116 : bias+=evaluateGaussianAndDerivatives(cv,hills_[i],der,dp);
1936 : }
1937 : } else {
1938 4297 : #pragma omp parallel num_threads(nt)
1939 : {
1940 : std::vector<double> omp_deriv(ncv,0.);
1941 : // for performance reasons and thread safety
1942 : std::vector<double> dp(ncv);
1943 : #pragma omp for reduction(+:bias) nowait
1944 : for(unsigned i=rank; i<hills_.size(); i+=stride) {
1945 : bias+=evaluateGaussianAndDerivatives(cv,hills_[i],omp_deriv,dp);
1946 : }
1947 : #pragma omp critical
1948 : for(unsigned i=0; i<ncv; i++) {
1949 : der[i]+=omp_deriv[i];
1950 : }
1951 : }
1952 : }
1953 : } else {
1954 5 : if(hills_.size()<2*nt*stride||nt==1) {
1955 : // for performance reasons and thread safety
1956 0 : std::vector<double> dp(ncv);
1957 0 : for(unsigned i=rank; i<nlist_hills_.size(); i+=stride) {
1958 0 : bias+=evaluateGaussianAndDerivatives(cv,nlist_hills_[i],der,dp);
1959 : }
1960 : } else {
1961 5 : #pragma omp parallel num_threads(nt)
1962 : {
1963 : std::vector<double> omp_deriv(ncv,0.);
1964 : // for performance reasons and thread safety
1965 : std::vector<double> dp(ncv);
1966 : #pragma omp for reduction(+:bias) nowait
1967 : for(unsigned i=rank; i<nlist_hills_.size(); i+=stride) {
1968 : bias+=evaluateGaussianAndDerivatives(cv,nlist_hills_[i],omp_deriv,dp);
1969 : }
1970 : #pragma omp critical
1971 : for(unsigned i=0; i<ncv; i++) {
1972 : der[i]+=omp_deriv[i];
1973 : }
1974 : }
1975 : }
1976 : }
1977 6889 : comm.Sum(bias);
1978 6889 : comm.Sum(der);
1979 : }
1980 :
1981 8395 : return bias;
1982 : }
1983 :
1984 0 : double MetaD::getGaussianNormalization(const Gaussian& hill) {
1985 : double norm=1;
1986 0 : unsigned ncv=hill.center.size();
1987 :
1988 0 : if(hill.multivariate) {
1989 : // recompose the full sigma from the upper diag cholesky
1990 : unsigned k=0;
1991 : Matrix<double> mymatrix(ncv,ncv);
1992 0 : for(unsigned i=0; i<ncv; i++) {
1993 0 : for(unsigned j=i; j<ncv; j++) {
1994 0 : mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k]; // recompose the full inverse matrix
1995 0 : k++;
1996 : }
1997 : double ldet;
1998 0 : logdet( mymatrix, ldet );
1999 0 : norm = std::exp( ldet ); // Not sure here if mymatrix is sigma or inverse
2000 : }
2001 : } else {
2002 0 : for(unsigned i=0; i<hill.sigma.size(); i++) {
2003 0 : norm*=hill.sigma[i];
2004 : }
2005 : }
2006 :
2007 0 : return norm*std::pow(2*pi,static_cast<double>(ncv)/2.0);
2008 : }
2009 :
2010 192 : double MetaD::evaluateGaussian(const std::vector<double>& cv, const Gaussian& hill) {
2011 192 : unsigned ncv=cv.size();
2012 :
2013 : // I use a pointer here because cv is const (and should be const)
2014 : // but when using doInt it is easier to locally replace cv[0] with
2015 : // the upper/lower limit in case it is out of range
2016 : double tmpcv[1];
2017 : const double *pcv=NULL; // pointer to cv
2018 192 : if(ncv>0) {
2019 : pcv=&cv[0];
2020 : }
2021 192 : if(doInt_) {
2022 0 : plumed_assert(ncv==1);
2023 0 : tmpcv[0]=cv[0];
2024 0 : if(cv[0]<lowI_) {
2025 0 : tmpcv[0]=lowI_;
2026 : }
2027 0 : if(cv[0]>uppI_) {
2028 0 : tmpcv[0]=uppI_;
2029 : }
2030 : pcv=&(tmpcv[0]);
2031 : }
2032 :
2033 : double dp2=0.0;
2034 192 : if(hill.multivariate) {
2035 : unsigned k=0;
2036 : // recompose the full sigma from the upper diag cholesky
2037 : Matrix<double> mymatrix(ncv,ncv);
2038 0 : for(unsigned i=0; i<ncv; i++) {
2039 0 : for(unsigned j=i; j<ncv; j++) {
2040 0 : mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k]; // recompose the full inverse matrix
2041 0 : k++;
2042 : }
2043 : }
2044 0 : for(unsigned i=0; i<ncv; i++) {
2045 0 : double dp_i=difference(i,hill.center[i],pcv[i]);
2046 0 : for(unsigned j=i; j<ncv; j++) {
2047 0 : if(i==j) {
2048 0 : dp2+=dp_i*dp_i*mymatrix(i,j)*0.5;
2049 : } else {
2050 0 : double dp_j=difference(j,hill.center[j],pcv[j]);
2051 0 : dp2+=dp_i*dp_j*mymatrix(i,j);
2052 : }
2053 : }
2054 : }
2055 : } else {
2056 576 : for(unsigned i=0; i<ncv; i++) {
2057 384 : double dp=difference(i,hill.center[i],pcv[i])*hill.invsigma[i];
2058 384 : dp2+=dp*dp;
2059 : }
2060 192 : dp2*=0.5;
2061 : }
2062 :
2063 : double bias=0.0;
2064 192 : if(dp2<dp2cutoff) {
2065 154 : bias=hill.height*(stretchA*std::exp(-dp2)+stretchB);
2066 : }
2067 :
2068 192 : return bias;
2069 : }
2070 :
2071 2408462 : double MetaD::evaluateGaussianAndDerivatives(const std::vector<double>& cv, const Gaussian& hill, std::vector<double>& der, std::vector<double>& dp_) {
2072 2408462 : unsigned ncv=cv.size();
2073 :
2074 : // I use a pointer here because cv is const (and should be const)
2075 : // but when using doInt it is easier to locally replace cv[0] with
2076 : // the upper/lower limit in case it is out of range
2077 : const double *pcv=NULL; // pointer to cv
2078 : double tmpcv[1]; // tmp array with cv (to be used with doInt_)
2079 2408462 : if(ncv>0) {
2080 : pcv=&cv[0];
2081 : }
2082 2408462 : if(doInt_) {
2083 602 : plumed_assert(ncv==1);
2084 602 : tmpcv[0]=cv[0];
2085 602 : if(cv[0]<lowI_) {
2086 118 : tmpcv[0]=lowI_;
2087 : }
2088 602 : if(cv[0]>uppI_) {
2089 360 : tmpcv[0]=uppI_;
2090 : }
2091 : pcv=&(tmpcv[0]);
2092 : }
2093 :
2094 : bool int_der=false;
2095 2408462 : if(doInt_) {
2096 602 : if(cv[0]<lowI_ || cv[0]>uppI_) {
2097 : int_der=true;
2098 : }
2099 : }
2100 :
2101 : double dp2=0.0;
2102 : double bias=0.0;
2103 2408462 : if(hill.multivariate) {
2104 : unsigned k=0;
2105 : // recompose the full sigma from the upper diag cholesky
2106 : Matrix<double> mymatrix(ncv,ncv);
2107 161635 : for(unsigned i=0; i<ncv; i++) {
2108 162513 : for(unsigned j=i; j<ncv; j++) {
2109 81476 : mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k]; // recompose the full inverse matrix
2110 81476 : k++;
2111 : }
2112 : }
2113 161635 : for(unsigned i=0; i<ncv; i++) {
2114 81037 : dp_[i]=difference(i,hill.center[i],pcv[i]);
2115 162513 : for(unsigned j=i; j<ncv; j++) {
2116 81476 : if(i==j) {
2117 81037 : dp2+=dp_[i]*dp_[i]*mymatrix(i,j)*0.5;
2118 : } else {
2119 439 : double dp_j=difference(j,hill.center[j],pcv[j]);
2120 439 : dp2+=dp_[i]*dp_j*mymatrix(i,j);
2121 : }
2122 : }
2123 : }
2124 80598 : if(dp2<dp2cutoff) {
2125 77683 : bias=hill.height*std::exp(-dp2);
2126 77683 : if(!int_der) {
2127 155673 : for(unsigned i=0; i<ncv; i++) {
2128 : double tmp=0.0;
2129 156594 : for(unsigned j=0; j<ncv; j++) {
2130 78604 : tmp += dp_[j]*mymatrix(i,j)*bias;
2131 : }
2132 77990 : der[i]-=tmp*stretchA;
2133 : }
2134 : } else {
2135 0 : for(unsigned i=0; i<ncv; i++) {
2136 0 : der[i]=0.;
2137 : }
2138 : }
2139 77683 : bias=stretchA*bias+hill.height*stretchB;
2140 : }
2141 : } else {
2142 6973950 : for(unsigned i=0; i<ncv; i++) {
2143 4646086 : dp_[i]=difference(i,hill.center[i],pcv[i])*hill.invsigma[i];
2144 4646086 : dp2+=dp_[i]*dp_[i];
2145 : }
2146 2327864 : dp2*=0.5;
2147 2327864 : if(dp2<dp2cutoff) {
2148 1355742 : bias=hill.height*std::exp(-dp2);
2149 1355742 : if(!int_der) {
2150 4058447 : for(unsigned i=0; i<ncv; i++) {
2151 2702944 : der[i]-=bias*dp_[i]*hill.invsigma[i]*stretchA;
2152 : }
2153 : } else {
2154 478 : for(unsigned i=0; i<ncv; i++) {
2155 239 : der[i]=0.;
2156 : }
2157 : }
2158 1355742 : bias=stretchA*bias+hill.height*stretchB;
2159 : }
2160 : }
2161 :
2162 2408462 : return bias;
2163 : }
2164 :
2165 2736 : double MetaD::getHeight(const std::vector<double>& cv) {
2166 2736 : double height=height0_;
2167 2736 : if(welltemp_) {
2168 275 : double vbias = getBias(cv);
2169 275 : if(biasf_>1.0) {
2170 259 : height = height0_*std::exp(-vbias/(kbt_*(biasf_-1.0)));
2171 : } else {
2172 : // notice that if gamma=1 we store directly -F
2173 16 : height = height0_*std::exp(-vbias/kbt_);
2174 : }
2175 : }
2176 2736 : if(dampfactor_>0.0) {
2177 18 : plumed_assert(BiasGrid_);
2178 18 : double m=BiasGrid_->getMaxValue();
2179 18 : height*=std::exp(-m/(kbt_*(dampfactor_)));
2180 : }
2181 2736 : if (tt_specs_.is_active) {
2182 60 : double vbarrier = transition_bias_;
2183 60 : temperHeight(height, tt_specs_, vbarrier);
2184 : }
2185 2736 : if(TargetGrid_) {
2186 18 : double f=TargetGrid_->getValue(cv)-TargetGrid_->getMaxValue();
2187 18 : height*=std::exp(f/kbt_);
2188 : }
2189 2736 : return height;
2190 : }
2191 :
2192 60 : void MetaD::temperHeight(double& height, const TemperingSpecs& t_specs, const double tempering_bias) {
2193 60 : if (t_specs.alpha == 1.0) {
2194 80 : height *= std::exp(-std::max(0.0, tempering_bias - t_specs.threshold) / (kbt_ * (t_specs.biasf - 1.0)));
2195 : } else {
2196 40 : height *= std::pow(1 + (1 - t_specs.alpha) / t_specs.alpha * std::max(0.0, tempering_bias - t_specs.threshold) / (kbt_ * (t_specs.biasf - 1.0)), - t_specs.alpha / (1 - t_specs.alpha));
2197 : }
2198 60 : }
2199 :
2200 8435 : void MetaD::calculate() {
2201 : // this is because presently there is no way to properly pass information
2202 : // on adaptive hills (diff) after exchanges:
2203 8435 : if(adaptive_==FlexibleBin::diffusion && getExchangeStep()) {
2204 0 : error("ADAPTIVE=DIFF is not compatible with replica exchange");
2205 : }
2206 :
2207 8435 : const unsigned ncv=getNumberOfArguments();
2208 8435 : std::vector<double> cv(ncv);
2209 21082 : for(unsigned i=0; i<ncv; ++i) {
2210 12647 : cv[i]=getArgument(i);
2211 : }
2212 :
2213 8435 : if(nlist_) {
2214 5 : nlist_steps_++;
2215 5 : if(getExchangeStep()) {
2216 0 : nlist_update_=true;
2217 : } else {
2218 11 : for(unsigned i=0; i<ncv; ++i) {
2219 8 : double d = difference(i, cv[i], nlist_center_[i]);
2220 8 : double nk_dist2 = d*d/nlist_dev2_[i];
2221 8 : if(nk_dist2>nlist_param_[1]) {
2222 2 : nlist_update_=true;
2223 2 : break;
2224 : }
2225 : }
2226 : }
2227 5 : if(nlist_update_) {
2228 4 : updateNlist();
2229 : }
2230 : }
2231 :
2232 : double ene = 0.;
2233 8435 : std::vector<double> der(ncv,0.);
2234 8435 : if(biasf_!=1.0) {
2235 8395 : ene = getBiasAndDerivatives(cv,der);
2236 : }
2237 : setBias(ene);
2238 21082 : for(unsigned i=0; i<ncv; i++) {
2239 12647 : setOutputForce(i,-der[i]);
2240 : }
2241 :
2242 8435 : if(calc_work_) {
2243 10 : getPntrToComponent("work")->set(work_);
2244 : }
2245 8435 : if(calc_rct_) {
2246 220 : getPntrToComponent("rbias")->set(ene - reweight_factor_);
2247 : }
2248 : // calculate the acceleration factor
2249 8435 : if(acceleration_&&!isFirstStep_) {
2250 329 : acc_ += static_cast<double>(getStride()) * std::exp(ene/(kbt_));
2251 329 : const double mean_acc = acc_/((double) getStep());
2252 329 : getPntrToComponent("acc")->set(mean_acc);
2253 8435 : } else if (acceleration_ && isFirstStep_ && acc_restart_mean_ > 0.0) {
2254 2 : acc_ = acc_restart_mean_ * static_cast<double>(getStep());
2255 2 : if(freq_adaptive_) {
2256 : // has to be done here if restarting, as the acc is not defined before
2257 1 : updateFrequencyAdaptiveStride();
2258 : }
2259 : }
2260 8435 : }
2261 :
2262 6239 : void MetaD::update() {
2263 : // adding hills criteria (could be more complex though)
2264 : bool nowAddAHill;
2265 6239 : if(getStep()%current_stride_==0 && !isFirstStep_) {
2266 : nowAddAHill=true;
2267 : } else {
2268 : nowAddAHill=false;
2269 3503 : isFirstStep_=false;
2270 : }
2271 :
2272 6239 : unsigned ncv=getNumberOfArguments();
2273 6239 : std::vector<double> cv(ncv);
2274 16690 : for(unsigned i=0; i<ncv; ++i) {
2275 10451 : cv[i] = getArgument(i);
2276 : }
2277 :
2278 : double vbias=0.;
2279 6239 : if(calc_work_) {
2280 5 : vbias=getBias(cv);
2281 : }
2282 :
2283 : // if you use adaptive, call the FlexibleBin
2284 : bool multivariate=false;
2285 6239 : if(adaptive_!=FlexibleBin::none) {
2286 778 : flexbin_->update(nowAddAHill);
2287 : multivariate=true;
2288 : }
2289 :
2290 : std::vector<double> thissigma;
2291 6239 : if(nowAddAHill) {
2292 : // add a Gaussian
2293 2736 : double height=getHeight(cv);
2294 : // returns upper diagonal inverse
2295 2736 : if(adaptive_!=FlexibleBin::none) {
2296 748 : thissigma=flexbin_->getInverseMatrix();
2297 : }
2298 : // returns normal sigma
2299 : else {
2300 2362 : thissigma=sigma0_;
2301 : }
2302 :
2303 : // In case we use walkers_mpi, it is now necessary to communicate with other replicas.
2304 2736 : if(walkers_mpi_) {
2305 : // Allocate arrays to store all walkers hills
2306 174 : std::vector<double> all_cv(mpi_nw_*ncv,0.0);
2307 174 : std::vector<double> all_sigma(mpi_nw_*thissigma.size(),0.0);
2308 174 : std::vector<double> all_height(mpi_nw_,0.0);
2309 174 : std::vector<int> all_multivariate(mpi_nw_,0);
2310 174 : if(comm.Get_rank()==0) {
2311 : // Communicate (only root)
2312 99 : multi_sim_comm.Allgather(cv,all_cv);
2313 99 : multi_sim_comm.Allgather(thissigma,all_sigma);
2314 : // notice that if gamma=1 we store directly -F so this scaling is not necessary:
2315 99 : multi_sim_comm.Allgather(height*(biasf_>1.0?biasf_/(biasf_-1.0):1.0),all_height);
2316 99 : multi_sim_comm.Allgather(int(multivariate),all_multivariate);
2317 : }
2318 : // Share info with group members
2319 174 : comm.Bcast(all_cv,0);
2320 174 : comm.Bcast(all_sigma,0);
2321 174 : comm.Bcast(all_height,0);
2322 174 : comm.Bcast(all_multivariate,0);
2323 :
2324 : // Flying Gaussian
2325 174 : if (flying_) {
2326 54 : hills_.clear();
2327 54 : comm.Barrier();
2328 : }
2329 :
2330 696 : for(unsigned i=0; i<mpi_nw_; i++) {
2331 : // actually add hills one by one
2332 522 : std::vector<double> cv_now(ncv);
2333 522 : std::vector<double> sigma_now(thissigma.size());
2334 1566 : for(unsigned j=0; j<ncv; j++) {
2335 1044 : cv_now[j]=all_cv[i*ncv+j];
2336 : }
2337 1674 : for(unsigned j=0; j<thissigma.size(); j++) {
2338 1152 : sigma_now[j]=all_sigma[i*thissigma.size()+j];
2339 : }
2340 : // notice that if gamma=1 we store directly -F so this scaling is not necessary:
2341 522 : double fact=(biasf_>1.0?(biasf_-1.0)/biasf_:1.0);
2342 522 : Gaussian newhill=Gaussian(all_multivariate[i],all_height[i]*fact,cv_now,sigma_now);
2343 522 : addGaussian(newhill);
2344 522 : if(!flying_) {
2345 360 : writeGaussian(newhill,hillsOfile_);
2346 : }
2347 522 : }
2348 : } else {
2349 2562 : Gaussian newhill=Gaussian(multivariate,height,cv,thissigma);
2350 2562 : addGaussian(newhill);
2351 2562 : writeGaussian(newhill,hillsOfile_);
2352 2562 : }
2353 :
2354 : // this is to update the hills neighbor list
2355 2736 : if(nlist_) {
2356 4 : nlist_update_=true;
2357 : }
2358 : }
2359 :
2360 : // this should be outside of the if block in case
2361 : // mw_rstride_ is not a multiple of stride_
2362 6239 : if(mw_n_>1 && getStep()%mw_rstride_==0) {
2363 3012 : hillsOfile_.flush();
2364 : }
2365 :
2366 6239 : if(calc_work_) {
2367 5 : if(nlist_) {
2368 0 : updateNlist();
2369 : }
2370 5 : double vbias1=getBias(cv);
2371 5 : work_+=vbias1-vbias;
2372 : }
2373 :
2374 : // dump grid on file
2375 6239 : if(wgridstride_>0&&(getStep()%wgridstride_==0||getCPT())) {
2376 : // in case old grids are stored, a sequence of grids should appear
2377 : // this call results in a repetition of the header:
2378 91 : if(storeOldGrids_) {
2379 40 : gridfile_.clearFields();
2380 : }
2381 : // in case only latest grid is stored, file should be rewound
2382 : // this will overwrite previously written grids
2383 : else {
2384 51 : int r = 0;
2385 51 : if(walkers_mpi_) {
2386 0 : if(comm.Get_rank()==0) {
2387 0 : r=multi_sim_comm.Get_rank();
2388 : }
2389 0 : comm.Bcast(r,0);
2390 : }
2391 51 : if(r==0) {
2392 51 : gridfile_.rewind();
2393 : }
2394 : }
2395 91 : BiasGrid_->writeToFile(gridfile_);
2396 : // if a single grid is stored, it is necessary to flush it, otherwise
2397 : // the file might stay empty forever (when a single grid is not large enough to
2398 : // trigger flushing from the operating system).
2399 : // on the other hand, if grids are stored one after the other this is
2400 : // no necessary, and we leave the flushing control to the user as usual
2401 : // (with FLUSH keyword)
2402 91 : if(!storeOldGrids_) {
2403 51 : gridfile_.flush();
2404 : }
2405 : }
2406 :
2407 : // if multiple walkers and time to read Gaussians
2408 6239 : if(mw_n_>1 && getStep()%mw_rstride_==0) {
2409 12048 : for(int i=0; i<mw_n_; ++i) {
2410 : // don't read your own Gaussians
2411 9036 : if(i==mw_id_) {
2412 3012 : continue;
2413 : }
2414 : // if the file is not open yet
2415 6024 : if(!(ifiles_[i]->isOpen())) {
2416 : // check if it exists now and open it!
2417 7 : if(ifiles_[i]->FileExist(ifilesnames_[i])) {
2418 7 : ifiles_[i]->open(ifilesnames_[i]);
2419 7 : ifiles_[i]->reset(false);
2420 : }
2421 : // otherwise read the new Gaussians
2422 : } else {
2423 6017 : log.printf(" Reading hills from %s:",ifilesnames_[i].c_str());
2424 6017 : readGaussians(ifiles_[i].get());
2425 6017 : ifiles_[i]->reset(false);
2426 : }
2427 : }
2428 : // this is to update the hills neighbor list
2429 3012 : if(nlist_) {
2430 0 : nlist_update_=true;
2431 : }
2432 : }
2433 :
2434 : // Recalculate special bias quantities whenever the bias has been changed by the update.
2435 6239 : bool bias_has_changed = (nowAddAHill || (mw_n_ > 1 && getStep() % mw_rstride_ == 0));
2436 6239 : if (calc_rct_ && bias_has_changed && getStep()%(stride_*rct_ustride_)==0) {
2437 102 : computeReweightingFactor();
2438 : }
2439 6239 : if (calc_max_bias_ && bias_has_changed) {
2440 0 : max_bias_ = BiasGrid_->getMaxValue();
2441 0 : getPntrToComponent("maxbias")->set(max_bias_);
2442 : }
2443 6239 : if (calc_transition_bias_ && bias_has_changed) {
2444 260 : transition_bias_ = getTransitionBarrierBias();
2445 520 : getPntrToComponent("transbias")->set(transition_bias_);
2446 : }
2447 :
2448 : // Frequency adaptive metadynamics - update hill addition frequency
2449 6239 : if(freq_adaptive_ && getStep()%fa_update_frequency_==0) {
2450 151 : updateFrequencyAdaptiveStride();
2451 : }
2452 6239 : }
2453 :
2454 : /// takes a pointer to the file and a template std::string with values v and gives back the next center, sigma and height
2455 8569 : bool MetaD::scanOneHill(IFile* ifile, std::vector<Value>& tmpvalues, std::vector<double>& center, std::vector<double>& sigma, double& height, bool& multivariate) {
2456 : double dummy;
2457 8569 : multivariate=false;
2458 17138 : if(ifile->scanField("time",dummy)) {
2459 2533 : unsigned ncv=tmpvalues.size();
2460 7553 : for(unsigned i=0; i<ncv; ++i) {
2461 5020 : ifile->scanField( &tmpvalues[i] );
2462 5020 : if( tmpvalues[i].isPeriodic() && ! getPntrToArgument(i)->isPeriodic() ) {
2463 0 : error("in hills file periodicity for variable " + tmpvalues[i].getName() + " does not match periodicity in input");
2464 5020 : } else if( tmpvalues[i].isPeriodic() ) {
2465 : std::string imin, imax;
2466 0 : tmpvalues[i].getDomain( imin, imax );
2467 : std::string rmin, rmax;
2468 0 : getPntrToArgument(i)->getDomain( rmin, rmax );
2469 0 : if( imin!=rmin || imax!=rmax ) {
2470 0 : error("in hills file periodicity for variable " + tmpvalues[i].getName() + " does not match periodicity in input");
2471 : }
2472 : }
2473 5020 : center[i]=tmpvalues[i].get();
2474 : }
2475 : // scan for kerneltype
2476 2533 : std::string ktype="stretched-gaussian";
2477 5066 : if( ifile->FieldExist("kerneltype") ) {
2478 5048 : ifile->scanField("kerneltype",ktype);
2479 : }
2480 2533 : if( ktype=="gaussian" ) {
2481 12 : noStretchWarning();
2482 2521 : } else if( ktype!="stretched-gaussian") {
2483 0 : error("non Gaussian kernels are not supported in MetaD");
2484 : }
2485 : // scan for multivariate label: record the actual file position so to eventually rewind
2486 : std::string sss;
2487 5066 : ifile->scanField("multivariate",sss);
2488 2533 : if(sss=="true") {
2489 0 : multivariate=true;
2490 2533 : } else if(sss=="false") {
2491 2533 : multivariate=false;
2492 : } else {
2493 0 : plumed_merror("cannot parse multivariate = "+ sss);
2494 : }
2495 2533 : if(multivariate) {
2496 0 : sigma.resize(ncv*(ncv+1)/2);
2497 : Matrix<double> upper(ncv,ncv);
2498 : Matrix<double> lower(ncv,ncv);
2499 0 : for(unsigned i=0; i<ncv; i++) {
2500 0 : for(unsigned j=0; j<ncv-i; j++) {
2501 0 : ifile->scanField("sigma_"+getPntrToArgument(j+i)->getName()+"_"+getPntrToArgument(j)->getName(),lower(j+i,j));
2502 0 : upper(j,j+i)=lower(j+i,j);
2503 : }
2504 : }
2505 : Matrix<double> mymult(ncv,ncv);
2506 : Matrix<double> invmatrix(ncv,ncv);
2507 0 : mult(lower,upper,mymult);
2508 : // now invert and get the sigmas
2509 0 : Invert(mymult,invmatrix);
2510 : // put the sigmas in the usual order: upper diagonal (this time in normal form and not in band form)
2511 : unsigned k=0;
2512 0 : for(unsigned i=0; i<ncv; i++) {
2513 0 : for(unsigned j=i; j<ncv; j++) {
2514 0 : sigma[k]=invmatrix(i,j);
2515 0 : k++;
2516 : }
2517 : }
2518 : } else {
2519 7553 : for(unsigned i=0; i<ncv; ++i) {
2520 10040 : ifile->scanField("sigma_"+getPntrToArgument(i)->getName(),sigma[i]);
2521 : }
2522 : }
2523 :
2524 2533 : ifile->scanField("height",height);
2525 2533 : ifile->scanField("biasf",dummy);
2526 5066 : if(ifile->FieldExist("clock")) {
2527 4700 : ifile->scanField("clock",dummy);
2528 : }
2529 5066 : if(ifile->FieldExist("lower_int")) {
2530 0 : ifile->scanField("lower_int",dummy);
2531 : }
2532 5066 : if(ifile->FieldExist("upper_int")) {
2533 0 : ifile->scanField("upper_int",dummy);
2534 : }
2535 2533 : ifile->scanField();
2536 : return true;
2537 : } else {
2538 : return false;
2539 : }
2540 : }
2541 :
2542 102 : void MetaD::computeReweightingFactor() {
2543 102 : if(biasf_==1.0) { // in this case we have no bias, so reweight factor is 0.0
2544 0 : getPntrToComponent("rct")->set(0.0);
2545 0 : return;
2546 : }
2547 :
2548 102 : double Z_0=0; //proportional to the integral of exp(-beta*F)
2549 102 : double Z_V=0; //proportional to the integral of exp(-beta*(F+V))
2550 102 : double minusBetaF=biasf_/(biasf_-1.)/kbt_;
2551 102 : double minusBetaFplusV=1./(biasf_-1.)/kbt_;
2552 102 : if (biasf_==-1.0) { //non well-tempered case
2553 0 : minusBetaF=1./kbt_;
2554 : minusBetaFplusV=0;
2555 : }
2556 102 : max_bias_=BiasGrid_->getMaxValue(); //to avoid exp overflow
2557 :
2558 102 : const unsigned rank=comm.Get_rank();
2559 102 : const unsigned stride=comm.Get_size();
2560 920504 : for (Grid::index_t t=rank; t<BiasGrid_->getSize(); t+=stride) {
2561 920402 : const double val=BiasGrid_->getValue(t);
2562 920402 : Z_0+=std::exp(minusBetaF*(val-max_bias_));
2563 920402 : Z_V+=std::exp(minusBetaFplusV*(val-max_bias_));
2564 : }
2565 102 : comm.Sum(Z_0);
2566 102 : comm.Sum(Z_V);
2567 :
2568 102 : reweight_factor_=kbt_*std::log(Z_0/Z_V)+max_bias_;
2569 204 : getPntrToComponent("rct")->set(reweight_factor_);
2570 : }
2571 :
2572 273 : double MetaD::getTransitionBarrierBias() {
2573 : // If there is only one well of interest, return the bias at that well point.
2574 273 : if (transitionwells_.size() == 1) {
2575 0 : double tb_bias = getBias(transitionwells_[0]);
2576 0 : return tb_bias;
2577 :
2578 : // Otherwise, check for the least barrier bias between all pairs of wells.
2579 : // Note that because the paths can be considered edges between the wells' nodes
2580 : // to make a graph and the path barriers satisfy certain cycle inequalities, it
2581 : // is sufficient to look at paths corresponding to a minimal spanning tree of the
2582 : // overall graph rather than examining every edge in the graph.
2583 : // For simplicity, I chose the star graph with center well 0 as the spanning tree.
2584 : // It is most efficient to start the path searches from the wells that are
2585 : // expected to be sampled last, so transitionwell_[0] should correspond to the
2586 : // starting well. With this choice the searches will terminate in one step until
2587 : // transitionwell_[1] is sampled.
2588 : } else {
2589 : double least_transition_bias;
2590 273 : std::vector<double> sink = transitionwells_[0];
2591 273 : std::vector<double> source = transitionwells_[1];
2592 273 : least_transition_bias = BiasGrid_->findMaximalPathMinimum(source, sink);
2593 273 : for (unsigned i = 2; i < transitionwells_.size(); i++) {
2594 0 : if (least_transition_bias == 0.0) {
2595 : break;
2596 : }
2597 0 : source = transitionwells_[i];
2598 0 : double curr_transition_bias = BiasGrid_->findMaximalPathMinimum(source, sink);
2599 0 : least_transition_bias = fmin(curr_transition_bias, least_transition_bias);
2600 : }
2601 : return least_transition_bias;
2602 : }
2603 : }
2604 :
2605 154 : void MetaD::updateFrequencyAdaptiveStride() {
2606 154 : plumed_massert(freq_adaptive_,"should only be used if frequency adaptive metadynamics is enabled");
2607 154 : plumed_massert(acceleration_,"frequency adaptive metadynamics can only be used if the acceleration factor is calculated");
2608 154 : const double mean_acc = acc_/((double) getStep());
2609 154 : int tmp_stride= stride_*floor((mean_acc/fa_min_acceleration_)+0.5);
2610 154 : if(mean_acc >= fa_min_acceleration_) {
2611 129 : if(tmp_stride > current_stride_) {
2612 6 : current_stride_ = tmp_stride;
2613 : }
2614 : }
2615 154 : if(fa_max_stride_!=0 && current_stride_>fa_max_stride_) {
2616 0 : current_stride_=fa_max_stride_;
2617 : }
2618 154 : getPntrToComponent("pace")->set(current_stride_);
2619 154 : }
2620 :
2621 8435 : bool MetaD::checkNeedsGradients()const {
2622 8435 : if(adaptive_==FlexibleBin::geometry) {
2623 192 : if(getStep()%stride_==0 && !isFirstStep_) {
2624 : return true;
2625 : } else {
2626 109 : return false;
2627 : }
2628 : } else {
2629 : return false;
2630 : }
2631 : }
2632 :
2633 4 : void MetaD::updateNlist() {
2634 : // no need to check for neighbors
2635 4 : if(hills_.size()==0) {
2636 0 : return;
2637 : }
2638 :
2639 : // here we generate the neighbor list
2640 4 : nlist_hills_.clear();
2641 : std::vector<Gaussian> local_flat_nl;
2642 4 : unsigned nt=OpenMP::getNumThreads();
2643 4 : if(hills_.size()<2*nt) {
2644 : nt=1;
2645 : }
2646 4 : #pragma omp parallel num_threads(nt)
2647 : {
2648 : std::vector<Gaussian> private_flat_nl;
2649 : #pragma omp for nowait
2650 : for(unsigned k=0; k<hills_.size(); k++) {
2651 : double dist2=0;
2652 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
2653 : const double d=difference(i,getArgument(i),hills_[k].center[i])/hills_[k].sigma[i];
2654 : dist2+=d*d;
2655 : }
2656 : if(dist2<=nlist_param_[0]*dp2cutoff) {
2657 : private_flat_nl.push_back(hills_[k]);
2658 : }
2659 : }
2660 : #pragma omp critical
2661 : local_flat_nl.insert(local_flat_nl.end(), private_flat_nl.begin(), private_flat_nl.end());
2662 : }
2663 4 : nlist_hills_ = local_flat_nl;
2664 :
2665 : // here we set some properties that are used to decide when to update it again
2666 12 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
2667 8 : nlist_center_[i]=getArgument(i);
2668 : }
2669 : std::vector<double> dev2;
2670 4 : dev2.resize(getNumberOfArguments(),0);
2671 46 : for(unsigned k=0; k<nlist_hills_.size(); k++) {
2672 126 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
2673 84 : const double d=difference(i,getArgument(i),nlist_hills_[k].center[i]);
2674 84 : dev2[i]+=d*d;
2675 : }
2676 : }
2677 12 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
2678 8 : if(dev2[i]>0.) {
2679 8 : nlist_dev2_[i]=dev2[i]/static_cast<double>(nlist_hills_.size());
2680 : } else {
2681 0 : nlist_dev2_[i]=hills_.back().sigma[i]*hills_.back().sigma[i];
2682 : }
2683 : }
2684 :
2685 : // we are done
2686 4 : getPntrToComponent("nlker")->set(nlist_hills_.size());
2687 4 : getPntrToComponent("nlsteps")->set(nlist_steps_);
2688 4 : nlist_steps_=0;
2689 4 : nlist_update_=false;
2690 4 : }
2691 :
2692 : }
2693 : }
|
