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 "MetaDKernel.h"
24 : #include "MetaDGrid.h"
25 : #include "MetaDHills.h"
26 : #include "MetaDRestart.h"
27 : #include "core/ActionRegister.h"
28 : #include "core/ActionSet.h"
29 : #include "core/PlumedMain.h"
30 : #include "core/FlexibleBin.h"
31 : #include "tools/Exception.h"
32 : #include "tools/Grid.h"
33 : #include "tools/OpenMP.h"
34 : #include "tools/File.h"
35 : #include "tools/Communicator.h"
36 : #include <algorithm>
37 : #include <cmath>
38 : #include <limits>
39 :
40 : namespace PLMD {
41 : namespace bias {
42 :
43 : //+PLUMEDOC BIAS METAD
44 : /*
45 : Used to perform 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 as discussed in the first paper cited below.
49 :
50 : $$
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 : $$
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 to 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 : $$
65 : V(\vec{s}) = -F(\vec{s})
66 : $$
67 :
68 : During post processing the free energy can be calculated in this way using the [sum_hills](sum_hills.md)
69 : utility.
70 :
71 : The following example provides an example input for a simple metadynamics calculation:
72 :
73 : ```plumed
74 : d1: DISTANCE ATOMS=3,5
75 : d2: DISTANCE ATOMS=2,4
76 : restraint: METAD ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3 PACE=500
77 : PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
78 : ```
79 :
80 : In this calculation the distance between atoms 3 and 5
81 : and the distance between atoms 2 and 4 are used as collective variables. Gaussian hills with a height of 0.3 kJ/mol
82 : and widths of 0.2 nm in the two CVs are added every 500 steps. The values
83 : of the CVs and the metadynamics bias potential are written to the COLVAR file every 100 steps.
84 :
85 : ## Using grids
86 :
87 : In the simplest possible implementation of a metadynamics calculation __the expense associated with calculating the
88 : bias potential increases with the length of the simulation as one has to, at every step, evaluate
89 : the values of a larger and larger number of Gaussian kernels. To avoid this issue you can
90 : store the bias on a grid.__ This approach is similar to that proposed in the second paper cited below but has the
91 : advantage that the grid spacing is independent of the Gaussian width. An example input for a calculation
92 : that uses a grid to store the bias is shown below:
93 :
94 :
95 : ```plumed
96 : t1: TORSION ATOMS=1,2,3,4
97 : t2: TORSION ATOMS=5,6,7,8
98 : m: METAD ARG=t1,t2 SIGMA=0.1,0.1 HEIGHT=1.2 PACE=500 GRID_MIN=-pi,-pi GRID_MAX=pi,pi GRID_BIN=100,100
99 : PRINT ARG=t1,t2,m.bias STRIDE=100 FILE=COLVAR
100 : ```
101 :
102 : GRID_MIN and GRID_MAX give the initial range of the grid. You can also specify the number
103 : of grid points for every collective variable (GRID_BIN) or the desired grid spacing (GRID_SPACING).
104 : If you specify the GRID_SPACING and the number of bins PLUMED uses the most conservative choice. In other words, it uses
105 : the choice that gives you the largest number of bins for each dimension.
106 :
107 : If you specify neither GRID_BIN nor GRID_SPACING PLUMED uses 1/5 of the Gaussian width (SIGMA) as grid spacing if the width is
108 : fixed or 1/5 of the minimum Gaussian width (SIGMA_MIN - see next section) if the width is variable.
109 : This default choice should be reasonable for most applications.
110 :
111 : The grid is extended automatically, along the variables that are not periodic, whenever a Gaussian does not fit in it:
112 : the grid always contains all the points where a Gaussian is not zero (Gaussians are truncated at a distance of
113 : $\sqrt{2\cdot 6.25}$ widths, and stretched so that they go to zero there), so that the
114 : bias is exactly zero outside it (with INTERVAL, the grid contains the interval, and beyond the limits the bias is
115 : the constant value at the limits). The grid spacing does not change and the new points are on the same lattice as the
116 : initial ones, so the bias is the same as with a grid that is large enough from the beginning. Every extension is
117 : written on the log. If you give GRID_SPACING without GRID_MIN and GRID_MAX the grid is placed around the first Gaussian and
118 : then extended as needed, so that you do not need to know in advance the range explored by the variables. The variables
119 : that are periodic are always covered by the grid over the whole period. GRID_SIZE_LIMIT sets the maximum number of points
120 : of the grid: the calculation stops with an error if a Gaussian would require a larger grid (e.g. with a distance that grows
121 : without limit); notice that every MPI process has its own copy of the grid. Quantities that are obtained from the whole
122 : grid (such as the c(t) reweighting factor, see CALC_RCT) are computed over the current range of the grid, so they change
123 : when the grid is extended. A TARGET grid is not extended: the variables must remain in its range.
124 :
125 : ```plumed
126 : d: DISTANCE ATOMS=1,10
127 : m: METAD ARG=d SIGMA=0.05 HEIGHT=1.2 PACE=500 GRID_SPACING=0.01
128 : PRINT ARG=d,m.bias STRIDE=100 FILE=COLVAR
129 : ```
130 :
131 : ## Using neighbour lists
132 :
133 : Another way to decrease the cost associated with evaluating the metadynamics bias is to use a neighbor list. If you
134 : wish to use this alternative to grids you use the NLIST keyword as illustrated below.
135 :
136 : ```plumed
137 : d1: DISTANCE ATOMS=1,2
138 : d2: DISTANCE ATOMS=3,4
139 : d3: DISTANCE ATOMS=5,6
140 : d4: DISTANCE ATOMS=7,8
141 : m: METAD ARG=d1,d2,d3,d4 SIGMA=0.1,0.1,0.1,0.1 HEIGHT=1.2 PACE=500 NLIST
142 : PRINT ARG=d1,d2,d3,d4,m.bias,m.nlker FILE=colvar
143 : ```
144 :
145 : Using NLIST is beneficial if you are using more than 2 collective variables as if you are using large numbers of CVs GRID becomes
146 : expensive and memory consuming. Furthermore, as indicated above, if you are using this option you can have access to a value called
147 : `nlker`, which tells you how many Gaussians were evaluated in each step.
148 :
149 : The neighbor list is updated everytime the CVs move farther than a cut-off value
150 : from the position they were at last neighbor list update. Gaussians are added to the neigbhor list if their center
151 : is within 6.*DP2CUTOFF*sigma*sigma. While the list is updated if the CVs are farther from the center than 0.5 of the
152 : standard deviation of the Gaussian center distribution of the list. These parameters (6 and 0.5) can be modified using
153 : NLIST_PARAMETERS. Note that the use of neighbor list does not provide the exact bias.
154 :
155 : ## Restarting a metadynamics simulation
156 :
157 : Metadynamics can be restarted from a HILLS file as illustrated in the following input:
158 :
159 : ```plumed
160 : #SETTINGS INPUTFILES=regtest/basic/rt-mpi6/HILLS_multi
161 :
162 : RESTART
163 : d: DISTANCE ATOMS=1,10
164 : c: COORDINATION GROUPA=1-108 GROUPB=1-108 R_0=0.5
165 : m: METAD ARG=c,d SIGMA=0.1,0.2 HEIGHT=0.1 PACE=500 FILE=regtest/basic/rt-mpi6/HILLS_multi
166 : PRINT ARG=d,c FILE=colvar STRIDE=10
167 : ```
168 :
169 : If you don't specify the location of the HILLS file then PLUMED looks for a file called HILLS to read in.
170 : Metadynamics can also be restarted from a GRID, as in this second example.
171 :
172 : ```plumed
173 : #SETTINGS INPUTFILES=extras/1d_bias.grid
174 :
175 : d1: DISTANCE ATOMS=1,2
176 : m: METAD ...
177 : ARG=d1 SIGMA=0.1
178 : HEIGHT=0.1 PACE=500
179 : GRID_MIN=1.14 GRID_MAX=1.32 GRID_BIN=5
180 : GRID_RFILE=extras/1d_bias.grid
181 : RESTART=YES
182 : ...
183 : PRINT ARG=d1,m.bias FILE=colvar STRIDE=10
184 : ```
185 :
186 : Notice that if you want to use this option you first need to save the GRID using GRID_WFILE (and GRID_WSTRIDE) in the
187 : previous simulation.
188 :
189 : No Gaussian is added at the first step of a simulation. When a simulation is continued with RESTART, a Gaussian is added
190 : at the first step unless the hills file already contains it, i.e. unless the previous run had already performed this step
191 : (some MD codes compute again the last step of the previous run when restarting, others do not). When the bias is read
192 : from a grid, the hills file (FILE) is still read, without adding its hills to the bias, to know the time of its last hill.
193 :
194 : Notice, also that in the first input above the
195 : [RESTART](RESTART.md) action was used, while the second input used METAD's RESTART keyword. If you use the first input
196 : every action in the input is restarted so the [PRINT](PRINT.md) command will append to the output file from the earlier calculation.
197 : In this second input, however, only the METAD action is restarted so the [PRINT](PRINT.md) file will output data to a new file and
198 : back up any old files it finds with the same name.
199 :
200 : ## Calculating the work done by the bias
201 :
202 : The work performed by the METAD bias can be calculated by using the CALC_WORK option. You can then output the work as shown below:
203 :
204 : ```plumed
205 : mu1: DISTANCE ATOMS=1,10
206 : vol: VOLUME
207 :
208 : md: METAD ARG=mu1,vol SIGMA=0.1,0.2 HEIGHT=1.0 PACE=10 GRID_MIN=-5,20 GRID_MAX=5,40 GRID_BIN=100,100 CALC_WORK
209 :
210 : PRINT ARG=mu1,vol,md.bias,md.work FILE=colvar
211 : ```
212 :
213 : Note that this calculation is expensive when not using grids.
214 :
215 : ## Well-tempered metadynamics
216 :
217 : Another option that is available in plumed is well-tempered metadynamics that is introduced in the third paper cited below. In this
218 : variant of metadynamics the heights of the Gaussian hills are scaled at each step so the bias is now
219 : given by:
220 :
221 : $$
222 : 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(
223 : -\sum_{i=1}^{d} \frac{(s_i({q})-s_i({q}(t')))^2}{2\sigma_i^2}
224 : \right),
225 : $$
226 :
227 : This method ensures that the bias converges more smoothly. The following example shows an example input for a well-tempered metadynamics
228 : calculation.
229 :
230 : ```plumed
231 : d1: DISTANCE ATOMS=3,5
232 : d2: DISTANCE ATOMS=2,4
233 : restraint: METAD ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3 BIASFACTOR=10 PACE=500
234 : PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
235 : ```
236 :
237 : It should be noted that, if you do well-tempered metadynamics,
238 : the hills file output contains the value of Gaussian height re-scaled by the bias factor.
239 : This rescaling is done to ensure that the HILLS file contains
240 : negative of the free-energy estimate and __not__ the bias potential. An added bonus of this choice is that
241 : one can restart a simulation using a different value for the $\Delta T$. The applied bias will be scaled accordingly.
242 :
243 : ## Metadynamics with flexible Gaussians
244 :
245 : You can also use the flexible Gaussian approach that is discussed in the fourth paper cited below.
246 : This method adapts the Gaussian to the extent of Cartesian space covered by a variable or
247 : to the space in collective variable covered in a given time. When you use these methods the width of the deposited
248 : Gaussian potential is denoted by one value only that is a Cartesian space (ADAPTIVE=GEOM) or a time
249 : (ADAPTIVE=DIFF). Note that a specific integration technique, which is discussed in the documentation for
250 : [sum_hills](sum_hills.md), must be used when calculating the free energy from the deposited Gaussian kernels in this case.
251 :
252 : The following example shows how to use adaptive Gaussian kernels, with the diffusion scheme. In this example
253 : Gaussians that cover the space of 20 time steps in collective variables are used.
254 :
255 : ```plumed
256 : d1: DISTANCE ATOMS=3,5
257 : d2: DISTANCE ATOMS=2,4
258 : restraint: METAD ARG=d1,d2 SIGMA=20 HEIGHT=0.3 PACE=500 ADAPTIVE=DIFF
259 : PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
260 : ```
261 :
262 : This second example shows you how to use adaptive Gaussian kernels, with the geometrical scheme. In this example
263 : Gaussians that cover the space of 0.05 nm in Cartesian space are used
264 :
265 : ```plumed
266 : d1: DISTANCE ATOMS=3,5
267 : d2: DISTANCE ATOMS=2,4
268 : restraint: METAD ARG=d1,d2 SIGMA=0.05 HEIGHT=0.3 PACE=500 ADAPTIVE=GEOM
269 : PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
270 : ```
271 :
272 : With both these methods you can limit how much the hills width can change by using the SIGMA_MIN and SIGMA_MAX keywords
273 : as shown below.
274 :
275 : ```plumed
276 : d1: DISTANCE ATOMS=3,5
277 : d2: DISTANCE ATOMS=2,4
278 : restraint: METAD ...
279 : ARG=d1,d2 SIGMA=0.05 HEIGHT=0.3 PACE=500 ADAPTIVE=GEOM
280 : SIGMA_MIN=0.2,0.1 SIGMA_MAX=0.5,1.0
281 : ...
282 : PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
283 : ```
284 :
285 : The sigma values that are provided to these two keywords are specified in the units of the CV.
286 : Note that if you use a negative number, this means that the limit is not set.
287 :
288 : With ADAPTIVE=DIFF the width of the hills is obtained from running averages of the CV values and of
289 : their fluctuations over a time given by SIGMA (in time steps). These statistics are accumulated when the bias
290 : is updated, so that with UPDATE_FROM they start only at the beginning of the update window and the first
291 : hills can be much narrower than the following ones. With the PRECAL_DIFF flag the statistics are
292 : accumulated from the beginning of the simulation: by using UPDATE_FROM with a time
293 : larger than the one given by SIGMA, all the hills have a width estimated over a sufficiently long time.
294 :
295 : ```plumed
296 : d1: DISTANCE ATOMS=3,5
297 : restraint: METAD ARG=d1 SIGMA=500 HEIGHT=0.3 PACE=50 ADAPTIVE=DIFF PRECAL_DIFF UPDATE_FROM=2.0
298 : ```
299 :
300 : When restarting, the statistics are initialized from the last hills of the hills file of the walker:
301 : the covariance of the fluctuations is the one of the last hill, and the running average of the CVs is
302 : estimated from the centers of the last hills, with the weights used by the running average.
303 : This is also done when the bias is read from a grid (GRID_RFILE), if the hills file is available,
304 : using only the hills deposited before the current time.
305 : In this way the hills deposited after a restart have the same width they would have had without the restart,
306 : apart from the effect of the limits set with SIGMA_MIN and SIGMA_MAX, which are applied to the hills.
307 : Only the last hills of the walker that belong to the same run and precede the current time are used.
308 : The hills file does not record how the hills were obtained: if it was written with a different
309 : SIGMA, or with ADAPTIVE=GEOM, use the DIFF_FROM_SCRATCH flag. With this flag the statistics are
310 : accumulated again from the restart, as in previous versions; no hill is then added before at least two
311 : CV values are available.
312 :
313 : ADAPTIVE=DIFF can be used with replica exchange. Since the MD code does not tell PLUMED
314 : whether an exchange has been accepted, after an exchange attempt the CV values are compared with those
315 : of the exchange partner: when the exchange is detected as accepted, the running average of the CVs is
316 : restarted from their new values, while the estimate of their fluctuations is kept.
317 : Note also that in this case the histogram correction is needed when summing up hills.
318 :
319 : ## Metadynamics with the interval keyword
320 :
321 : When you use the keyword INTERVAL the metadynamics bias force is set equal to zero
322 : outside of a boundary as discussed in the fifth paper cited below. This feature is useful if, for example, metadynamics is performed on one or more CVs and if one is only interested
323 : in the free energy for s > boundary. When you use this option the history dependent potential is still updated in accordance with the formulas above.
324 : The difference is simply that the metadynamics force is set to zero for s < boundary. This means that Gaussian kernels are added also
325 : if s < boundary, as the tails of these Gaussian kernels influence VG in the relevant region s > boundary. In this way, the
326 : force on the system in the region s > boundary comes from both metadynamics and the force field, while in the region
327 : s < boundary the only force the system experiences is from the forcefield. This approach allows one to obtain a history-dependent bias potential VG that
328 : fluctuates around a stable estimator, equal to the negative of the free energy far enough from the
329 : boundaries. Notice that:
330 :
331 : - This method only works for one-dimensional biases;
332 : - It works both with and without GRID;
333 : - The interval limit boundary should be set in a region where the derivative of the free energy is not large;
334 : - If in the region outside the limit boundary the system has a free energy minimum, the INTERVAL keyword should
335 : be used together with a [UPPER_WALLS](UPPER_WALLS.md) or [LOWER_WALLS](LOWER_WALLS.md) that act when the system is beyond the boundary as illustrated in the following example.
336 :
337 : ```plumed
338 : r2g: DISTANCE ATOMS=1,2
339 : m2g: METAD ARG=r2g HEIGHT=0.1 PACE=500 SIGMA=1.0 INTERVAL=4,8 FILE=H2G GRID_MIN=0 GRID_MAX=20 GRID_BIN=100
340 : lw: LOWER_WALLS ARG=r2g AT=4 KAPPA=100
341 : uw: UPPER_WALLS ARG=r2g AT=8 KAPPA=100
342 : PRINT ARG=r2g,*.bias FILE=colvar STRIDE=10
343 : ```
344 :
345 : As a final note, since version 2.0.2 the force is set to zero and the bias value to the value at the corresponding boundary
346 : when the system is outside of the selected interval. This choice allows acceptances for replica exchange methods to be computed correctly.
347 :
348 : ## Multiple walkers metadynamics
349 :
350 : The multiple walkers metadynamics method that is described in the sixth paper cited below can be used by using an input
351 : like the one shown below:
352 :
353 : ```plumed
354 : d1: DISTANCE ATOMS=3,5
355 : restraint: METAD ...
356 : ARG=d1 SIGMA=0.05 HEIGHT=0.3 PACE=500
357 : WALKERS_N=10
358 : WALKERS_ID=3
359 : WALKERS_DIR=../
360 : WALKERS_RSTRIDE=100
361 : ...
362 : ```
363 :
364 : To use this method you need to set the number of walkers (WALKERS_N), the ID of the
365 : current walker (WALKERS-ID), the directory where the files containing the hills files reside
366 : (WALKERS_DIR) and the frequency to read the other walkers (WALKERS_RSTRIDE).
367 : The bias potential experienced by all the walkers is synched between walkers every
368 : WALKERS_RSTRIDE steps. In addition, since version 2.2.5, hills files are automatically
369 : flushed every WALKERS_RSTRIDE steps.
370 : The grid written by a walker (GRID_WFILE) contains also the hills it has read from the files of the other walkers: their
371 : number is written in a comment line after the grid. When restarting from this grid (GRID_RFILE), only the hills of the
372 : other walkers that are not in it are added. For a grid written by an older version of PLUMED, the grid is assumed to
373 : contain the hills of the other walkers up to the time of the last hill in the hills file of this walker.
374 : With multiple walkers, the hills files and the grid files are read and written without the suffix of the replica
375 : (e.g. with `--multi`): give a different GRID_WFILE to each walker (e.g. with `@replicas:`).
376 :
377 : ## Reweighting a metadynamics simulation
378 :
379 : The $c(t)$ reweighting factor can also be calculated on the fly using the equations
380 : presented in the seventh paper in the references section by using the keyword CALC_RCT as illustrated in the example input below:
381 :
382 : ```plumed
383 : phi: TORSION ATOMS=1,2,3,4
384 : psi: TORSION ATOMS=5,6,7,8
385 :
386 : m: METAD ...
387 : ARG=phi,psi SIGMA=0.20,0.20 HEIGHT=1.20 BIASFACTOR=5 TEMP=300.0 PACE=500
388 : GRID_MIN=-pi,-pi GRID_MAX=pi,pi GRID_BIN=150,150
389 : CALC_RCT
390 : RCT_USTRIDE=10
391 : ...
392 :
393 : rw: REWEIGHT_METAD TEMP=300
394 :
395 : h: HISTOGRAM ...
396 : ARG=phi,psi
397 : GRID_MIN=-pi,-pi GRID_MAX=pi,pi GRID_BIN=150,150
398 : BANDWIDTH=0.1,0.1
399 : LOGWEIGHTS=rw
400 : NORMALIZATION=true
401 : ...
402 :
403 : ff: CONVERT_TO_FES ARG=h TEMP=300
404 : DUMPGRID ARG=ff FILE=fes.dat
405 : ```
406 :
407 :
408 : The expression used to calculate $c(t)$ follows directly from Eq. 3 in the seventh paper cited below,
409 : where $F(\vec{s})=-\gamma/(\gamma-1) V(\vec{s})$.
410 : This gives smoother results than Eqs. 13 and Eqs. 14 in that paper.
411 :
412 : Notice that you can only use this option if you are accumulating the metadynamics bias on a grid.
413 : $c(t)$ is given by the rct component of the METAD action above while the bias
414 : normalized by $c(t)$ is given by the rbias component (rbias=bias-rct). As illustrated above this
415 : value can be used to obtain a reweighted histogram.
416 :
417 : By default $c(t)$ is updated every time the bias changes. If this slows down the simulation, you can use the keyword
418 : RCT_USTRIDE to update it only every RCT_USTRIDE*PACE steps (every RCT_USTRIDE Gaussian depositions, unless the pace
419 : changes with FREQUENCY_ADAPTIVE): in the example input above $c(t)$ is calculated every 10 Gaussian depositions.
420 :
421 :
422 : ## Concurrent metadynamics
423 :
424 : To implement the method of concurrent metadynamics that is discussed in the eighth paper cited below you insert multiple
425 : METAD actions in your plumed input file.
426 :
427 : ## Infrequent metadynamics
428 :
429 : The kinetics of the transitions between basins can also be analyzed on the fly as
430 : in discussed in the 9th paper cited below. The flag ACCELERATION turn on accumulation of the acceleration
431 : factor that can then be used to determine the rate as illustrated by the following input.
432 :
433 : ```plumed
434 : psi: ANGLE ATOMS=7,9,15
435 :
436 : md1: METAD ...
437 : ARG=psi SIGMA=0.20
438 : HEIGHT=1.20 TEMP=300.0 BIASFACTOR=10
439 : PACE=500 GRID_MIN=0 GRID_MAX=pi GRID_BIN=10
440 : ACCELERATION
441 : ...
442 :
443 : PRINT ARG=psi,md1.bias,md1.acc FILE=colvar STRIDE=10
444 : ```
445 :
446 : This method can be used together with [COMMITTOR](COMMITTOR.md) analysis to stop the simulation when the system get to the target basin.
447 : It is meant to be used together with Well-Tempered Metadynamics. If restarting from a previous
448 : metadynamics you need to use the ACCELERATION_RFILE keyword to give the name of the
449 : data file from which the previous value of the acceleration factor should be read, otherwise the
450 : calculation of the acceleration factor will be wrong.
451 :
452 : ## Changing the deposition frequency adaptively
453 :
454 : By using the flag FREQUENCY_ADAPTIVE the frequency adaptive scheme introduced in the 10th paper cited below
455 : is turned on. The frequency for hill addition then changes dynamically based on the acceleration factor
456 : according to the following equation
457 :
458 : $$
459 : \tau_{\mathrm{dep}}(t) =
460 : \min\left[
461 : \tau_0 \cdot
462 : \max\left[\frac{\alpha(t)}{\theta},1\right]
463 : ,\tau_{c}
464 : \right]
465 : $$
466 :
467 : where $\tau_0$ is the initial hill addition frequency given by the PACE keyword,
468 : $\tau_{c}$ is the maximum allowed frequency given by the FA_MAX_PACE keyword,
469 : $\alpha(t)$ is the instantaneous acceleration factor at time $t$,
470 : and $\theta$ is a threshold value that acceleration factor has to reach before
471 : triggering a change in the hill addition frequency given by the FA_MIN_ACCELERATION keyword.
472 : The frequency for updating the hill addition frequency according to this equation is
473 : given by the FA_UPDATE_FREQUENCY keyword, by default it is the same as the value given
474 : in PACE. The hill addition frequency increases monotonically such that if the
475 : instantaneous acceleration factor is lower than in the previous updating step the
476 : previous $\tau_{\mathrm{dep}}$ is kept rather than updating it to a lower value.
477 : The instantaneous hill addition frequency $\tau_{\mathrm{dep}}(t)$ is outputted
478 : to pace component. Note that if restarting from a previous metadynamics run you need to
479 : use the ACCELERATION_RFILE keyword to read in the acceleration factors from the
480 : previous run (the acceleration factor and the pace, if they are in the file), otherwise the hill addition frequency
481 : will start from the initial frequency.
482 :
483 : ## Targeted metadynamics
484 :
485 : You can also provide a target distribution for your metadynamics simulations by using the keyword TARGET.
486 : This idea is discussed in the last four papers cited in the references section.
487 :
488 : The TARGET should be a grid containing a free-energy (i.e. the -$k_B$T*log for the desired target distribution).
489 : Gaussian kernels will then be scaled by a factor
490 :
491 : $$
492 : e^{\beta(\tilde{F}(s)-\tilde{F}_{max})}
493 : $$
494 :
495 : Here $\tilde{F}(s)$ is the free energy defined on the grid and $\tilde{F}_{max}$ its maximum value.
496 : Notice that we here used the maximum value as was done in the penultimate paper in the reference section.
497 : This choice allows one to avoid the addition of exceedingly large Gaussian kernels. However,
498 : it can also make the Gaussian too small. You should always choose the HEIGHT parameter carefully
499 : in this case.
500 : The grid file should be similar to other PLUMED grid files in that it should contain
501 : both the target free-energy and its derivatives.
502 :
503 : Notice that if you wish your simulation to converge to the target free energy you should use
504 : the DAMPFACTOR command to provide a global tempering as is discussed in the last paper from the reference section.
505 : Alternatively, if you use a BIASFACTOR your simulation will converge to a free
506 : energy that is a linear combination of the target free energy and of the intrinsic free energy
507 : determined by the original force field.
508 :
509 : ```plumed
510 : #SETTINGS INPUTFILES=extras/target.fes
511 :
512 : d1: DISTANCE ATOMS=3,5
513 : t1: METAD ...
514 : ARG=d1 SIGMA=0.05 TAU=200 DAMPFACTOR=100 PACE=250
515 : GRID_MIN=1.14 GRID_MAX=1.32 GRID_BIN=6
516 : TARGET=extras/target.fes
517 : ...
518 :
519 : PRINT ARG=d1,t1.bias STRIDE=100 FILE=COLVAR
520 : ```
521 :
522 : Notice that the BIASFACTOR can also be set equal to 1. In this case one will perform
523 : unbiased sampling: no bias is applied, and the hills files and the grid contain the estimate of the negative of the
524 : free energy (the heights of the hills are those of the limit of well-tempered metadynamics for a bias factor going to 1,
525 : as in [PBMETAD](PBMETAD.md)). Instead of using HEIGHT, one should provide the TAU parameter as illustrated below.
526 :
527 : ```plumed
528 : d: DISTANCE ATOMS=3,5
529 : METAD ARG=d SIGMA=0.1 TAU=4.0 TEMP=300 PACE=100 BIASFACTOR=1.0
530 : ```
531 :
532 : The HILLS file obtained will still work with `plumed sum_hills` so as to plot a free-energy.
533 : The case where this makes sense is probably that of RECT simulations.
534 :
535 : Regarding RECT simulations, you can also use the RECT keyword so as to avoid using multiple input files.
536 : For instance, a single input file will be
537 :
538 : ```plumed
539 : d: DISTANCE ATOMS=3,5
540 : METAD ARG=d SIGMA=0.1 TAU=4.0 TEMP=300 PACE=100 RECT=1.0,1.5,2.0,3.0
541 : ```
542 :
543 : The number of elements in the RECT array should be equal to the number of replicas.
544 :
545 : */
546 : //+ENDPLUMEDOC
547 :
548 : class MetaD : public Bias {
549 :
550 : private:
551 : using Gaussian = MetaDGaussian;
552 15 : struct TemperingSpecs {
553 : bool is_active;
554 : std::string name_stem;
555 : std::string name;
556 : double biasf;
557 : double threshold;
558 : double alpha;
559 : //no need for a specialized ctor: data is assigned in order of declaration (but needs an initializer list)
560 : };
561 : // general setup
562 : double kbt_=0.0;
563 : int stride_=0;
564 : /// current hill deposition pace (equal to stride_ unless FREQUENCY_ADAPTIVE is used)
565 : int current_stride_=0;
566 : bool calc_work_=false;
567 : // well-tempered MetaD
568 : bool welltemp_=false;
569 : double biasf_=-1.0;
570 : // output files format
571 : std::string fmt_;
572 : /// the steps when hills are added (the first step of the run, the last hill read when restarting)
573 : MetaDHillSchedule schedule_;
574 : // Gaussian starting parameters
575 : double height0_=std::numeric_limits<double>::max();
576 : std::vector<double> sigma0_;
577 : std::vector<double> sigma0min_;
578 : std::vector<double> sigma0max_;
579 : // Gaussians
580 : std::vector<Gaussian> hills_;
581 : /// the hills without grid also as contiguous arrays (center and inverse widths of diagonal hills),
582 : /// for the loops that sum them
583 : struct HillArrays {
584 : std::vector<double> center;
585 : std::vector<double> invsigma;
586 : std::vector<double> height;
587 : std::vector<char> diagonal;
588 54 : void clear() {
589 54 : center.clear();
590 54 : invsigma.clear();
591 54 : height.clear();
592 54 : diagonal.clear();
593 54 : }
594 : } hillarrays_;
595 : /// value of hill h, and its derivatives (subtracted to der); dp is a work array
596 : double evaluateHill(const std::vector<double>& cv, unsigned h) const;
597 : double evaluateHillWithDerivatives(const std::vector<double>& cv, unsigned h, double* der, double* dp) const;
598 : std::unique_ptr<FlexibleBin> flexbin_;
599 : int adaptive_=FlexibleBin::none;
600 : /// ADAPTIVE=DIFF: accumulate the CV statistics at every step, also outside the update window (PRECAL_DIFF)
601 : bool precal_diff_=false;
602 : /// ADAPTIVE=DIFF: when restarting, the statistics of the fluctuations are initialized from the last hills
603 : /// read from the hills file of this walker (unless DIFF_FROM_SCRATCH). They are kept here until the
604 : /// first update, when the current time is known.
605 : MetaDDiffusionRestart diffrestart_;
606 : /// ADAPTIVE=DIFF with replica exchange: CVs of the exchange partner and CVs at the previous update,
607 : /// used to detect whether an attempted exchange was accepted
608 270 : struct ExchangeDetection {
609 : bool attempted=false;
610 : std::vector<double> partner_cv;
611 : std::vector<double> last_cv;
612 : } exchange_;
613 : OFile hillsOfile_;
614 : /// hills files that are read (own file when restarting, files of the other walkers), and writing of the hills
615 : std::unique_ptr<MetaDHillsFiles> hillsfiles_;
616 : // Grids
617 : bool grid_=false;
618 : /// the bias grid, extended when a hill does not fit in it
619 : std::unique_ptr<MetaDBiasGrid> BiasGrid_;
620 : /// maximum number of points of the grid (GRID_SIZE_LIMIT)
621 : double gridsizelimit_=1.0e8;
622 : /// WALKERS_N, restart from a grid: the number of hills of each of the other walkers in the grid (-1 for this walker);
623 : /// empty if not known (grid files written by older versions)
624 : OFile gridfile_;
625 : bool storeOldGrids_=false;
626 : int wgridstride_=0;
627 : /// multiple walkers, sharing hills either through files (WALKERS_N) or through MPI (WALKERS_MPI)
628 : MetaDWalkers mw_;
629 : /// flying Gaussians (MPI walkers only)
630 : bool flying_=false;
631 : // kinetics from metadynamics
632 : bool acceleration_=false;
633 : double acc_=0.0;
634 : double acc_restart_mean_=0.0;
635 : /// the step of the acceleration factor read from ACCELERATION_RFILE
636 : long long int acc_restart_step_=0;
637 : // transition-tempering metadynamics
638 : bool calc_max_bias_=false;
639 : double max_bias_=0.0;
640 : bool calc_transition_bias_=false;
641 : double transition_bias_=0.0;
642 : std::vector<std::vector<double> > transitionwells_;
643 : static const size_t n_tempering_options_ = 1;
644 : static const std::string tempering_names_[1][2];
645 : double dampfactor_=0.0;
646 : TemperingSpecs tt_specs_{false,"TT","Transition Tempered", -1.0, 0.0, 1.0};
647 : std::string targetfilename_;
648 : std::unique_ptr<GridBase> TargetGrid_;
649 : /// frequency adaptive metadynamics
650 255 : struct FrequencyAdaptive {
651 : bool enabled=false;
652 : /// how often the pace is updated
653 : int update_frequency=0;
654 : /// maximum pace (0 means no maximum)
655 : int max_stride=0;
656 : /// the pace is changed only when the acceleration is larger than this
657 : double min_acceleration=1.0;
658 : } fa_;
659 : /// INTERVAL (1D only): outside [lower,upper] the bias force is not felt
660 255 : struct Interval {
661 : bool enabled=false;
662 : double lower=-1;
663 : double upper=-1;
664 : } interval_;
665 : // reweighting
666 : bool calc_rct_=false;
667 : double reweight_factor_=0.0;
668 : unsigned rct_ustride_=1;
669 : // work
670 : double work_=0;
671 : /// neighbor list of the Gaussians that are close to the current CV values (NLIST)
672 255 : struct NeighborList {
673 : bool enabled=false;
674 : /// the list has to be rebuilt
675 : bool update=false;
676 : /// steps since the last update
677 : unsigned steps=0;
678 : /// cutoff (in units of dp2cutoff) and rebuild threshold (in units of dev2)
679 : std::array<double,2> param;
680 : /// indices of the listed hills in hills_
681 : std::vector<unsigned> index;
682 : /// CV values at the last update
683 : std::vector<double> center;
684 : /// mean square distance of the listed hills from center, per CV
685 : std::vector<double> dev2;
686 : } nlist_;
687 :
688 : /// evaluation of the Gaussians
689 : MetaDKernel kernel_;
690 :
691 : /// values parsed in the constructor that are needed by more than one setup step
692 : struct SetupInput {
693 : std::string hillsfname;
694 : std::vector<double> rect_biasf;
695 : double tau=0.0;
696 : std::vector<std::string> gmin;
697 : std::vector<std::string> gmax;
698 : std::vector<unsigned> gbin;
699 : /// GRID_SPACING without GRID_MIN and GRID_MAX: the grid is placed around the first hill
700 : bool gridnobounds=false;
701 : std::vector<double> gspacing;
702 : bool sparsegrid=false;
703 : bool spline=true;
704 : std::string gridfilename;
705 : std::string gridreadfilename;
706 : std::string acc_rfilename;
707 : bool restartedFromGrid=false;
708 : /// WALKERS_N: the numbers of hills of the other walkers in the grid read from GRID_RFILE
709 : std::vector<long long> gridWalkerHills;
710 : };
711 : // setup steps called by the constructor, in this order
712 : /// ADAPTIVE, FMT, SIGMA, SIGMA_MIN, SIGMA_MAX
713 : void parseGaussianWidth();
714 : /// HEIGHT, PACE, FILE, RECT, BIASFACTOR, DAMPFACTOR, CALC_WORK, transition tempering, TARGET, TAU
715 : void parseHillsAndTempering(SetupInput& in);
716 : /// GRID_* keywords
717 : void parseGrid(SetupInput& in);
718 : /// NLIST, NLIST_PARAMETERS
719 : void parseNeighborList();
720 : /// CALC_RCT, RCT_USTRIDE
721 : void parseReweightingFactor();
722 : /// WALKERS_*, FLYING_GAUSSIAN
723 : void parseMultipleWalkers();
724 : /// INTERVAL
725 : void parseInterval();
726 : /// ACCELERATION, ACCELERATION_RFILE, FREQUENCY_ADAPTIVE and FA_*
727 : void parseAcceleration(SetupInput& in);
728 : /// log Gaussian parameters
729 : void logHillsSetup(const SetupInput& in);
730 : /// check consistency and stability of the tempering options
731 : void checkTempering();
732 : /// log interval and grid parameters
733 : void logIntervalAndGrid(const SetupInput& in);
734 : /// check and log multiple walkers options
735 : void setupMultipleWalkers();
736 : /// add output components (and read ACCELERATION_RFILE)
737 : void setupComponents(const SetupInput& in);
738 : /// create the bias grid, or read it from GRID_RFILE
739 : void setupBiasGrid(SetupInput& in);
740 : /// open the hills files of all walkers and read them when restarting
741 : void readHillsFiles(const SetupInput& in);
742 : /// read the TARGET grid
743 : void setupTargetGrid();
744 : /// initialize the quantities that depend on the bias when restarting
745 : /// the quantities that depend on the bias read when restarting (or from GRID_RFILE)
746 : void initializeRestartQuantities(const SetupInput& in);
747 : /// open the grid and hills files for writing
748 : void openOutputFiles(SetupInput& in);
749 : /// log the bibliography
750 : void logBibliography(const SetupInput& in);
751 :
752 : static void registerTemperingKeywords(const std::string &name_stem, const std::string &name, Keywords &keys);
753 : void readTemperingSpecs(TemperingSpecs &t_specs);
754 : void logTemperingSpecs(const TemperingSpecs &t_specs);
755 : /// ADAPTIVE=DIFF: initialize the statistics of the fluctuations from the hills kept in diffrestart_
756 : void initializeDiffusionStatistics();
757 : void writeGaussian(const Gaussian&,OFile&);
758 : void addGaussian(const Gaussian&);
759 : double getHeight(const std::vector<double>&);
760 : void temperHeight(double &height, const TemperingSpecs &t_specs, const double tempering_bias);
761 : double getBias(const std::vector<double>&);
762 : double getBiasAndDerivatives(const std::vector<double>&, std::vector<double>&);
763 : /// update the adaptive hill widths; with ADAPTIVE=DIFF also handle replica exchanges
764 : void updateFlexibleBin(bool nowAddAHill);
765 : /// true if a hill has to be added at the current step
766 : bool isHillStep() const;
767 : /// the c(t) reweighting factor (CALC_RCT)
768 : void computeReweightingFactor();
769 : double getTransitionBarrierBias();
770 : /// add a hill in cv (and the hills of the other MPI walkers), write them to the hills file
771 : void addHill(const std::vector<double>& cv, bool multivariate);
772 : /// write the bias grid to GRID_WFILE
773 : void writeBiasGrid();
774 : /// read the new hills deposited by the other walkers (WALKERS_N)
775 : void readOtherWalkersHills();
776 : /// recompute rct, maxbias and transbias if the bias has changed in this step
777 : void updateBiasQuantities(bool nowAddAHill);
778 : void updateFrequencyAdaptiveStride();
779 : void updateNlist();
780 : /// the position of the neighbor list: the CVs (with INTERVAL, within the limits, where the hills are evaluated)
781 : std::vector<double> getNlistPosition() const;
782 :
783 : public:
784 : explicit MetaD(const ActionOptions&);
785 : void calculate() override;
786 : void update() override;
787 : void beforeUpdate() override;
788 : static void registerKeywords(Keywords& keys);
789 : bool checkNeedsGradients()const override;
790 : };
791 :
792 : PLUMED_REGISTER_ACTION(MetaD,"METAD")
793 :
794 260 : void MetaD::registerKeywords(Keywords& keys) {
795 260 : Bias::registerKeywords(keys);
796 520 : keys.addOutputComponent("rbias","CALC_RCT","scalar","the instantaneous value of the bias normalized using the c(t) reweighting factor [rbias=bias-rct]."
797 : "This component can be used to obtain a reweighted histogram.");
798 520 : keys.addOutputComponent("rct","CALC_RCT","scalar","the reweighting factor c(t).");
799 520 : keys.addOutputComponent("work","CALC_WORK","scalar","accumulator for work");
800 520 : keys.addOutputComponent("acc","ACCELERATION","scalar","the metadynamics acceleration factor");
801 520 : keys.addOutputComponent("maxbias", "CALC_MAX_BIAS", "scalar","the maximum of the metadynamics V(s, t)");
802 520 : keys.addOutputComponent("transbias", "CALC_TRANSITION_BIAS", "scalar","the metadynamics transition bias V*(t)");
803 520 : keys.addOutputComponent("pace","FREQUENCY_ADAPTIVE","scalar","the hill addition frequency when employing frequency adaptive metadynamics");
804 260 : keys.addOutputComponent("nlker","NLIST","scalar","number of hills in the neighbor list");
805 520 : keys.addOutputComponent("nlsteps","NLIST","scalar","number of steps from last neighbor list update");
806 260 : keys.add("compulsory","SIGMA","the widths of the Gaussian hills");
807 260 : keys.add("compulsory","PACE","the frequency for hill addition");
808 260 : keys.add("compulsory","FILE","HILLS","a file in which the list of added hills is stored");
809 260 : keys.add("optional","HEIGHT","the heights of the Gaussian hills. Compulsory unless TAU and either BIASFACTOR or DAMPFACTOR are given");
810 260 : keys.add("optional","FMT","specify format for HILLS files (useful for decrease the number of digits in regtests)");
811 260 : keys.add("optional","BIASFACTOR","use well tempered metadynamics and use this bias factor. Please note you must also specify TEMP. With BIASFACTOR=1 no bias is applied, and the hills files and grids contain the estimate of -F");
812 260 : keys.addFlag("CALC_WORK",false,"calculate the total accumulated work done by the bias since last restart");
813 260 : keys.add("optional","RECT","list of bias factors for all the replicas");
814 260 : keys.add("optional","DAMPFACTOR","damp hills with exp(-max(V)/(kT*DAMPFACTOR)");
815 520 : for (size_t i = 0; i < n_tempering_options_; i++) {
816 260 : registerTemperingKeywords(tempering_names_[i][0], tempering_names_[i][1], keys);
817 : }
818 260 : keys.add("optional","TARGET","target to a predefined distribution");
819 260 : keys.add("optional","TEMP","the system temperature - this is needed (unless the MD engine passes it to PLUMED) with well-tempered metadynamics, DAMPFACTOR, TARGET, ACCELERATION, CALC_RCT and transition-tempered metadynamics");
820 260 : keys.add("optional","TAU","in well tempered metadynamics, sets height to (k_B Delta T*pace*timestep)/tau");
821 260 : keys.addFlag("CALC_RCT",false,"calculate the c(t) reweighting factor and use that to obtain the normalized bias [rbias=bias-rct]."
822 : "This method is not compatible with metadynamics not on a grid.");
823 260 : keys.add("optional","RCT_USTRIDE","the update stride for calculating the c(t) reweighting factor, in units of PACE steps. "
824 : "The default is 1, so c(t) is updated every time the bias is updated.");
825 260 : keys.add("optional","GRID_MIN","the lower bounds for the grid (initial range: the grid is extended when needed)");
826 260 : keys.add("optional","GRID_MAX","the upper bounds for the grid (initial range: the grid is extended when needed)");
827 260 : keys.add("optional","GRID_BIN","the number of bins for the grid");
828 260 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN). "
829 : "Without GRID_MIN and GRID_MAX it gives the grid spacing, and the grid is placed around the first Gaussian");
830 260 : keys.add("compulsory","GRID_SIZE_LIMIT","100000000","the maximum number of points of the grid, which is extended when a Gaussian does not fit in it");
831 260 : keys.addFlag("GRID_SPARSE",false,"use a sparse grid to store hills");
832 260 : keys.addFlag("GRID_NOSPLINE",false,"don't use spline interpolation with grids");
833 260 : keys.add("optional","GRID_WSTRIDE","write the grid to a file every N steps");
834 260 : keys.add("optional","GRID_WFILE","the file on which to write the grid");
835 260 : keys.add("optional","GRID_RFILE","a grid file from which the bias should be read at the initial step of the simulation (PBMETAD: GRID_RFILES)");
836 260 : keys.addFlag("STORE_GRIDS",false,"append every grid that is written to GRID_WFILE, instead of overwriting the previous one");
837 260 : keys.addFlag("NLIST",false,"Use neighbor list for kernels summation, faster but experimental");
838 260 : keys.add("optional", "NLIST_PARAMETERS","(default=6.,0.5) the two cutoff parameters for the Gaussians neighbor list");
839 260 : 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");
840 260 : keys.addFlag("DIFF_FROM_SCRATCH",false,"with ADAPTIVE=DIFF, when restarting accumulate the statistics of the CV fluctuations from scratch, instead of initializing them from the last hills of the hills file");
841 260 : keys.addFlag("PRECAL_DIFF",false,"with ADAPTIVE=DIFF, accumulate the statistics of the CV fluctuations at every step from the beginning of the simulation, also before UPDATE_FROM, after UPDATE_UNTIL and when UPDATE_IF is false, so that the first hills have a width estimated over a sufficiently long time");
842 260 : keys.add("optional","SIGMA_MAX","the upper bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
843 260 : keys.add("optional","SIGMA_MIN","the lower bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
844 260 : MetaDWalkers::registerKeywords(keys);
845 260 : keys.add("optional","INTERVAL","one dimensional lower and upper limits, outside the limits the system will not feel the biasing force (PBMETAD: INTERVAL_MIN and INTERVAL_MAX).");
846 260 : keys.addFlag("FLYING_GAUSSIAN",false,"Switch on flying Gaussian method, must be used with WALKERS_MPI");
847 260 : keys.addFlag("ACCELERATION",false,"Set to TRUE if you want to compute the metadynamics acceleration factor.");
848 260 : keys.add("optional","ACCELERATION_RFILE","a data file from which the acceleration should be read at the initial step of the simulation");
849 260 : keys.addFlag("CALC_MAX_BIAS", false, "Set to TRUE if you want to compute the maximum of the metadynamics V(s, t)");
850 260 : keys.addFlag("CALC_TRANSITION_BIAS", false, "Set to TRUE if you want to compute a metadynamics transition bias V*(t)");
851 260 : 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.");
852 260 : 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.");
853 260 : 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");
854 260 : keys.add("optional","FA_MAX_PACE","the maximum hill addition frequency allowed in frequency adaptive metadynamics. By default there is no maximum value.");
855 260 : 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.");
856 260 : keys.use("RESTART");
857 260 : keys.use("UPDATE_FROM");
858 260 : keys.use("UPDATE_UNTIL");
859 260 : keys.addDOI("10.1073/pnas.202427399");
860 260 : keys.addDOI("10.1063/1.2844595");
861 260 : keys.addDOI("10.1103/PhysRevLett.100.020603");
862 260 : keys.addDOI("10.1021/ct3002464");
863 260 : keys.addDOI("10.2174/1877946811202010079");
864 260 : keys.addDOI("10.1021/jp054359r");
865 260 : keys.addDOI("10.1021/jp504920s");
866 260 : keys.addDOI("10.1021/ct5009087");
867 260 : keys.addDOI("10.1103/PhysRevLett.111.230602");
868 260 : keys.addDOI("10.1063/1.5024679");
869 260 : keys.addDOI("10.1021/acs.jctc.5b00178");
870 260 : keys.addDOI("10.1016/j.bpj.2015.05.024");
871 260 : keys.addDOI("10.1021/acs.jctc.6b00299");
872 260 : keys.addDOI("10.1103/PhysRevLett.112.240602");
873 260 : }
874 :
875 : const std::string MetaD::tempering_names_[1][2] = {{"TT", "transition tempered"}};
876 :
877 260 : void MetaD::registerTemperingKeywords(const std::string &name_stem, const std::string &name, Keywords &keys) {
878 520 : keys.add("optional", name_stem + "BIASFACTOR", "use " + name + " metadynamics with this bias factor. Please note you must also specify temp");
879 780 : keys.add("optional", name_stem + "BIASTHRESHOLD", "use " + name + " metadynamics with this bias threshold. Please note you must also specify " + name_stem + "BIASFACTOR");
880 780 : 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");
881 260 : }
882 :
883 256 : MetaD::MetaD(const ActionOptions& ao):
884 256 : PLUMED_BIAS_INIT(ao) {
885 255 : kernel_.setArguments(getArguments());
886 255 : if(!dp2cutoffNoStretch()) {
887 : kernel_.setStretch(dp2cutoffA,dp2cutoffB);
888 : }
889 255 : SetupInput in;
890 255 : parseGaussianWidth();
891 249 : parseHillsAndTempering(in);
892 247 : parseGrid(in);
893 245 : parseNeighborList();
894 245 : parseReweightingFactor();
895 245 : if(dampfactor_>0.0) {
896 2 : if(!grid_) {
897 0 : error("With DAMPFACTOR you should use grids");
898 : }
899 : }
900 245 : parseMultipleWalkers();
901 242 : parseInterval();
902 242 : parseAcceleration(in);
903 :
904 242 : checkRead();
905 :
906 242 : logHillsSetup(in);
907 242 : checkTempering();
908 242 : logIntervalAndGrid(in);
909 242 : setupMultipleWalkers();
910 242 : setupComponents(in);
911 242 : setupBiasGrid(in);
912 241 : readHillsFiles(in);
913 241 : setupTargetGrid();
914 241 : initializeRestartQuantities(in);
915 241 : openOutputFiles(in);
916 240 : logBibliography(in);
917 495 : }
918 :
919 255 : void MetaD::parseGaussianWidth() {
920 : // parse the flexible hills
921 : std::string adaptiveoption;
922 : adaptiveoption="NONE";
923 255 : parse("ADAPTIVE",adaptiveoption);
924 255 : if(adaptiveoption=="GEOM") {
925 27 : log.printf(" Uses Geometry-based hills width: sigma must be in distance units and only one sigma is needed\n");
926 27 : adaptive_=FlexibleBin::geometry;
927 228 : } else if(adaptiveoption=="DIFF") {
928 25 : log.printf(" Uses Diffusion-based hills width: sigma must be in time steps and only one sigma is needed\n");
929 25 : adaptive_=FlexibleBin::diffusion;
930 203 : } else if(adaptiveoption=="NONE") {
931 202 : adaptive_=FlexibleBin::none;
932 : } else {
933 1 : error("I do not know this type of adaptive scheme");
934 : }
935 :
936 254 : parseFlag("PRECAL_DIFF",precal_diff_);
937 254 : if(precal_diff_) {
938 1 : if(adaptive_!=FlexibleBin::diffusion) {
939 0 : error("PRECAL_DIFF can only be used with ADAPTIVE=DIFF");
940 : }
941 1 : log.printf(" The statistics for the hills width are accumulated at every step, also outside the update window (PRECAL_DIFF)\n");
942 : }
943 254 : bool diff_from_scratch=false;
944 254 : parseFlag("DIFF_FROM_SCRATCH",diff_from_scratch);
945 254 : if(diff_from_scratch && adaptive_!=FlexibleBin::diffusion) {
946 0 : error("DIFF_FROM_SCRATCH can only be used with ADAPTIVE=DIFF");
947 : }
948 : // the hills are kept only if needed (restart with ADAPTIVE=DIFF), see MetaDHillsFiles::readHills
949 254 : diffrestart_.pending=(adaptive_==FlexibleBin::diffusion && getRestart() && !diff_from_scratch);
950 254 : if(adaptive_==FlexibleBin::diffusion && getRestart()) {
951 15 : if(diff_from_scratch) {
952 3 : log.printf(" The statistics for the hills width are accumulated from scratch after the restart (DIFF_FROM_SCRATCH)\n");
953 : } else {
954 12 : log.printf(" The statistics for the hills width are initialized from the last hills of the hills file, if available\n");
955 : }
956 : }
957 :
958 254 : parse("FMT",fmt_);
959 :
960 : // parse the sigma
961 508 : parseVector("SIGMA",sigma0_);
962 656 : for(const auto s : sigma0_) {
963 403 : if(!(s>0.0)) {
964 1 : error("SIGMA should be positive");
965 : }
966 : }
967 253 : if(adaptive_==FlexibleBin::none) {
968 : // if you use normal sigma you need one sigma per argument
969 201 : if( sigma0_.size()!=getNumberOfArguments() ) {
970 0 : error("number of arguments does not match number of SIGMA parameters");
971 : }
972 201 : std::vector<double> smin,smax;
973 201 : parseVector("SIGMA_MIN",smin);
974 402 : parseVector("SIGMA_MAX",smax);
975 201 : if(smin.size()>0 || smax.size()>0) {
976 2 : error("SIGMA_MIN and SIGMA_MAX can only be used with ADAPTIVE");
977 : }
978 : } else {
979 : // if you use flexible hills you need one sigma
980 52 : if(sigma0_.size()!=1) {
981 1 : error("If you choose ADAPTIVE you need only one sigma according to your choice of type (GEOM/DIFF)");
982 : }
983 : // if adaptive then the number must be an integer
984 51 : if(adaptive_==FlexibleBin::diffusion) {
985 25 : if(int(sigma0_[0])-sigma0_[0]>1.e-9 || int(sigma0_[0])-sigma0_[0] <-1.e-9 || int(sigma0_[0])<1 ) {
986 0 : error("In case of adaptive hills with diffusion, the sigma must be an integer which is the number of time steps\n");
987 : }
988 : }
989 : // the limits of the widths
990 102 : parseVector("SIGMA_MIN",sigma0min_);
991 51 : if(sigma0min_.size()>0 && sigma0min_.size()!=getNumberOfArguments()) {
992 1 : error("the number of SIGMA_MIN values should be the same as the number of arguments");
993 : } else if(sigma0min_.size()==0) {
994 41 : sigma0min_.resize(getNumberOfArguments());
995 117 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
996 76 : sigma0min_[i]=-1.;
997 : }
998 : }
999 :
1000 100 : parseVector("SIGMA_MAX",sigma0max_);
1001 50 : if(sigma0max_.size()>0 && sigma0max_.size()!=getNumberOfArguments()) {
1002 1 : error("the number of SIGMA_MAX values should be the same as the number of arguments");
1003 49 : } else if(sigma0max_.size()==0) {
1004 47 : sigma0max_.resize(getNumberOfArguments());
1005 135 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
1006 88 : sigma0max_[i]=-1.;
1007 : }
1008 : }
1009 :
1010 104 : flexbin_=Tools::make_unique<FlexibleBin>(adaptive_,this,sigma0_[0],sigma0min_,sigma0max_);
1011 : }
1012 249 : }
1013 :
1014 249 : void MetaD::parseHillsAndTempering(SetupInput& in) {
1015 : // note: HEIGHT is not compulsory, since one could use the TAU keyword, see below
1016 249 : parse("HEIGHT",height0_);
1017 249 : if(height0_<0.0) {
1018 0 : error("HEIGHT should not be negative");
1019 : }
1020 249 : parse("PACE",stride_);
1021 248 : if(stride_<=0) {
1022 0 : error("frequency for hill addition is nonsensical");
1023 : }
1024 248 : current_stride_ = stride_;
1025 248 : in.hillsfname="HILLS";
1026 248 : parse("FILE",in.hillsfname);
1027 :
1028 : // Manually set to calculate special bias quantities
1029 : // throughout the course of simulation. (These are chosen due to
1030 : // relevance for tempering and event-driven logic as well.)
1031 248 : parseFlag("CALC_MAX_BIAS", calc_max_bias_);
1032 248 : parseFlag("CALC_TRANSITION_BIAS", calc_transition_bias_);
1033 :
1034 496 : parseVector("RECT",in.rect_biasf);
1035 248 : if(in.rect_biasf.size()>0) {
1036 18 : int r=getReplicaIndex();
1037 18 : if(r>=static_cast<int>(in.rect_biasf.size())) {
1038 0 : error("RECT needs one bias factor for each replica");
1039 : }
1040 18 : biasf_=in.rect_biasf[r];
1041 18 : log<<" You are using RECT\n";
1042 : } else {
1043 460 : parse("BIASFACTOR",biasf_);
1044 : }
1045 248 : if( biasf_<1.0 && biasf_!=-1.0) {
1046 0 : error("well tempered bias factor is nonsensical");
1047 : }
1048 248 : parse("DAMPFACTOR",dampfactor_);
1049 248 : kbt_=getkBT();
1050 248 : if(biasf_>=1.0) {
1051 92 : if(kbt_==0.0) {
1052 0 : error("Unless the MD engine passes the temperature to plumed, with well-tempered metad you must specify it using TEMP");
1053 : }
1054 92 : welltemp_=true;
1055 : }
1056 248 : if(dampfactor_>0.0) {
1057 2 : if(kbt_==0.0) {
1058 0 : error("Unless the MD engine passes the temperature to plumed, with damped metad you must specify it using TEMP");
1059 : }
1060 : }
1061 :
1062 248 : parseFlag("CALC_WORK",calc_work_);
1063 :
1064 : // Set transition tempering parameters.
1065 : // Transition wells are read later via calc_transition_bias_.
1066 248 : readTemperingSpecs(tt_specs_);
1067 248 : if (tt_specs_.is_active) {
1068 3 : calc_transition_bias_ = true;
1069 : }
1070 :
1071 : // If any previous option specified to calculate a transition bias,
1072 : // now read the transition wells for that quantity.
1073 248 : if (calc_transition_bias_) {
1074 15 : std::vector<double> tempcoords(getNumberOfArguments());
1075 30 : for (unsigned i = 0; ; i++) {
1076 90 : if (!parseNumberedVector("TRANSITIONWELL", i, tempcoords) ) {
1077 : break;
1078 : }
1079 30 : if (tempcoords.size() != getNumberOfArguments()) {
1080 0 : error("incorrect number of coordinates for transition tempering well");
1081 : }
1082 30 : transitionwells_.push_back(tempcoords);
1083 : }
1084 : }
1085 :
1086 496 : parse("TARGET",targetfilename_);
1087 248 : if(targetfilename_.length()>0 && kbt_==0.0) {
1088 0 : error("with TARGET temperature must be specified");
1089 : }
1090 248 : in.tau=0.0;
1091 248 : parse("TAU",in.tau);
1092 248 : if(in.tau<0.0) {
1093 0 : error("TAU should be positive");
1094 : }
1095 248 : if(in.tau==0.0) {
1096 226 : if(height0_==std::numeric_limits<double>::max()) {
1097 0 : error("At least one between HEIGHT and TAU should be specified");
1098 : }
1099 : // if tau is not set, we compute it here from the other input parameters
1100 226 : if(welltemp_) {
1101 73 : in.tau=(kbt_*(biasf_-1.0))/height0_*getTimeStep()*stride_;
1102 153 : } else if(dampfactor_>0.0) {
1103 0 : in.tau=(kbt_*dampfactor_)/height0_*getTimeStep()*stride_;
1104 : }
1105 : } else {
1106 22 : if(height0_!=std::numeric_limits<double>::max()) {
1107 0 : error("At most one between HEIGHT and TAU should be specified");
1108 : }
1109 22 : if(welltemp_) {
1110 19 : if(biasf_!=1.0) {
1111 15 : height0_=(kbt_*(biasf_-1.0))/in.tau*getTimeStep()*stride_;
1112 : } else {
1113 4 : height0_=kbt_/in.tau*getTimeStep()*stride_; // special case for gamma=1
1114 : }
1115 3 : } else if(dampfactor_>0.0) {
1116 2 : height0_=(kbt_*dampfactor_)/in.tau*getTimeStep()*stride_;
1117 : } else {
1118 1 : error("TAU only makes sense in well-tempered or damped metadynamics");
1119 : }
1120 : }
1121 247 : }
1122 :
1123 247 : void MetaD::parseGrid(SetupInput& in) {
1124 : // Grid Stuff
1125 : // (parseVector checks the number of values: the vectors are empty if the keywords are not given)
1126 247 : in.gmin.assign(getNumberOfArguments(),std::string());
1127 247 : parseVector("GRID_MIN",in.gmin);
1128 247 : in.gmax.assign(getNumberOfArguments(),std::string());
1129 247 : parseVector("GRID_MAX",in.gmax);
1130 247 : in.gbin.assign(getNumberOfArguments(),0);
1131 247 : parseVector("GRID_BIN",in.gbin);
1132 247 : std::vector<double> gspacing(getNumberOfArguments());
1133 247 : parseVector("GRID_SPACING",gspacing);
1134 247 : in.gridnobounds=MetaDBiasGrid::setupBins(*this,getArguments(),in.gmin,in.gmax,in.gbin,gspacing,adaptive_!=FlexibleBin::none,sigma0_,sigma0min_);
1135 246 : if(in.gridnobounds) {
1136 10 : in.gspacing=gspacing;
1137 : }
1138 246 : if(in.gbin.size()>0) {
1139 121 : grid_=true;
1140 : }
1141 246 : parse("GRID_SIZE_LIMIT",gridsizelimit_);
1142 246 : MetaDBiasGrid::checkSizeLimit(*this,gridsizelimit_);
1143 :
1144 246 : in.sparsegrid=false;
1145 246 : parseFlag("GRID_SPARSE",in.sparsegrid);
1146 246 : bool nospline=false;
1147 246 : parseFlag("GRID_NOSPLINE",nospline);
1148 246 : in.spline=!nospline;
1149 246 : parse("GRID_WSTRIDE",wgridstride_);
1150 246 : if(wgridstride_<0) {
1151 0 : error("GRID_WSTRIDE should not be negative");
1152 : }
1153 246 : parse("GRID_WFILE",in.gridfilename);
1154 246 : parseFlag("STORE_GRIDS",storeOldGrids_);
1155 246 : if(grid_ && in.gridfilename.length()>0) {
1156 22 : if(wgridstride_==0 ) {
1157 0 : error("frequency with which to output grid not specified use GRID_WSTRIDE");
1158 : }
1159 : }
1160 246 : if(grid_ && wgridstride_>0) {
1161 22 : if(in.gridfilename.length()==0) {
1162 0 : error("grid filename not specified use GRID_WFILE");
1163 : }
1164 : }
1165 :
1166 246 : parse("GRID_RFILE",in.gridreadfilename);
1167 :
1168 246 : if(!grid_&&(in.gridfilename.length()> 0 || wgridstride_>0)) {
1169 1 : error("To write a grid you need first to define it!");
1170 : }
1171 245 : if(!grid_&&in.gridreadfilename.length()>0) {
1172 2 : error("To read a grid you need first to define it!");
1173 : }
1174 245 : }
1175 :
1176 245 : void MetaD::parseNeighborList() {
1177 : /*setup neighbor list stuff*/
1178 245 : parseFlag("NLIST", nlist_.enabled);
1179 245 : nlist_.center.resize(getNumberOfArguments());
1180 245 : nlist_.dev2.resize(getNumberOfArguments());
1181 245 : if(nlist_.enabled&&grid_) {
1182 0 : error("NLIST and GRID cannot be combined!");
1183 : }
1184 245 : std::vector<double> nlist_param;
1185 490 : parseVector("NLIST_PARAMETERS",nlist_param);
1186 245 : if(nlist_param.size()==0) {
1187 245 : nlist_.param[0]=6.0;//*DP2CUTOFF -> max distance of neighbors
1188 245 : nlist_.param[1]=0.5;//*nlist_dev2_[i] -> condition for rebuilding
1189 : } else {
1190 0 : if(nlist_param.size()!=2) {
1191 0 : error("two cutoff parameters are needed for the neighbor list");
1192 : }
1193 0 : if(!(nlist_param[0]>1.0)) {
1194 0 : error("the first of NLIST_PARAMETERS must be greater than 1. The smaller the first, the smaller should be the second as well");
1195 : }
1196 0 : const double min_PARAM_1=(1.-1./std::sqrt(nlist_param[0]/2))+0.16;
1197 0 : if(!(nlist_param[1]>0)) {
1198 0 : error("the second of NLIST_PARAMETERS must be greater than 0");
1199 : }
1200 0 : if(nlist_param[1]>min_PARAM_1) {
1201 0 : error("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));
1202 : }
1203 0 : nlist_.param[0]=nlist_param[0];
1204 0 : nlist_.param[1]=nlist_param[1];
1205 : }
1206 245 : }
1207 :
1208 245 : void MetaD::parseReweightingFactor() {
1209 : // Reweighting factor rct
1210 245 : parseFlag("CALC_RCT",calc_rct_);
1211 245 : if (calc_rct_) {
1212 15 : if(!grid_) {
1213 0 : error("CALC_RCT is supported only if bias is on a grid");
1214 : }
1215 15 : if(kbt_==0.0) {
1216 0 : error("Unless the MD engine passes the temperature to plumed, with CALC_RCT you must specify it using TEMP");
1217 : }
1218 : }
1219 245 : parse("RCT_USTRIDE",rct_ustride_);
1220 245 : if(rct_ustride_==0) {
1221 0 : error("RCT_USTRIDE should be positive");
1222 : }
1223 245 : }
1224 :
1225 245 : void MetaD::parseMultipleWalkers() {
1226 : // Multiple walkers
1227 245 : mw_.parse(*this);
1228 :
1229 : // Flying Gaussian
1230 242 : parseFlag("FLYING_GAUSSIAN", flying_);
1231 242 : if(flying_ && grid_) {
1232 : // the hills of the flying Gaussians are replaced at every deposition, which cannot be done on a grid
1233 0 : error("FLYING_GAUSSIAN cannot be used with a grid");
1234 : }
1235 242 : }
1236 :
1237 242 : void MetaD::parseInterval() {
1238 : // Interval keyword
1239 : // (parseVector checks that there are two values: the vector is empty if the keyword is not given)
1240 242 : std::vector<double> tmpI(2);
1241 484 : parseVector("INTERVAL",tmpI);
1242 242 : if(tmpI.size()==2) {
1243 9 : interval_.lower=tmpI.at(0);
1244 9 : interval_.upper=tmpI.at(1);
1245 9 : if(getNumberOfArguments()!=1) {
1246 0 : error("INTERVAL limits correction works only for monodimensional metadynamics!");
1247 : }
1248 9 : if(interval_.upper<interval_.lower) {
1249 0 : error("The Upper limit must be greater than the Lower limit!");
1250 : }
1251 9 : if(getPntrToArgument(0)->isPeriodic()) {
1252 0 : error("INTERVAL cannot be used with periodic variables!");
1253 : }
1254 9 : interval_.enabled=true;
1255 9 : kernel_.setInterval(interval_.lower,interval_.upper);
1256 : }
1257 242 : }
1258 :
1259 242 : void MetaD::parseAcceleration(SetupInput& in) {
1260 242 : parseFlag("ACCELERATION",acceleration_);
1261 : // Check for a restart acceleration if acceleration is active.
1262 242 : if(acceleration_) {
1263 10 : parse("ACCELERATION_RFILE", in.acc_rfilename);
1264 : }
1265 :
1266 242 : fa_.enabled=false;
1267 242 : parseFlag("FREQUENCY_ADAPTIVE",fa_.enabled);
1268 : //
1269 242 : fa_.update_frequency=0;
1270 242 : parse("FA_UPDATE_FREQUENCY",fa_.update_frequency);
1271 242 : if(fa_.update_frequency!=0 && !fa_.enabled) {
1272 0 : error("It doesn't make sense to use the FA_UPDATE_FREQUENCY keyword if frequency adaptive METAD hasn't been activated by using the FREQUENCY_ADAPTIVE flag");
1273 : }
1274 242 : if(fa_.update_frequency==0 && fa_.enabled) {
1275 0 : fa_.update_frequency=stride_;
1276 : }
1277 : //
1278 242 : fa_.max_stride=0;
1279 242 : parse("FA_MAX_PACE",fa_.max_stride);
1280 242 : if(fa_.max_stride!=0 && !fa_.enabled) {
1281 0 : error("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");
1282 : }
1283 : //
1284 242 : fa_.min_acceleration=1.0;
1285 242 : parse("FA_MIN_ACCELERATION",fa_.min_acceleration);
1286 242 : if(fa_.min_acceleration!=1.0 && !fa_.enabled) {
1287 0 : error("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");
1288 : }
1289 242 : }
1290 :
1291 242 : void MetaD::logHillsSetup(const SetupInput& in) {
1292 242 : log.printf(" Gaussian width ");
1293 242 : if (adaptive_==FlexibleBin::diffusion) {
1294 24 : log.printf(" (Note: The units of sigma are in timesteps) ");
1295 : }
1296 242 : if (adaptive_==FlexibleBin::geometry) {
1297 24 : log.printf(" (Note: The units of sigma are in dist units) ");
1298 : }
1299 631 : for(unsigned i=0; i<sigma0_.size(); ++i) {
1300 389 : log.printf(" %f",sigma0_[i]);
1301 : }
1302 242 : log.printf("\n Gaussian height %f\n",height0_);
1303 242 : log.printf(" Gaussian deposition pace %d\n",stride_);
1304 242 : log.printf(" Gaussian file %s\n",in.hillsfname.c_str());
1305 242 : if(biasf_==1.0) {
1306 4 : log.printf(" Bias factor 1: no bias is applied, the hills files and grids contain the estimate of -F\n");
1307 : }
1308 242 : if(welltemp_) {
1309 92 : log.printf(" Well-Tempered Bias Factor %f\n",biasf_);
1310 92 : log.printf(" Hills relaxation time (tau) %f\n",in.tau);
1311 92 : log.printf(" KbT %f\n",kbt_);
1312 : }
1313 242 : }
1314 :
1315 242 : void MetaD::checkTempering() {
1316 : // Transition tempered metadynamics options
1317 242 : if (tt_specs_.is_active) {
1318 3 : logTemperingSpecs(tt_specs_);
1319 : }
1320 :
1321 : // Overall tempering sanity check (this gets tricky when multiple are active).
1322 : // When multiple temperings are active, it's fine to have one tempering attempt
1323 : // to increase hill size with increasing bias, so long as the others can shrink
1324 : // the hills faster than it increases their size in the long-time limit.
1325 : // This set of checks ensures that the hill sizes eventually decay to zero as c(t)
1326 : // diverges to infinity.
1327 : // The alpha parameter allows hills to decay as 1/t^alpha instead of 1/t,
1328 : // a slower decay, so as t -> infinity, only the temperings with the largest
1329 : // alphas govern the final asymptotic decay. (Alpha helps prevent false convergence.)
1330 242 : if (welltemp_ || dampfactor_ > 0.0 || tt_specs_.is_active) {
1331 : // Determine the number of active temperings.
1332 : int n_active = 0;
1333 97 : if (welltemp_) {
1334 : n_active++;
1335 : }
1336 97 : if (dampfactor_ > 0.0) {
1337 2 : n_active++;
1338 : }
1339 97 : if (tt_specs_.is_active) {
1340 3 : n_active++;
1341 : }
1342 : // Find the greatest alpha.
1343 97 : double greatest_alpha = 0.0;
1344 97 : if (welltemp_) {
1345 92 : greatest_alpha = std::max(greatest_alpha, 1.0);
1346 : }
1347 97 : if (dampfactor_ > 0.0) {
1348 4 : greatest_alpha = std::max(greatest_alpha, 1.0);
1349 : }
1350 97 : if (tt_specs_.is_active) {
1351 6 : greatest_alpha = std::max(greatest_alpha, tt_specs_.alpha);
1352 : }
1353 : // Find the least alpha.
1354 97 : double least_alpha = 1.0;
1355 97 : if (welltemp_) {
1356 : least_alpha = std::min(least_alpha, 1.0);
1357 : }
1358 97 : if (dampfactor_ > 0.0) {
1359 2 : least_alpha = std::min(least_alpha, 1.0);
1360 : }
1361 97 : if (tt_specs_.is_active) {
1362 4 : least_alpha = std::min(least_alpha, tt_specs_.alpha);
1363 : }
1364 : // Find the inverse harmonic average of the delta T parameters for all
1365 : // of the temperings with the greatest alpha values.
1366 : double total_governing_deltaT_inv = 0.0;
1367 97 : if (welltemp_ && 1.0 == greatest_alpha && biasf_ != 1.0) {
1368 88 : total_governing_deltaT_inv += 1.0 / (biasf_ - 1.0);
1369 : }
1370 97 : if (dampfactor_ > 0.0 && 1.0 == greatest_alpha) {
1371 2 : total_governing_deltaT_inv += 1.0 / (dampfactor_);
1372 : }
1373 97 : if (tt_specs_.is_active && tt_specs_.alpha == greatest_alpha) {
1374 3 : total_governing_deltaT_inv += 1.0 / (tt_specs_.biasf - 1.0);
1375 : }
1376 : // Give a newbie-friendly error message for people using one tempering if
1377 : // only one is active.
1378 97 : if (n_active == 1 && total_governing_deltaT_inv < 0.0) {
1379 0 : error("for stable tempering, the bias factor must be greater than one");
1380 : // Give a slightly more complex error message to users stacking multiple
1381 : // tempering options at a time, but all with uniform alpha values.
1382 0 : } else if (total_governing_deltaT_inv < 0.0 && greatest_alpha == least_alpha) {
1383 0 : error("for stable tempering, the sum of the inverse Delta T parameters must be greater than zero!");
1384 : // Give the most technical error message to users stacking multiple tempering
1385 : // options with different alpha parameters.
1386 97 : } else if (total_governing_deltaT_inv < 0.0 && greatest_alpha != least_alpha) {
1387 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!");
1388 : }
1389 : }
1390 242 : }
1391 :
1392 242 : void MetaD::logIntervalAndGrid(const SetupInput& in) {
1393 242 : if(interval_.enabled) {
1394 9 : log.printf(" Upper and Lower limits boundaries for the bias are activated at %f - %f\n", interval_.lower, interval_.upper);
1395 : }
1396 :
1397 242 : if(grid_ && in.gridnobounds) {
1398 10 : log.printf(" Grid spacing");
1399 25 : for(unsigned i=0; i<in.gspacing.size(); ++i) {
1400 30 : log.printf(" %s",exactString(in.gspacing[i]).c_str());
1401 : }
1402 10 : log.printf(" (periodic variables: grid bin");
1403 25 : for(unsigned i=0; i<in.gbin.size(); ++i) {
1404 15 : if(getPntrToArgument(i)->isPeriodic()) {
1405 5 : log.printf(" %u",in.gbin[i]);
1406 : } else {
1407 10 : log.printf(" -");
1408 : }
1409 : }
1410 10 : log.printf(")\n The grid is placed around the first Gaussian\n");
1411 242 : } else if(grid_) {
1412 111 : log.printf(" Grid min");
1413 304 : for(unsigned i=0; i<in.gmin.size(); ++i) {
1414 193 : log.printf(" %s",in.gmin[i].c_str() );
1415 : }
1416 111 : log.printf("\n");
1417 111 : log.printf(" Grid max");
1418 304 : for(unsigned i=0; i<in.gmax.size(); ++i) {
1419 193 : log.printf(" %s",in.gmax[i].c_str() );
1420 : }
1421 111 : log.printf("\n");
1422 111 : log.printf(" Grid bin");
1423 304 : for(unsigned i=0; i<in.gbin.size(); ++i) {
1424 193 : log.printf(" %u",in.gbin[i]);
1425 : }
1426 111 : log.printf("\n");
1427 : }
1428 242 : if(grid_) {
1429 121 : if(in.spline) {
1430 120 : log.printf(" Grid uses spline interpolation\n");
1431 : }
1432 121 : if(in.sparsegrid) {
1433 12 : log.printf(" Grid uses sparse grid\n");
1434 : }
1435 121 : if(wgridstride_>0) {
1436 22 : log.printf(" Grid is written on file %s with stride %d\n",in.gridfilename.c_str(),wgridstride_);
1437 : }
1438 242 : log.printf(" The grid is extended when a Gaussian does not fit in it, up to %s points\n",exactString(gridsizelimit_).c_str());
1439 : }
1440 242 : }
1441 :
1442 242 : void MetaD::setupMultipleWalkers() {
1443 242 : mw_.setup(*this);
1444 :
1445 242 : if(flying_) {
1446 6 : if(!mw_.mpi) {
1447 0 : error("Flying Gaussian method must be used with MPI version of multiple walkers");
1448 : }
1449 6 : log.printf(" Flying Gaussian method with %d walkers active\n",mw_.mpi_nw);
1450 : }
1451 242 : }
1452 :
1453 242 : void MetaD::setupComponents(const SetupInput& in) {
1454 242 : if(nlist_.enabled) {
1455 8 : addComponent("nlker");
1456 8 : componentIsNotPeriodic("nlker");
1457 8 : addComponent("nlsteps");
1458 8 : componentIsNotPeriodic("nlsteps");
1459 : }
1460 :
1461 242 : if(calc_rct_) {
1462 30 : addComponent("rbias");
1463 30 : componentIsNotPeriodic("rbias");
1464 30 : addComponent("rct");
1465 15 : componentIsNotPeriodic("rct");
1466 15 : log.printf(" The c(t) reweighting factor will be calculated every %u hills\n",rct_ustride_);
1467 30 : getPntrToComponent("rct")->set(reweight_factor_);
1468 : }
1469 :
1470 242 : if(calc_work_) {
1471 4 : addComponent("work");
1472 4 : componentIsNotPeriodic("work");
1473 : }
1474 :
1475 242 : if(acceleration_) {
1476 5 : if (kbt_ == 0.0) {
1477 0 : error("The calculation of the acceleration works only if simulation temperature has been defined");
1478 : }
1479 5 : log.printf(" calculation on the fly of the acceleration factor\n");
1480 10 : addComponent("acc");
1481 10 : componentIsNotPeriodic("acc");
1482 : // Set the initial value of the acceleration.
1483 : // If this is not a restart, set to 1.0.
1484 5 : if (in.acc_rfilename.length() == 0) {
1485 3 : getPntrToComponent("acc")->set(1.0);
1486 3 : if(getRestart()) {
1487 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");
1488 1 : log.printf(" You should use the ACCELERATION_RFILE keyword.\n");
1489 : }
1490 : // Otherwise, read and set the restart value.
1491 : } else {
1492 : // Open the ACCELERATION_RFILE.
1493 2 : IFile acc_rfile;
1494 2 : acc_rfile.link(*this);
1495 2 : if(acc_rfile.FileExist(in.acc_rfilename)) {
1496 2 : acc_rfile.open(in.acc_rfilename);
1497 : } else {
1498 0 : error("The ACCELERATION_RFILE file you want to read: " + in.acc_rfilename + ", cannot be found!");
1499 : }
1500 : // Read the file to find the restart acceleration.
1501 2 : double acc_rmean=0.0;
1502 2 : double acc_rtime=0.0;
1503 : bool found=false;
1504 2 : std::string acclabel = getLabel() + ".acc";
1505 : // with FREQUENCY_ADAPTIVE also the pace (it never decreases, so it is not a function of the current acceleration)
1506 2 : const std::string pacelabel = getLabel() + ".pace";
1507 2 : double acc_rpace=0.0;
1508 2 : acc_rfile.allowIgnoredFields();
1509 248 : while(acc_rfile.scanField("time", acc_rtime)) {
1510 122 : acc_rfile.scanField(acclabel, acc_rmean);
1511 122 : if(fa_.enabled && acc_rfile.FieldExist(pacelabel)) {
1512 101 : acc_rfile.scanField(pacelabel, acc_rpace);
1513 : }
1514 122 : acc_rfile.scanField();
1515 : found=true;
1516 : }
1517 2 : if(fa_.enabled && acc_rpace>current_stride_) {
1518 1 : current_stride_=static_cast<int>(std::lround(acc_rpace));
1519 1 : log.printf(" hill addition pace read from file %s: %d\n",in.acc_rfilename.c_str(),current_stride_);
1520 : }
1521 2 : acc_restart_step_=std::llround(acc_rtime/getTimeStep());
1522 2 : if(!found) {
1523 0 : error("The ACCELERATION_RFILE file you want to read: " + in.acc_rfilename + ", does not contain a time field!");
1524 : }
1525 2 : acc_restart_mean_ = acc_rmean;
1526 : // Set component based on the read values.
1527 2 : getPntrToComponent("acc")->set(acc_rmean);
1528 2 : log.printf(" initial acceleration factor read from file %s: value of %f at time %f\n",in.acc_rfilename.c_str(),acc_rmean,acc_rtime);
1529 2 : }
1530 : }
1531 :
1532 242 : if (calc_max_bias_) {
1533 1 : if (!grid_) {
1534 0 : error("Calculating the maximum bias on the fly works only with a grid");
1535 : }
1536 1 : log.printf(" calculation on the fly of the maximum bias max(V(s,t)) \n");
1537 2 : addComponent("maxbias");
1538 2 : componentIsNotPeriodic("maxbias");
1539 : }
1540 :
1541 242 : if (calc_transition_bias_) {
1542 15 : if (!grid_) {
1543 0 : error("Calculating the transition bias on the fly works only with a grid");
1544 : }
1545 15 : log.printf(" calculation on the fly of the transition bias V*(t)\n");
1546 30 : addComponent("transbias");
1547 15 : componentIsNotPeriodic("transbias");
1548 15 : log<<" Number of transition wells "<<transitionwells_.size()<<"\n";
1549 15 : if (transitionwells_.size() == 0) {
1550 0 : error("Calculating the transition bias on the fly requires definition of at least one transition well");
1551 : }
1552 : // Log the wells (the grid is extended to contain them in setupBiasGrid)
1553 45 : for (unsigned i = 0; i < transitionwells_.size(); i++) {
1554 : // Log the coordinate.
1555 30 : log.printf(" Transition well %d at coordinate ", i);
1556 76 : for (unsigned j = 0; j < getNumberOfArguments(); j++) {
1557 46 : log.printf("%f ", transitionwells_[i][j]);
1558 : }
1559 30 : log.printf("\n");
1560 : }
1561 : }
1562 :
1563 242 : if(fa_.enabled) {
1564 3 : if(!acceleration_) {
1565 0 : error("Frequency adaptive metadynamics only works if the calculation of the acceleration factor is enabled with the ACCELERATION keyword");
1566 : }
1567 3 : if(mw_.mpi) {
1568 0 : error("Combining frequency adaptive metadynamics with MPI multiple walkers is not allowed");
1569 : }
1570 :
1571 3 : log.printf(" Frequency adaptive metadynamics enabled\n");
1572 3 : if(getRestart() && in.acc_rfilename.length() == 0) {
1573 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");
1574 0 : log.printf(" You should use the ACCELERATION_RFILE keyword.\n");
1575 : }
1576 3 : log.printf(" The frequency for hill addition will change dynamically based on the metadynamics acceleration factor\n");
1577 3 : log.printf(" The hill addition frequency will be updated every %d steps\n",fa_.update_frequency);
1578 3 : if(fa_.min_acceleration>1.0) {
1579 2 : log.printf(" The hill addition frequency will only be updated once the metadynamics acceleration factor becomes larger than %.1f \n",fa_.min_acceleration);
1580 : }
1581 3 : if(fa_.max_stride!=0) {
1582 3 : log.printf(" The hill addition frequency will not become larger than %d steps\n",fa_.max_stride);
1583 : }
1584 6 : addComponent("pace");
1585 3 : componentIsNotPeriodic("pace");
1586 3 : updateFrequencyAdaptiveStride();
1587 : }
1588 242 : }
1589 :
1590 242 : void MetaD::setupBiasGrid(SetupInput& in) {
1591 : // initializing and checking grid
1592 242 : in.restartedFromGrid=false; // restart from grid file
1593 242 : if(grid_) {
1594 121 : std::string funcl=getLabel() + ".bias";
1595 121 : BiasGrid_=Tools::make_unique<MetaDBiasGrid>(*this);
1596 : // check for mesh and sigma size (for a grid read from file, on the lattice of the grid in input)
1597 121 : MetaDBiasGrid::checkSpacing(*this,in.gmin,in.gmax,in.gbin,adaptive_!=FlexibleBin::none,sigma0_,sigma0min_);
1598 121 : if(!(in.gridreadfilename.length()>0)) {
1599 88 : BiasGrid_->create(funcl,getArguments(),in.gmin,in.gmax,in.gbin,in.sparsegrid,in.spline,in.gridnobounds);
1600 : } else {
1601 : // read the grid in input
1602 65 : in.gridWalkerHills=readBiasGrid(*this,mw_,*BiasGrid_,in.gridreadfilename,funcl,getArguments(),in.gmin,in.gmax,in.gbin,
1603 65 : in.sparsegrid,in.spline,in.gridnobounds);
1604 32 : if(getRestart()) {
1605 32 : in.restartedFromGrid=true;
1606 : }
1607 : }
1608 120 : BiasGrid_->setSizeLimit(gridsizelimit_);
1609 : // with INTERVAL the bias beyond the limits is constant: the grid contains the interval and at least one point
1610 : // beyond each limit, which holds this constant value
1611 120 : if(interval_.enabled) {
1612 6 : BiasGrid_->setInterval(interval_.lower,interval_.upper,getTime());
1613 : }
1614 : // the transition wells must be in the grid (it is extended along the CVs that are not periodic)
1615 120 : if(calc_transition_bias_) {
1616 : const GridBase& grid=BiasGrid_->getGrid();
1617 16 : const std::vector<double> zero(getNumberOfArguments(),0.0);
1618 45 : for(const auto & well : transitionwells_) {
1619 76 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
1620 46 : if(getPntrToArgument(i)->isPeriodic()) {
1621 20 : const double k=grid.getLatticeCoordinate(i,well[i]);
1622 40 : if(!(k>=0.0 && k<static_cast<double>(grid.getNbin()[i]))) {
1623 0 : error(" transition well is not in grid");
1624 : }
1625 : }
1626 : }
1627 30 : BiasGrid_->extend(well,zero,getTime());
1628 : }
1629 : }
1630 : }
1631 241 : }
1632 :
1633 241 : void MetaD::readHillsFiles(const SetupInput& in) {
1634 : // open all files at the beginning and read Gaussians if restarting
1635 482 : hillsfiles_=Tools::make_unique<MetaDHillsFiles>(*this,comm,getArguments());
1636 : MetaDHillsFiles::Bias bias;
1637 241 : bias.add=[this](const Gaussian& hill) {
1638 3393 : addGaussian(hill);
1639 241 : };
1640 241 : bias.grid=BiasGrid_.get();
1641 241 : bias.diffrestart=&diffrestart_;
1642 : // note that for gamma=1 we store directly -F
1643 241 : if(welltemp_ && biasf_>1.0) {
1644 88 : bias.heightscale=(biasf_-1.0)/biasf_;
1645 : }
1646 241 : bias.tolerance=0.5*stride_*getTimeStep();
1647 : hillsfiles_->setBias(bias);
1648 241 : if(diffrestart_.pending) {
1649 : // hills of the last ~6 averaging times (older ones have weight < exp(-6));
1650 : // with MPI walkers the file contains the hills of all walkers
1651 14 : diffrestart_.maxhills=(static_cast<unsigned>(std::ceil(6.0*sigma0_[0]/stride_))+2)*std::max(1u,mw_.mpi_nw);
1652 : }
1653 : MetaDHillsInput input;
1654 241 : input.files=hillsfiles_.get();
1655 241 : input.name=in.hillsfname;
1656 241 : input.gridWalkerHills=in.gridWalkerHills;
1657 241 : readHillsAtStart(*this,mw_,std::vector<MetaDHillsInput>(1,input),in.restartedFromGrid,schedule_);
1658 241 : }
1659 :
1660 241 : void MetaD::setupTargetGrid() {
1661 241 : if(targetfilename_.length()>0) {
1662 2 : IFile gridfile;
1663 2 : gridfile.open(targetfilename_);
1664 2 : std::string funcl=getLabel() + ".target";
1665 4 : TargetGrid_=GridBase::create(funcl,getArguments(),gridfile,false,false,true);
1666 2 : if(TargetGrid_->getDimension()!=getNumberOfArguments()) {
1667 0 : error("mismatch between dimensionality of input grid and number of arguments");
1668 : }
1669 4 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1670 4 : if( getPntrToArgument(i)->isPeriodic()!=TargetGrid_->getIsPeriodic()[i] ) {
1671 0 : error("periodicity mismatch between arguments and input bias");
1672 : }
1673 : }
1674 2 : }
1675 241 : }
1676 :
1677 241 : void MetaD::initializeRestartQuantities(const SetupInput& in) {
1678 : // (a grid can be read also without RESTART)
1679 241 : if(getRestart() || in.gridreadfilename.length()>0) {
1680 : // if this is a restart the neighbor list should be immediately updated
1681 72 : if(nlist_.enabled) {
1682 1 : nlist_.update=true;
1683 : }
1684 : // Calculate the Tiwary-Parrinello reweighting factor if we are restarting from previous hills
1685 72 : if(calc_rct_) {
1686 0 : computeReweightingFactor();
1687 : }
1688 : // Calculate all special bias quantities desired if restarting with nonzero bias.
1689 72 : if(calc_max_bias_) {
1690 0 : max_bias_ = BiasGrid_->getGrid().getMaxValue();
1691 0 : getPntrToComponent("maxbias")->set(max_bias_);
1692 : }
1693 72 : if(calc_transition_bias_) {
1694 13 : transition_bias_ = getTransitionBarrierBias();
1695 26 : getPntrToComponent("transbias")->set(transition_bias_);
1696 : }
1697 : }
1698 241 : }
1699 :
1700 241 : void MetaD::openOutputFiles(SetupInput& in) {
1701 : // open grid file for writing
1702 241 : if(wgridstride_>0) {
1703 22 : mw_.openOutput(*this,gridfile_,in.gridfilename);
1704 : }
1705 :
1706 : // open hills file for writing
1707 241 : mw_.openOutput(*this,hillsOfile_,hillsfiles_->getInputName(mw_.id));
1708 240 : if(fmt_.length()>0) {
1709 199 : hillsOfile_.fmtField(fmt_);
1710 : }
1711 240 : hillsOfile_.addConstantField("multivariate");
1712 240 : hillsOfile_.addConstantField("kerneltype");
1713 240 : if(interval_.enabled) {
1714 18 : hillsOfile_.addConstantField("lower_int").printField("lower_int",interval_.lower);
1715 18 : hillsOfile_.addConstantField("upper_int").printField("upper_int",interval_.upper);
1716 : }
1717 : hillsOfile_.setHeavyFlush();
1718 : // output periodicities of variables
1719 669 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1720 429 : hillsOfile_.setupPrintValue( getPntrToArgument(i) );
1721 : }
1722 240 : }
1723 :
1724 240 : void MetaD::logBibliography(const SetupInput& in) {
1725 : bool concurrent=false;
1726 240 : const ActionSet&actionSet(plumed.getActionSet());
1727 3156 : for(const auto & p : actionSet)
1728 3051 : if(dynamic_cast<MetaD*>(p.get())) {
1729 : concurrent=true;
1730 : break;
1731 : }
1732 240 : if(concurrent) {
1733 135 : log<<" You are using concurrent metadynamics\n";
1734 : }
1735 240 : if(in.rect_biasf.size()>0) {
1736 18 : if(mw_.mpi) {
1737 12 : log<<" You are using RECT in its 'altruistic' implementation\n";
1738 : } else {
1739 6 : log<<" You are using RECT\n";
1740 : }
1741 : }
1742 :
1743 240 : log<<" Bibliography "<<plumed.cite("Laio and Parrinello, PNAS 99, 12562 (2002)");
1744 240 : if(welltemp_) {
1745 184 : log<<plumed.cite("Barducci, Bussi, and Parrinello, Phys. Rev. Lett. 100, 020603 (2008)");
1746 : }
1747 240 : if(tt_specs_.is_active) {
1748 3 : log << plumed.cite("Dama, Rotskoff, Parrinello, and Voth, J. Chem. Theory Comput. 10, 3626 (2014)");
1749 6 : log << plumed.cite("Dama, Parrinello, and Voth, Phys. Rev. Lett. 112, 240602 (2014)");
1750 : }
1751 240 : if(mw_.n>1||mw_.mpi) {
1752 136 : log<<plumed.cite("Raiteri, Laio, Gervasio, Micheletti, and Parrinello, J. Phys. Chem. B 110, 3533 (2006)");
1753 : }
1754 240 : if(adaptive_!=FlexibleBin::none) {
1755 96 : log<<plumed.cite("Branduardi, Bussi, and Parrinello, J. Chem. Theory Comput. 8, 2247 (2012)");
1756 : }
1757 240 : if(interval_.enabled) {
1758 18 : log<<plumed.cite("Baftizadeh, Cossio, Pietrucci, and Laio, Curr. Phys. Chem. 2, 79 (2012)");
1759 : }
1760 240 : if(acceleration_) {
1761 10 : log<<plumed.cite("Pratyush and Parrinello, Phys. Rev. Lett. 111, 230602 (2013)");
1762 : }
1763 240 : if(calc_rct_) {
1764 30 : log<<plumed.cite("Pratyush and Parrinello, J. Phys. Chem. B, 119, 736 (2015)");
1765 : }
1766 240 : if(concurrent || in.rect_biasf.size()>0) {
1767 282 : log<<plumed.cite("Gil-Ley and Bussi, J. Chem. Theory Comput. 11, 1077 (2015)");
1768 : }
1769 240 : if(in.rect_biasf.size()>0 && mw_.mpi) {
1770 24 : log<<plumed.cite("Hosek, Toulcova, Bortolato, and Spiwok, J. Phys. Chem. B 120, 2209 (2016)");
1771 : }
1772 240 : if(targetfilename_.length()>0) {
1773 2 : log<<plumed.cite("White, Dama, and Voth, J. Chem. Theory Comput. 11, 2451 (2015)");
1774 2 : log<<plumed.cite("Marinelli and Faraldo-Gómez, Biophys. J. 108, 2779 (2015)");
1775 4 : log<<plumed.cite("Gil-Ley, Bottaro, and Bussi, J. Chem. Theory Comput. 12, 2790 (2016)");
1776 : }
1777 240 : if(fa_.enabled) {
1778 6 : log<<plumed.cite("Wang, Valsson, Tiwary, Parrinello, and Lindorff-Larsen, J. Chem. Phys. 149, 072309 (2018)");
1779 : }
1780 240 : log<<"\n";
1781 240 : }
1782 :
1783 248 : void MetaD::readTemperingSpecs(TemperingSpecs &t_specs) {
1784 : // Set global tempering parameters.
1785 248 : parse(t_specs.name_stem + "BIASFACTOR", t_specs.biasf);
1786 248 : if (t_specs.biasf != -1.0) {
1787 3 : if (kbt_ == 0.0) {
1788 0 : error("Unless the MD engine passes the temperature to plumed, with tempered metad you must specify it using TEMP");
1789 : }
1790 3 : if (t_specs.biasf == 1.0) {
1791 0 : error("A bias factor of 1 corresponds to zero delta T and zero hill size, so it is not allowed.");
1792 : }
1793 3 : if (t_specs.biasf < 1.0) {
1794 0 : error(t_specs.name + " bias factor is nonsensical");
1795 : }
1796 3 : t_specs.is_active = true;
1797 3 : parse(t_specs.name_stem + "BIASTHRESHOLD", t_specs.threshold);
1798 3 : if (t_specs.threshold < 0.0) {
1799 0 : error(t_specs.name + " bias threshold is nonsensical");
1800 : }
1801 3 : parse(t_specs.name_stem + "ALPHA", t_specs.alpha);
1802 3 : if (t_specs.alpha <= 0.0 || t_specs.alpha > 1.0) {
1803 0 : error(t_specs.name + " decay shape parameter alpha is nonsensical");
1804 : }
1805 : }
1806 248 : }
1807 :
1808 3 : void MetaD::logTemperingSpecs(const TemperingSpecs &t_specs) {
1809 3 : log.printf(" %s bias factor %f\n", t_specs.name.c_str(), t_specs.biasf);
1810 3 : log.printf(" KbT %f\n", kbt_);
1811 3 : if (t_specs.threshold != 0.0) {
1812 2 : log.printf(" %s bias threshold %f\n", t_specs.name.c_str(), t_specs.threshold);
1813 : }
1814 3 : if (t_specs.alpha != 1.0) {
1815 1 : log.printf(" %s decay shape parameter alpha %f\n", t_specs.name.c_str(), t_specs.alpha);
1816 : }
1817 3 : }
1818 :
1819 10293 : void MetaD::writeGaussian(const Gaussian& hill, OFile&file) {
1820 10293 : double height=hill.height;
1821 : // note that for gamma=1 we store directly -F
1822 10293 : if(welltemp_ && biasf_>1.0) {
1823 5273 : height*=biasf_/(biasf_-1.0);
1824 : }
1825 10293 : hillsfiles_->write(file,getTimeStep()*getStep(),hill,height,biasf_,mw_.n>1);
1826 10293 : }
1827 :
1828 13848 : void MetaD::addGaussian(const Gaussian& hill) {
1829 13848 : if(grid_) {
1830 6090 : BiasGrid_->addHill(hill,kernel_,getTime());
1831 : } else {
1832 7758 : hills_.push_back(hill);
1833 7758 : const unsigned ncv=getNumberOfArguments();
1834 22570 : for(unsigned i=0; i<ncv; i++) {
1835 14812 : hillarrays_.center.push_back(hill.center[i]);
1836 14812 : hillarrays_.invsigma.push_back(hill.multivariate ? 0.0 : hill.invsigma[i]);
1837 : }
1838 7758 : hillarrays_.height.push_back(hill.height);
1839 7758 : hillarrays_.diagonal.push_back(!hill.multivariate);
1840 : }
1841 13848 : }
1842 :
1843 28251 : double MetaD::evaluateHill(const std::vector<double>& cv, unsigned h) const {
1844 28251 : if(hillarrays_.diagonal[h] && !kernel_.hasInterval()) {
1845 : const std::size_t ncv=cv.size();
1846 2512 : return kernel_.evaluateDiagonal(cv,&hillarrays_.center[h*ncv],&hillarrays_.invsigma[h*ncv],hillarrays_.height[h]);
1847 : }
1848 25739 : return kernel_.evaluate(cv,hills_[h]);
1849 : }
1850 :
1851 2926381 : double MetaD::evaluateHillWithDerivatives(const std::vector<double>& cv, unsigned h, double* der, double* dp) const {
1852 2926381 : if(hillarrays_.diagonal[h] && !kernel_.hasInterval()) {
1853 : const std::size_t ncv=cv.size();
1854 2470894 : return kernel_.evaluateDiagonalWithDerivatives(cv,&hillarrays_.center[h*ncv],&hillarrays_.invsigma[h*ncv],hillarrays_.height[h],der,dp);
1855 : }
1856 455487 : return kernel_.evaluateWithDerivatives(cv,hills_[h],der,dp);
1857 : }
1858 :
1859 5683 : double MetaD::getBias(const std::vector<double>& cv) {
1860 5683 : double bias=0.0;
1861 5683 : if(grid_) {
1862 5122 : bias=BiasGrid_->getBias(cv);
1863 : } else {
1864 : #ifdef _OPENMP
1865 : //nt is unused if openmp is not declared
1866 561 : const unsigned nt=OpenMP::getNumThreads();
1867 : #endif //_OPENMP
1868 561 : unsigned stride=comm.Get_size();
1869 561 : unsigned rank=comm.Get_rank();
1870 : // with the neighbor list only the listed hills are summed
1871 561 : const unsigned nhills=(nlist_.enabled ? nlist_.index.size() : hills_.size());
1872 :
1873 561 : #pragma omp parallel num_threads(nt)
1874 : {
1875 : #pragma omp for reduction(+:bias) nowait
1876 : for(unsigned i=rank; i<nhills; i+=stride) {
1877 : bias+=evaluateHill(cv,nlist_.enabled ? nlist_.index[i] : i);
1878 : }
1879 : }
1880 561 : comm.Sum(bias);
1881 : }
1882 :
1883 5683 : return bias;
1884 : }
1885 :
1886 40512 : double MetaD::getBiasAndDerivatives(const std::vector<double>& cv, std::vector<double>& der) {
1887 40512 : unsigned ncv=getNumberOfArguments();
1888 : double bias=0.0;
1889 40512 : if(grid_) {
1890 22681 : bias=BiasGrid_->getBiasAndDerivatives(cv,der);
1891 : } else {
1892 17831 : unsigned nt=OpenMP::getNumThreads();
1893 17831 : unsigned stride=comm.Get_size();
1894 17831 : unsigned rank=comm.Get_rank();
1895 : // with the neighbor list only the listed hills are summed
1896 17831 : const unsigned nhills=(nlist_.enabled ? nlist_.index.size() : hills_.size());
1897 :
1898 : // NOTE: the choice between the serial and the parallel loop is based on the total
1899 : // number of hills also when the neighbor list is used
1900 17831 : if(hills_.size()<2*static_cast<std::size_t>(nt)*stride||nt==1) {
1901 : // for performance reasons and thread safety
1902 3498 : std::vector<double> dp(ncv);
1903 8158 : for(unsigned i=rank; i<nhills; i+=stride) {
1904 4660 : bias+=evaluateHillWithDerivatives(cv,nlist_.enabled ? nlist_.index[i] : i,der.data(),dp.data());
1905 : }
1906 : } else {
1907 14333 : ThreadWorkArrays<double> derwork(nt,ncv);
1908 14333 : ThreadWorkArrays<double> dpwork(nt,ncv);
1909 14333 : #pragma omp parallel num_threads(nt)
1910 : {
1911 : double* omp_deriv=derwork.get(OpenMP::getThreadNum());
1912 : // for performance reasons and thread safety
1913 : double* dp=dpwork.get(OpenMP::getThreadNum());
1914 : #pragma omp for reduction(+:bias) nowait
1915 : for(unsigned i=rank; i<nhills; i+=stride) {
1916 : bias+=evaluateHillWithDerivatives(cv,nlist_.enabled ? nlist_.index[i] : i,omp_deriv,dp);
1917 : }
1918 : #pragma omp critical
1919 : for(unsigned i=0; i<ncv; i++) {
1920 : der[i]+=omp_deriv[i];
1921 : }
1922 : }
1923 : }
1924 17831 : if(stride>1) {
1925 : // a single reduction for bias and derivatives
1926 850 : std::vector<double> buffer(ncv+1);
1927 850 : buffer[0]=bias;
1928 2550 : for(unsigned i=0; i<ncv; i++) {
1929 1700 : buffer[1+i]=der[i];
1930 : }
1931 850 : comm.Sum(buffer);
1932 850 : bias=buffer[0];
1933 2550 : for(unsigned i=0; i<ncv; i++) {
1934 1700 : der[i]=buffer[1+i];
1935 : }
1936 : }
1937 : }
1938 :
1939 40512 : return bias;
1940 : }
1941 :
1942 9439 : double MetaD::getHeight(const std::vector<double>& cv) {
1943 9439 : double height=height0_;
1944 9439 : if(welltemp_) {
1945 4581 : double vbias = getBias(cv);
1946 4581 : if(biasf_>1.0) {
1947 4565 : height = height0_*std::exp(-vbias/(kbt_*(biasf_-1.0)));
1948 : } else {
1949 : // notice that if gamma=1 we store directly -F
1950 16 : height = height0_*std::exp(-vbias/kbt_);
1951 : }
1952 : }
1953 9439 : if(dampfactor_>0.0) {
1954 18 : plumed_assert(BiasGrid_);
1955 18 : double m=BiasGrid_->getGrid().getMaxValue();
1956 18 : height*=std::exp(-m/(kbt_*(dampfactor_)));
1957 : }
1958 9439 : if (tt_specs_.is_active) {
1959 60 : double vbarrier = transition_bias_;
1960 60 : temperHeight(height, tt_specs_, vbarrier);
1961 : }
1962 9439 : if(TargetGrid_) {
1963 18 : if(!TargetGrid_->isInside(cv)) {
1964 0 : error("the variables are outside the range of the TARGET grid");
1965 : }
1966 18 : double f=TargetGrid_->getValue(cv)-TargetGrid_->getMaxValue();
1967 18 : height*=std::exp(f/kbt_);
1968 : }
1969 9439 : return height;
1970 : }
1971 :
1972 60 : void MetaD::temperHeight(double& height, const TemperingSpecs& t_specs, const double tempering_bias) {
1973 60 : if (t_specs.alpha == 1.0) {
1974 80 : height *= std::exp(-std::max(0.0, tempering_bias - t_specs.threshold) / (kbt_ * (t_specs.biasf - 1.0)));
1975 : } else {
1976 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));
1977 : }
1978 60 : }
1979 :
1980 40552 : void MetaD::calculate() {
1981 :
1982 40552 : schedule_.start(getStep());
1983 :
1984 40552 : const unsigned ncv=getNumberOfArguments();
1985 40552 : std::vector<double> cv(ncv);
1986 107926 : for(unsigned i=0; i<ncv; ++i) {
1987 67374 : cv[i]=getArgument(i);
1988 : }
1989 :
1990 : // during a replica exchange attempt the bias is computed on the coordinates of the partner.
1991 : // The MD code does not tell whether the exchange is accepted: with ADAPTIVE=DIFF this is
1992 : // detected at the next update, comparing the CVs with those of the partner (see updateFlexibleBin)
1993 40552 : if(adaptive_==FlexibleBin::diffusion && getExchangeStep()) {
1994 6 : exchange_.attempted=true;
1995 6 : exchange_.partner_cv=cv;
1996 : }
1997 :
1998 40552 : if(nlist_.enabled) {
1999 1643 : nlist_.steps++;
2000 1643 : if(getExchangeStep()) {
2001 0 : nlist_.update=true;
2002 : } else {
2003 1643 : const std::vector<double> position=getNlistPosition();
2004 3787 : for(unsigned i=0; i<ncv; ++i) {
2005 2469 : double d = difference(i, position[i], nlist_.center[i]);
2006 2469 : double nk_dist2 = d*d/nlist_.dev2[i];
2007 2469 : if(nk_dist2>nlist_.param[1]) {
2008 325 : nlist_.update=true;
2009 325 : break;
2010 : }
2011 : }
2012 : }
2013 1643 : if(nlist_.update) {
2014 856 : updateNlist();
2015 : }
2016 : }
2017 :
2018 : double ene = 0.;
2019 40552 : std::vector<double> der(ncv,0.);
2020 40552 : if(biasf_!=1.0) {
2021 40512 : ene = getBiasAndDerivatives(cv,der);
2022 : }
2023 40552 : setBias(ene);
2024 107926 : for(unsigned i=0; i<ncv; i++) {
2025 67374 : setOutputForce(i,-der[i]);
2026 : }
2027 :
2028 40552 : if(calc_work_) {
2029 1102 : getPntrToComponent("work")->set(work_);
2030 : }
2031 40552 : if(calc_rct_) {
2032 10048 : getPntrToComponent("rbias")->set(ene - reweight_factor_);
2033 : }
2034 : // calculate the acceleration factor
2035 : // (only once per STRIDE steps: the action can be calculated more often when other actions need it)
2036 41426 : if(acceleration_&&!schedule_.isRunFirstStep(getStep())&&!getExchangeStep()&&getStep()%getStride()==0) {
2037 874 : acc_ += static_cast<double>(getStride()) * std::exp(ene/(kbt_));
2038 874 : const double mean_acc = acc_/static_cast<double>(getStep());
2039 1748 : getPntrToComponent("acc")->set(mean_acc);
2040 39683 : } else if (acceleration_ && schedule_.isRunFirstStep(getStep()) && acc_restart_mean_ > 0.0) {
2041 2 : if(getStep()>acc_restart_step_ && !getExchangeStep() && getStep()%getStride()==0) {
2042 : // the previous run did not compute this step (the steps in between, if any, are assumed to have the mean value)
2043 1 : acc_ = acc_restart_mean_ * static_cast<double>(getStep()-getStride()) + static_cast<double>(getStride()) * std::exp(ene/(kbt_));
2044 2 : getPntrToComponent("acc")->set(acc_/static_cast<double>(getStep()));
2045 : } else {
2046 : // (with STRIDE, the acceleration is accumulated at the steps that are multiples of STRIDE)
2047 1 : acc_ = acc_restart_mean_ * static_cast<double>(getStep()-getStep()%getStride());
2048 : }
2049 2 : if(fa_.enabled) {
2050 : // has to be done here if restarting, as the acc is not defined before
2051 1 : updateFrequencyAdaptiveStride();
2052 : }
2053 : }
2054 40552 : }
2055 :
2056 40722 : bool MetaD::isHillStep() const {
2057 40722 : return schedule_.isHillStep(getStep(),current_stride_,getTime(),0.5*stride_*getTimeStep());
2058 : }
2059 :
2060 38856 : void MetaD::beforeUpdate() {
2061 : // called at every step, also outside the update window and when the action is not active
2062 : // (e.g. with STRIDE), in which case the arguments have not been computed
2063 38856 : if(precal_diff_ && isActive()) {
2064 546 : updateFlexibleBin(false);
2065 : }
2066 38856 : }
2067 :
2068 8488 : void MetaD::updateFlexibleBin(bool nowAddAHill) {
2069 8488 : if(adaptive_==FlexibleBin::diffusion && diffrestart_.pending) {
2070 12 : initializeDiffusionStatistics();
2071 : }
2072 8488 : if(adaptive_==FlexibleBin::diffusion) {
2073 5566 : const unsigned ncv=getNumberOfArguments();
2074 5566 : std::vector<double> cv(ncv);
2075 14912 : for(unsigned i=0; i<ncv; ++i) {
2076 9346 : cv[i]=getArgument(i);
2077 : }
2078 5566 : if(exchange_.attempted && exchange_.last_cv.size()==ncv) {
2079 : // the exchange has been accepted if the CVs are now closer to those of the partner than
2080 : // to their own value at the previous update (distances in units of the current fluctuations)
2081 6 : const std::vector<double> var=flexbin_->getMatrix();
2082 : double dpartner=0.0;
2083 : double down=0.0;
2084 : unsigned k=0;
2085 12 : for(unsigned i=0; i<ncv; i++) {
2086 6 : const double s2=(var.size()>k && var[k]>0.0 ? var[k] : 1.0);
2087 6 : const double dp=difference(i,exchange_.partner_cv[i],cv[i]);
2088 6 : const double dl=difference(i,exchange_.last_cv[i],cv[i]);
2089 6 : dpartner+=dp*dp/s2;
2090 6 : down+=dl*dl/s2;
2091 6 : k+=ncv-i; // index of the next diagonal element of the upper triangle
2092 : }
2093 6 : if(dpartner<down) {
2094 : // the CVs jumped: restart the running average, but keep the variance
2095 4 : flexbin_->restartAverage();
2096 : }
2097 : }
2098 5566 : exchange_.attempted=false;
2099 5566 : exchange_.last_cv=cv;
2100 : }
2101 8488 : flexbin_->update(nowAddAHill);
2102 8488 : }
2103 :
2104 12 : void MetaD::initializeDiffusionStatistics() {
2105 12 : diffrestart_.pending=false;
2106 12 : const double dt=getTimeStep();
2107 : // times are written with a finite precision (FMT): times closer than half a PACE are the same time.
2108 : // With MPI walkers each time appears once per walker, in the order of the replicas.
2109 12 : std::vector<double> average;
2110 12 : std::vector<double> variance;
2111 12 : double lasttime=0.0;
2112 24 : const unsigned n=diffrestart_.compute(0.5*stride_*dt,getTime(),mw_.mpi ? mw_.mpi_nw : 0,mw_.mpi_id,
2113 12 : 1.0/sigma0_[0],dt*getStride(),
2114 958 : [&](unsigned i,double a,double b) {
2115 958 : return difference(i,a,b);
2116 : },
2117 18 : [&](unsigned i,double x) {
2118 18 : return bringBackInPbc(i,x);
2119 : },
2120 : average,variance,lasttime);
2121 12 : diffrestart_.clear();
2122 12 : if(n==0) {
2123 3 : log.printf(" ADAPTIVE=DIFF: no hills available, the statistics for the hills width start from scratch\n");
2124 : return;
2125 : }
2126 9 : flexbin_->setStatistics(average,variance,2);
2127 9 : log.printf(" ADAPTIVE=DIFF: the statistics for the hills width are initialized from %u hills of %s (last one at time %f)\n",
2128 : n,diffrestart_.getFile().c_str(),lasttime);
2129 : }
2130 :
2131 37800 : void MetaD::update() {
2132 : // adding hills criteria (could be more complex though)
2133 37800 : bool nowAddAHill=isHillStep();
2134 :
2135 37800 : unsigned ncv=getNumberOfArguments();
2136 37800 : std::vector<double> cv(ncv);
2137 101872 : for(unsigned i=0; i<ncv; ++i) {
2138 64072 : cv[i] = getArgument(i);
2139 : }
2140 :
2141 : // (with BIASFACTOR=1 there is no bias: the hills contain -F)
2142 : double vbias=0.;
2143 37800 : if(calc_work_ && biasf_!=1.0) {
2144 551 : vbias=getBias(cv);
2145 : }
2146 :
2147 : // if you use adaptive, call the FlexibleBin
2148 : // (with PRECAL_DIFF the statistics have already been updated in beforeUpdate)
2149 : bool multivariate=false;
2150 37800 : if(adaptive_!=FlexibleBin::none) {
2151 8238 : if(!precal_diff_) {
2152 7942 : updateFlexibleBin(nowAddAHill);
2153 : }
2154 : multivariate=true;
2155 : }
2156 :
2157 : // with ADAPTIVE=DIFF the width of a hill requires the fluctuations of at least two CV values
2158 : // (this happens at the first update, e.g. after UPDATE_FROM or when restarting)
2159 37800 : if(nowAddAHill && adaptive_==FlexibleBin::diffusion && flexbin_->getNumberOfSamples()<2) {
2160 : nowAddAHill=false;
2161 : }
2162 :
2163 37795 : if(nowAddAHill) {
2164 9439 : addHill(cv,multivariate);
2165 : }
2166 :
2167 : // this should be outside of the if block in case
2168 : // WALKERS_RSTRIDE is not a multiple of PACE
2169 37800 : if(mw_.isReadStep(getStep())) {
2170 3012 : hillsOfile_.flush();
2171 : }
2172 :
2173 37800 : if(calc_work_ && biasf_!=1.0) {
2174 551 : if(nlist_.enabled) {
2175 0 : updateNlist();
2176 : }
2177 551 : double vbias1=getBias(cv);
2178 551 : work_+=vbias1-vbias;
2179 : }
2180 :
2181 : // dump grid on file
2182 37800 : if(wgridstride_>0&&(getStep()%wgridstride_==0||getCPT())) {
2183 98 : writeBiasGrid();
2184 : }
2185 :
2186 : // if multiple walkers and time to read Gaussians
2187 37800 : if(mw_.isReadStep(getStep())) {
2188 3012 : readOtherWalkersHills();
2189 : }
2190 :
2191 : // Recalculate special bias quantities whenever the bias has been changed by the update.
2192 37800 : updateBiasQuantities(nowAddAHill);
2193 :
2194 : // Frequency adaptive metadynamics - update hill addition frequency
2195 37800 : if(fa_.enabled && getStep()%fa_.update_frequency==0) {
2196 179 : updateFrequencyAdaptiveStride();
2197 : }
2198 37800 : }
2199 :
2200 9439 : void MetaD::addHill(const std::vector<double>& cv, bool multivariate) {
2201 9439 : unsigned ncv=getNumberOfArguments();
2202 9439 : std::vector<double> thissigma;
2203 : // add a Gaussian
2204 9439 : double height=getHeight(cv);
2205 : // returns upper diagonal inverse
2206 9439 : if(adaptive_!=FlexibleBin::none) {
2207 3908 : thissigma=flexbin_->getInverseMatrix();
2208 : }
2209 : // returns normal sigma
2210 : else {
2211 7485 : thissigma=sigma0_;
2212 : }
2213 :
2214 : // In case we use walkers_mpi, it is now necessary to communicate with other replicas.
2215 9439 : if(mw_.mpi) {
2216 : // Allocate arrays to store all walkers hills
2217 : // (packed in a single array: center, sigma, height, multivariate, so that one collective is enough)
2218 842 : const unsigned nsigma=thissigma.size();
2219 842 : const std::size_t recordsize=ncv+nsigma+2;
2220 842 : std::vector<double> record(recordsize);
2221 2526 : for(unsigned j=0; j<ncv; j++) {
2222 1684 : record[j]=cv[j];
2223 : }
2224 2602 : for(unsigned j=0; j<nsigma; j++) {
2225 1760 : record[ncv+j]=thissigma[j];
2226 : }
2227 : // notice that if gamma=1 we store directly -F so this scaling is not necessary:
2228 842 : record[ncv+nsigma]=height*(biasf_>1.0?biasf_/(biasf_-1.0):1.0);
2229 1608 : record[ncv+nsigma+1]=(multivariate ? 1.0 : 0.0);
2230 842 : std::vector<double> all_records(mw_.mpi_nw*recordsize,0.0);
2231 842 : if(comm.Get_rank()==0) {
2232 : // Communicate (only root)
2233 473 : multi_sim_comm.Allgather(record,all_records);
2234 : }
2235 : // Share info with group members
2236 842 : comm.Bcast(all_records,0);
2237 :
2238 : // Flying Gaussian
2239 842 : if (flying_) {
2240 54 : hills_.clear();
2241 54 : hillarrays_.clear();
2242 54 : comm.Barrier();
2243 : }
2244 :
2245 2700 : for(unsigned i=0; i<mw_.mpi_nw; i++) {
2246 : // actually add hills one by one
2247 1858 : const double* rec=&all_records[i*recordsize];
2248 1858 : std::vector<double> cv_now(rec,rec+ncv);
2249 1858 : std::vector<double> sigma_now(rec+ncv,rec+ncv+nsigma);
2250 : // notice that if gamma=1 we store directly -F so this scaling is not necessary:
2251 1858 : double fact=(biasf_>1.0?(biasf_-1.0)/biasf_:1.0);
2252 1858 : Gaussian newhill=Gaussian(rec[ncv+nsigma+1]!=0.0,rec[ncv+nsigma]*fact,cv_now,sigma_now);
2253 1858 : addGaussian(newhill);
2254 1858 : if(!flying_) {
2255 1696 : writeGaussian(newhill,hillsOfile_);
2256 : }
2257 1858 : }
2258 : } else {
2259 8597 : Gaussian newhill=Gaussian(multivariate,height,cv,thissigma);
2260 8597 : addGaussian(newhill);
2261 8597 : writeGaussian(newhill,hillsOfile_);
2262 8597 : }
2263 :
2264 : // this is to update the hills neighbor list
2265 9439 : if(nlist_.enabled) {
2266 684 : nlist_.update=true;
2267 : }
2268 9439 : }
2269 :
2270 98 : void MetaD::writeBiasGrid() {
2271 : if(mw_.writesFiles()) {
2272 : // the hills of the other walkers in the grid (see readHillsAtStart), in a comment after the grid
2273 294 : BiasGrid_->write(gridfile_,storeOldGrids_,mw_.n>1 ? hillsfiles_->getGridComment(mw_.id) : std::string());
2274 : }
2275 98 : }
2276 :
2277 3012 : void MetaD::readOtherWalkersHills() {
2278 3012 : bias::readOtherWalkersHills(*this,mw_,std::vector<MetaDHillsFiles*>(1,hillsfiles_.get()));
2279 : // this is to update the hills neighbor list
2280 3012 : if(nlist_.enabled) {
2281 0 : nlist_.update=true;
2282 : }
2283 3012 : }
2284 :
2285 37800 : void MetaD::updateBiasQuantities(bool nowAddAHill) {
2286 37800 : bool bias_has_changed = (nowAddAHill || mw_.isReadStep(getStep()));
2287 37800 : if (calc_rct_ && bias_has_changed && getStep()%(stride_*rct_ustride_)==0) {
2288 1001 : computeReweightingFactor();
2289 : }
2290 37800 : if (calc_max_bias_ && bias_has_changed) {
2291 109 : max_bias_ = BiasGrid_->getGrid().getMaxValue();
2292 218 : getPntrToComponent("maxbias")->set(max_bias_);
2293 : }
2294 37800 : if (calc_transition_bias_ && bias_has_changed) {
2295 478 : transition_bias_ = getTransitionBarrierBias();
2296 956 : getPntrToComponent("transbias")->set(transition_bias_);
2297 : }
2298 37800 : }
2299 :
2300 1001 : void MetaD::computeReweightingFactor() {
2301 1001 : if(biasf_==1.0) { // in this case we have no bias, so reweight factor is 0.0
2302 0 : getPntrToComponent("rct")->set(0.0);
2303 0 : return;
2304 : }
2305 :
2306 1001 : double Z_0=0; //proportional to the integral of exp(-beta*F)
2307 1001 : double Z_V=0; //proportional to the integral of exp(-beta*(F+V))
2308 1001 : double minusBetaF=biasf_/(biasf_-1.)/kbt_;
2309 1001 : double minusBetaFplusV=1./(biasf_-1.)/kbt_;
2310 1001 : if (biasf_==-1.0) { //non well-tempered case
2311 0 : minusBetaF=1./kbt_;
2312 0 : minusBetaFplusV=0;
2313 : }
2314 1001 : max_bias_=BiasGrid_->getGrid().getMaxValue(); //to avoid exp overflow
2315 :
2316 : // the integrals are over all the points of the grid (in its current range); the points of a sparse grid that are
2317 : // not stored have zero bias
2318 1001 : const unsigned rank=comm.Get_rank();
2319 1001 : const unsigned stride=comm.Get_size();
2320 1001 : Grid::index_t nstored=0;
2321 2002 : BiasGrid_->getGrid().forEachStoredPoint([&](Grid::index_t,double val) {
2322 10507134 : if(nstored%stride==rank) {
2323 5125168 : Z_0+=std::exp(minusBetaF*(val-max_bias_));
2324 5125168 : Z_V+=std::exp(minusBetaFplusV*(val-max_bias_));
2325 : }
2326 10507134 : nstored++;
2327 10507134 : });
2328 : double npoints=1.0;
2329 3003 : for(const auto n : BiasGrid_->getGrid().getNbin()) {
2330 2002 : npoints*=n;
2331 : }
2332 1001 : const double nzero=npoints-static_cast<double>(nstored);
2333 1001 : if(rank==0 && nzero>0.0) {
2334 218 : Z_0+=nzero*std::exp(-minusBetaF*max_bias_);
2335 218 : Z_V+=nzero*std::exp(-minusBetaFplusV*max_bias_);
2336 : }
2337 1001 : comm.Sum(Z_0);
2338 1001 : comm.Sum(Z_V);
2339 :
2340 1001 : reweight_factor_=kbt_*std::log(Z_0/Z_V)+max_bias_;
2341 2002 : getPntrToComponent("rct")->set(reweight_factor_);
2342 : }
2343 :
2344 491 : double MetaD::getTransitionBarrierBias() {
2345 : // If there is only one well of interest, return the bias at that well point.
2346 491 : if (transitionwells_.size() == 1) {
2347 0 : double tb_bias = getBias(transitionwells_[0]);
2348 0 : return tb_bias;
2349 :
2350 : // Otherwise, check for the least barrier bias between all pairs of wells.
2351 : // Note that because the paths can be considered edges between the wells' nodes
2352 : // to make a graph and the path barriers satisfy certain cycle inequalities, it
2353 : // is sufficient to look at paths corresponding to a minimal spanning tree of the
2354 : // overall graph rather than examining every edge in the graph.
2355 : // For simplicity, I chose the star graph with center well 0 as the spanning tree.
2356 : // It is most efficient to start the path searches from the wells that are
2357 : // expected to be sampled last, so transitionwell_[0] should correspond to the
2358 : // starting well. With this choice the searches will terminate in one step until
2359 : // transitionwell_[1] is sampled.
2360 : } else {
2361 : double least_transition_bias;
2362 491 : std::vector<double> sink = transitionwells_[0];
2363 491 : std::vector<double> source = transitionwells_[1];
2364 491 : least_transition_bias = BiasGrid_->getGrid().findMaximalPathMinimum(source, sink);
2365 491 : for (unsigned i = 2; i < transitionwells_.size(); i++) {
2366 0 : if (least_transition_bias == 0.0) {
2367 : break;
2368 : }
2369 0 : source = transitionwells_[i];
2370 0 : double curr_transition_bias = BiasGrid_->getGrid().findMaximalPathMinimum(source, sink);
2371 0 : least_transition_bias = fmin(curr_transition_bias, least_transition_bias);
2372 : }
2373 : return least_transition_bias;
2374 : }
2375 : }
2376 :
2377 183 : void MetaD::updateFrequencyAdaptiveStride() {
2378 183 : plumed_massert(fa_.enabled,"should only be used if frequency adaptive metadynamics is enabled");
2379 183 : plumed_massert(acceleration_,"frequency adaptive metadynamics can only be used if the acceleration factor is calculated");
2380 183 : if(getStep()==0) {
2381 : // no acceleration factor yet (e.g. when called from the constructor)
2382 5 : getPntrToComponent("pace")->set(current_stride_);
2383 5 : return;
2384 : }
2385 178 : const double mean_acc = acc_/static_cast<double>(getStep());
2386 178 : int tmp_stride= stride_*floor((mean_acc/fa_.min_acceleration)+0.5);
2387 178 : if(mean_acc >= fa_.min_acceleration) {
2388 156 : if(tmp_stride > current_stride_) {
2389 32 : current_stride_ = tmp_stride;
2390 : }
2391 : }
2392 178 : if(fa_.max_stride!=0 && current_stride_>fa_.max_stride) {
2393 26 : current_stride_=fa_.max_stride;
2394 : }
2395 356 : getPntrToComponent("pace")->set(current_stride_);
2396 : }
2397 :
2398 40558 : bool MetaD::checkNeedsGradients()const {
2399 : // with ADAPTIVE=GEOM gradients are needed when a hill is added (the pace can change with FREQUENCY_ADAPTIVE)
2400 40558 : return adaptive_==FlexibleBin::geometry && isHillStep();
2401 : }
2402 :
2403 2488 : std::vector<double> MetaD::getNlistPosition() const {
2404 2488 : std::vector<double> cv(getNumberOfArguments());
2405 6772 : for(unsigned i=0; i<cv.size(); i++) {
2406 4284 : cv[i]=getArgument(i);
2407 : }
2408 2488 : if(interval_.enabled) {
2409 1384 : cv[0]=std::min(std::max(cv[0],interval_.lower),interval_.upper);
2410 : }
2411 2488 : return cv;
2412 : }
2413 :
2414 856 : void MetaD::updateNlist() {
2415 : // no need to check for neighbors
2416 856 : if(hills_.size()==0) {
2417 11 : return;
2418 : }
2419 :
2420 : // here we generate the neighbor list:
2421 : // first the hills close to the current CVs are flagged, then they are collected in their
2422 : // original order, so that the list (and thus the bias) does not depend on the number of threads
2423 845 : std::vector<double> cv=getNlistPosition();
2424 845 : std::vector<char> inlist(hills_.size(),0); // char rather than bool, so that threads can write concurrently
2425 845 : unsigned nt=OpenMP::getNumThreads();
2426 845 : if(hills_.size()<2*nt) {
2427 : nt=1;
2428 : }
2429 845 : #pragma omp parallel for num_threads(nt)
2430 : for(unsigned k=0; k<hills_.size(); k++) {
2431 : inlist[k]=(kernel_.getDistance2(cv,hills_[k])<=nlist_.param[0]*dp2cutoff);
2432 : }
2433 845 : nlist_.index.clear();
2434 93750 : for(unsigned k=0; k<hills_.size(); k++) {
2435 92905 : if(inlist[k]) {
2436 77958 : nlist_.index.push_back(k);
2437 : }
2438 : }
2439 :
2440 : // here we set some properties that are used to decide when to update it again
2441 2389 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
2442 1544 : nlist_.center[i]=cv[i];
2443 : }
2444 845 : std::vector<double> dev2;
2445 845 : dev2.resize(getNumberOfArguments(),0);
2446 78803 : for(unsigned k=0; k<nlist_.index.size(); k++) {
2447 224443 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
2448 146485 : const double d=difference(i,cv[i],hills_[nlist_.index[k]].center[i]);
2449 146485 : dev2[i]+=d*d;
2450 : }
2451 : }
2452 2389 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
2453 1544 : if(dev2[i]>0.) {
2454 1544 : nlist_.dev2[i]=dev2[i]/static_cast<double>(nlist_.index.size());
2455 : } else {
2456 0 : nlist_.dev2[i]=kernel_.getVariance(hills_.back(),i);
2457 : }
2458 : }
2459 :
2460 : // we are done
2461 845 : getPntrToComponent("nlker")->set(nlist_.index.size());
2462 845 : getPntrToComponent("nlsteps")->set(nlist_.steps);
2463 845 : nlist_.steps=0;
2464 845 : nlist_.update=false;
2465 845 : }
2466 :
2467 : }
2468 : }
|
