Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-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 "MetaDGrid.h"
24 : #include "MetaDHills.h"
25 : #include "MetaDKernel.h"
26 : #include "MetaDRestart.h"
27 : #include "core/ActionRegister.h"
28 : #include "core/PlumedMain.h"
29 : #include "core/FlexibleBin.h"
30 : #include "tools/Exception.h"
31 : #include "tools/OpenMP.h"
32 : #include "tools/File.h"
33 : #include "tools/Communicator.h"
34 : #include <algorithm>
35 : #include <cmath>
36 : #include <limits>
37 :
38 : namespace PLMD {
39 : namespace bias {
40 :
41 : //+PLUMEDOC BIAS PBMETAD
42 : /*
43 : Used to performed Parallel Bias metadynamics.
44 :
45 : This action provides an implementation of the Parallel Bias Metadynamics (PBMetaD) method that was introduced in the first paper cited below.
46 : This variant of metadynamics (see [METAD](METAD.md)) uses multiple low-dimensional bias potentials that are applied in parallel.
47 : In the current implementation, these have the form of mono-dimensional metadynamics bias
48 : potentials:
49 :
50 : $$
51 : {V(s_1,t), ..., V(s_N,t)}
52 : $$
53 :
54 : where:
55 :
56 : $$
57 : V(s_i,t) = \sum_{ k \tau < t} W_i(k \tau)
58 : \exp\left(
59 : - \frac{(s_i-s_i^{(0)}(k \tau))^2}{2\sigma_i^2}
60 : \right).
61 : $$
62 :
63 : To ensure the convergence of each mono-dimensional bias potential to the corresponding free energy,
64 : at each deposition step the Gaussian heights are multiplied by the so-called conditional term:
65 :
66 : $$
67 : W_i(k \tau)=W_0 \frac{\exp\left(
68 : - \frac{V(s_i,k \tau)}{k_B T}
69 : \right)}{\sum_{i=1}^N
70 : \exp\left(
71 : - \frac{V(s_i,k \tau)}{k_B T}
72 : \right)}
73 : $$
74 :
75 : where $W_0$ is the initial Gaussian height.
76 :
77 : The PBMetaD bias potential is defined by:
78 :
79 : $$
80 : V_{PB}(\vec{s},t) = -k_B T \log{\frac{1}{N}\sum_{i=1}^N
81 : \exp\left(
82 : - \frac{V(s_i,t)}{k_B T}
83 : \right)}.
84 : $$
85 :
86 : PBMETAD needs the temperature: give it with TEMP, unless the MD engine passes it to PLUMED.
87 : As in [METAD](METAD.md), the Gaussians are truncated and stretched so that they go to zero at the cutoff
88 : (the hills files contain stretched Gaussians, which is what [sum_hills](sum_hills.md) assumes).
89 :
90 :
91 : Information on the Gaussian functions that build each bias potential are printed to
92 : multiple HILLS files, which
93 : are used both to restart the calculation and to reconstruct the mono-dimensional
94 : free energies as a function of the corresponding CVs.
95 : These can be reconstructed using the [sum_hills](sum_hills.md) utility because the final bias is given
96 : by:
97 :
98 : $$
99 : V(s_i) = -F(s_i)
100 : $$
101 :
102 : Currently, only a subset of the [METAD](METAD.md) options are available in PBMetaD.
103 :
104 : The bias potentials can be stored on a grid to increase performances of long PBMetaD simulations.
105 : You should
106 : provide either the number of bins for every collective variable (GRID_BIN) or
107 : the desired grid spacing (GRID_SPACING). In case you provide both PLUMED will use
108 : the most conservative choice (highest number of bins) for each dimension.
109 : In case you do not provide any information about bin size (neither GRID_BIN nor GRID_SPACING)
110 : and if Gaussian width is fixed PLUMED will use 1/5 of the Gaussian width as grid spacing.
111 : This default choice should be reasonable for most applications.
112 :
113 : As in [METAD](METAD.md), the grids are extended automatically, along the variables that are not periodic, whenever a
114 : Gaussian does not fit in them: a grid always contains all the points where a Gaussian is not zero, so that the bias is
115 : exactly zero outside it (with INTERVAL_MIN and INTERVAL_MAX, the grid contains the interval, and beyond the limits the bias
116 : is the constant value at the limits). The grid spacing does not change and the new points are on the same lattice as the initial ones.
117 : Every extension is written on the log. If you give GRID_SPACING without GRID_MIN and GRID_MAX each grid is placed
118 : around its first Gaussian and then extended as needed. GRID_SIZE_LIMIT sets the maximum number of points of each grid.
119 : A grid read with GRID_RFILES can have any range, but it must be on the lattice of the grid given in input.
120 :
121 : As in [METAD](METAD.md), no Gaussian is added at the first step of a simulation. When a simulation is continued with
122 : RESTART and the bias is read from the hills files, a Gaussian is added at the first step unless they already contain
123 : it (i.e. the previous run had already performed this step). When the bias is read from the grids (GRID_RFILES), the hills files are still read,
124 : without adding their hills to the bias, to know the time of their last hill. With WALKERS_N, the grids written by a walker
125 : contain also the hills it has read from the files of the other walkers (their number is written in a comment line after
126 : the grid), and only the hills of the other walkers that are not in the grids are added.
127 :
128 : With ADAPTIVE=DIFF, as in [METAD](METAD.md), when restarting the statistics of the fluctuations of each variable
129 : (running average and variance) are initialized from its last hills in the hills file of its family, so that the hills
130 : deposited after the restart have the same width they would have had without the restart (use DIFF_FROM_SCRATCH to
131 : accumulate the statistics from scratch). With replica exchange, the running averages of the variables are
132 : restarted when an exchange is accepted, and the variances are kept.
133 : With PRECAL_DIFF the statistics are accumulated at every step, also outside the update window (UPDATE_FROM,
134 : UPDATE_UNTIL), so that the first hills have a width estimated over a sufficiently long time.
135 :
136 : With GRID_WSTRIDE the grids are written to the GRID_WFILES, overwriting the previous ones, or, with STORE_GRIDS,
137 : appending them one after the other.
138 :
139 : Another option that is available is the well-tempered metadynamics that is discussed in the second paper cited below. In this
140 : variant of PBMetaD the heights of the Gaussian hills are scaled at each step by the
141 : additional well-tempered metadynamics term.
142 : This ensures that each bias converges more smoothly. It should be noted that, in the case of well-tempered metadynamics, in
143 : the output printed the Gaussian height is re-scaled using the bias factor.
144 : Also notice that with well-tempered metadynamics the HILLS files do not contain the bias,
145 : but the negative of the free-energy estimate. This choice has the advantage that
146 : one can restart a simulation using a different value for the $\Delta T$. The applied bias will be scaled accordingly.
147 : As in [METAD](METAD.md), BIASFACTOR can also be set equal to 1: in this case no bias is applied, and the hills files and
148 : the grids contain the estimate of the negative of the free energy (the heights of the hills are those of the limit
149 : of well-tempered metadynamics for a bias factor going to 1). Instead of HEIGHT, one can give TAU.
150 :
151 : Note that you can also use here the flexible Gaussian approach that is discussed in the third paper cited below.
152 : This method allows you to adapt the Gaussian to the extent of Cartesian space covered by a variable or
153 : to the space in collective variable covered in a given time. In this case the width of the deposited
154 : Gaussian potential is denoted by one value only that is a Cartesian space (ADAPTIVE=GEOM) or a time
155 : (ADAPTIVE=DIFF), and SIGMA_MIN and SIGMA_MAX (one value for each family) limit the widths. Note that a specific
156 : integration technique for the deposited Gaussian kernels should be used in this case. Check the documentation for utility
157 : sum_hills.
158 :
159 : With the keywords INTERVAL_MIN and INTERVAL_MAX (one value for each family; the values for periodic variables are
160 : ignored) one changes the metadynamics algorithm and sets the bias force equal to zero
161 : outside the boundaries as discussed in the fourth paper cited below. If, for example, metadynamics is performed on a CV s and one is interested only
162 : to the free energy for s > boundary, the history dependent potential is still updated according to the above
163 : equations but the metadynamics force is set to zero for s < boundary. Notice that Gaussian kernels are added also
164 : if s < boundary, as the tails of these Gaussian kernels influence VG in the relevant region s > boundary. In this way, the
165 : force on the system in the region s > boundary comes from both metadynamics and the force field, in the region
166 : s < boundary only from the latter. This approach allows obtaining a history-dependent bias potential VG that
167 : fluctuates around a stable estimator, equal to the negative of the free energy far enough from the
168 : boundaries. Note that:
169 :
170 : - It works only for one-dimensional biases;
171 : - It works both with and without GRID;
172 : - The interval limit boundary should be set in a region where the free energy derivative is not large;
173 : - If in the region outside the limit boundary the system has a free energy minimum, INTERVAL_MIN and INTERVAL_MAX should
174 : be used together with a [UPPER_WALLS](UPPER_WALLS.md) or [LOWER_WALLS](LOWER_WALLS.md) at boundary.
175 :
176 : For systems with multiple CVs that share identical properties, PBMetaD with partitioned families can be used
177 : to group them under one bias potential, to which each of them contributes, as discussed in the penultimate paper in the
178 : references below. This is done with a list of numbered PF keywords (PF0, PF1, ...), each of which contains the list of
179 : CVs from ARG to be placed in that family. Once invoked, each CV in ARG must be placed in exactly one PF, even if it
180 : results in families containing only one CV, and the CVs in a family must have the same periodicity.
181 : Keywords that take one value for each argument without families (SIGMA, SIGMA_MIN, SIGMA_MAX, FILE, the GRID keywords,
182 : INTERVAL_MIN, INTERVAL_MAX) take one value for each family instead. The hills of all the CVs of a family are written in
183 : the hills file of the family, as hills of its first CV.
184 :
185 : With SELECTOR, the bias is applied and the hills are added only when the value of the selector (see
186 : [SELECTOR](SELECTOR.md)) is equal to SELECTOR_ID.
187 :
188 : The multiple walkers that is discussed in the final paper cited below can also be used. See below the examples.
189 :
190 : ## Examples
191 :
192 : The following input is for PBMetaD calculation using as
193 : collective variables the distance between atoms 3 and 5
194 : and the distance between atoms 2 and 4. The value of the CVs and
195 : the PBMetaD bias potential are written to the COLVAR file every 100 steps.
196 :
197 : ```plumed
198 : DISTANCE ATOMS=3,5 LABEL=d1
199 : DISTANCE ATOMS=2,4 LABEL=d2
200 : PBMETAD ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3 PACE=500 TEMP=300 LABEL=pb FILE=HILLS_d1,HILLS_d2
201 : PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
202 : ```
203 :
204 : (See also [DISTANCE](DISTANCE.md) and [PRINT](PRINT.md)).
205 :
206 : If you use well-tempered metadynamics, you should specify a single bias factor and initial Gaussian height.
207 :
208 : ```plumed
209 : DISTANCE ATOMS=3,5 LABEL=d1
210 : DISTANCE ATOMS=2,4 LABEL=d2
211 : PBMETAD ...
212 : ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3
213 : PACE=500 BIASFACTOR=8 TEMP=300 LABEL=pb
214 : FILE=HILLS_d1,HILLS_d2
215 : ... PBMETAD
216 : PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
217 : ```
218 :
219 : Using partitioned families, each CV in ARG must be in exactly one family. Here,
220 : the distance between atoms 1,2 is degenerate with 2,4, but not with the
221 : distance between 3,5. Note that two SIGMA are provided to match the two PF.
222 :
223 : ```plumed
224 : DISTANCE ATOMS=3,5 LABEL=d1
225 : DISTANCE ATOMS=2,4 LABEL=d2
226 : DISTANCE ATOMS=1,2 LABEL=d3
227 : PBMETAD ...
228 : ARG=d1,d2,d3 SIGMA=0.2,0.2 HEIGHT=0.3
229 : PF0=d1 PF1=d2,d3
230 : PACE=500 BIASFACTOR=8 TEMP=300 LABEL=pb
231 : FILE=HILLS_d1,HILLS_d2
232 : ... PBMETAD
233 : PRINT ARG=d1,d2,d3,pb.bias STRIDE=100 FILE=COLVAR
234 : ```
235 :
236 : The following input enables the MPI version of multiple-walkers.
237 :
238 : ```plumed
239 : DISTANCE ATOMS=3,5 LABEL=d1
240 : DISTANCE ATOMS=2,4 LABEL=d2
241 : PBMETAD ...
242 : ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3
243 : PACE=500 BIASFACTOR=8 TEMP=300 LABEL=pb
244 : FILE=HILLS_d1,HILLS_d2
245 : WALKERS_MPI
246 : ... PBMETAD
247 : PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
248 : ```
249 :
250 : The disk version of multiple-walkers can be
251 : enabled by setting the number of walker used, the id of the
252 : current walker which interprets the input file, the directory where the
253 : hills containing files resides, and the frequency to read the other walkers.
254 : Here is an example
255 :
256 : ```plumed
257 : DISTANCE ATOMS=3,5 LABEL=d1
258 : DISTANCE ATOMS=2,4 LABEL=d2
259 : PBMETAD ...
260 : ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3
261 : PACE=500 BIASFACTOR=8 TEMP=300 LABEL=pb
262 : FILE=HILLS_d1,HILLS_d2
263 : WALKERS_N=10
264 : WALKERS_ID=3
265 : WALKERS_DIR=../
266 : WALKERS_RSTRIDE=100
267 : ... PBMETAD
268 : PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
269 : ```
270 :
271 : where WALKERS_N is the total number of walkers, WALKERS_ID is the
272 : id of the present walker (starting from 0 ) and the WALKERS_DIR is the directory
273 : where all the walkers are located. WALKERS_RSTRIDE is the number of step between
274 : one update and the other.
275 :
276 : */
277 : //+ENDPLUMEDOC
278 :
279 : class PBMetaD : public Bias {
280 :
281 : private:
282 : using Gaussian = MetaDGaussian;
283 : /// A partitioned family: its arguments share a one dimensional bias, with its hills file and grid
284 : /// (without PF, there is one family per argument)
285 : struct Family {
286 : /// the first argument of the family, used for the hills file and the grid
287 : Value* arg=nullptr;
288 : /// evaluation of the hills (periodicity of the argument, INTERVAL)
289 : MetaDKernel kernel;
290 : /// the hills, without grid
291 : std::vector<Gaussian> hills;
292 : /// the hills, without grid, in contiguous arrays (for speed): center, inverse width (0 for multivariate hills),
293 : /// height, and whether they can be evaluated with these arrays (not multivariate, no INTERVAL)
294 : std::vector<double> center;
295 : std::vector<double> invsigma;
296 : std::vector<double> height;
297 : std::vector<char> diagonal;
298 : /// the value of hill h at cv (one dimensional), and its derivative (added to der) if der is given
299 806003 : double evaluate(const std::vector<double>& cv, unsigned h, double* der) const {
300 : double dp[1];
301 806003 : if(diagonal[h]) {
302 365484 : return der ? kernel.evaluateDiagonalWithDerivatives(cv,¢er[h],&invsigma[h],height[h],der,dp) :
303 66858 : kernel.evaluateDiagonal(cv,¢er[h],&invsigma[h],height[h]);
304 : }
305 440519 : return der ? kernel.evaluateWithDerivatives(cv,hills[h],der,dp) : kernel.evaluate(cv,hills[h]);
306 : }
307 : /// the bias, with grid
308 : std::unique_ptr<MetaDBiasGrid> grid;
309 : /// INTERVAL
310 : bool interval=false;
311 : double lower=0.0;
312 : double upper=0.0;
313 : /// hills files that are read (own file when restarting, files of the other walkers), and writing of the hills
314 : std::unique_ptr<MetaDHillsFiles> hillsfiles;
315 : std::unique_ptr<OFile> hillsOfile;
316 : std::unique_ptr<OFile> gridOfile;
317 : /// ADAPTIVE=DIFF: the last hills of the own hills file, for the statistics of the fluctuations at restart
318 : MetaDDiffusionRestart diffrestart;
319 : /// the arguments in the family
320 : std::vector<unsigned> args;
321 : };
322 : // general setup
323 : double kbt_=0.0;
324 : int stride_=0;
325 : // well-tempered MetaD
326 : bool welltemp_=false;
327 : /// BIASFACTOR (1 without well-tempering)
328 : double biasf_=1.0;
329 : /// BIASFACTOR=1: no bias is applied, the hills (which contain -F) are added with the heights of the limit of
330 : /// well-tempered metadynamics for a bias factor going to 1
331 : bool nobias_=false;
332 : /// ADAPTIVE=DIFF: accumulate the CV statistics at every step, also outside the update window (PRECAL_DIFF)
333 : bool precal_diff_=false;
334 : /// append every grid that is written (STORE_GRIDS)
335 : bool storeOldGrids_=false;
336 : // output files format
337 : std::string fmt_;
338 : /// the steps when hills are added (the first step of the run, the last hill read when restarting)
339 : MetaDHillSchedule schedule_;
340 : // Gaussian starting parameters
341 : double height0_=std::numeric_limits<double>::max();
342 : std::vector<double> sigma0_;
343 : std::vector<double> sigma0min_;
344 : std::vector<double> sigma0max_;
345 : int adaptive_=FlexibleBin::none;
346 : /// with ADAPTIVE, the widths of the hills of each argument
347 : std::vector<FlexibleBin> flexbin_;
348 : // Grids
349 : bool grid_=false;
350 : int wgridstride_=0;
351 : // Partitioned Families
352 : std::vector<Family> families_;
353 : /// the family of each argument
354 : std::vector<unsigned> pfs_;
355 : /// if partitioned families are enabled
356 : bool do_pf_=false;
357 : /// multiple walkers, sharing hills either through files (WALKERS_N) or through MPI (WALKERS_MPI)
358 : MetaDWalkers mw_;
359 : /// ADAPTIVE=DIFF and replica exchange: an exchange has been attempted at the last step, with the CVs of
360 : /// the partner (partner_cv) and the CVs at the previous update (last_cv)
361 98 : struct {
362 : bool attempted=false;
363 : std::vector<double> partner_cv;
364 : std::vector<double> last_cv;
365 : } exchange_;
366 : // variable for selector
367 : std::string selector_;
368 : bool do_select_=false;
369 : unsigned select_value_=0;
370 : unsigned current_value_=0;
371 :
372 : /// ADAPTIVE=DIFF: initialize the statistics of the fluctuations from the hills read when restarting
373 : void initializeDiffusionStatistics();
374 : /// update the adaptive hill widths; with ADAPTIVE=DIFF also handle replica exchanges
375 : void updateFlexibleBins(bool nowAddAHill);
376 : /// true if a hill has to be added at the current step
377 : bool isHillStep() const;
378 : void writeGaussian(unsigned f, const Gaussian&);
379 : void addGaussian(unsigned f, const Gaussian&);
380 : /// the bias of each argument, and its derivative if der is given
381 : void getBiases(const std::vector<double>& cv, std::vector<double>& bias, std::vector<double>* der);
382 :
383 : public:
384 : explicit PBMetaD(const ActionOptions&);
385 : void calculate() override;
386 : void update() override;
387 : void beforeUpdate() override;
388 : static void registerKeywords(Keywords& keys);
389 : bool checkNeedsGradients()const override;
390 : };
391 :
392 : PLUMED_REGISTER_ACTION(PBMetaD,"PBMETAD")
393 :
394 99 : void PBMetaD::registerKeywords(Keywords& keys) {
395 99 : Bias::registerKeywords(keys);
396 99 : keys.add("compulsory","SIGMA","the widths of the Gaussian hills");
397 99 : keys.add("compulsory","PACE","the frequency for hill addition, one for all biases");
398 99 : keys.add("optional","FILE","files in which the lists of added hills are stored, default names are assigned using arguments if FILE is not found");
399 99 : keys.add("optional","HEIGHT","the height of the Gaussian hills, one for all biases. Compulsory unless TAU, TEMP and BIASFACTOR are given");
400 99 : keys.add("optional","FMT","specify format for HILLS files (useful for decrease the number of digits in regtests)");
401 99 : keys.add("optional","BIASFACTOR","use well tempered metadynamics with this bias factor, one for all biases. With BIASFACTOR=1 no bias is applied, and the hills files and grids contain the estimate of -F");
402 99 : keys.add("optional","TEMP","the system temperature - it is needed unless the MD engine passes it to PLUMED");
403 99 : keys.add("optional","TAU","in well tempered metadynamics, sets height to (k_B Delta T*pace*timestep)/tau");
404 99 : keys.add("optional","GRID_MIN","the lower bounds for the grid (initial range: the grid is extended when needed)");
405 99 : keys.add("optional","GRID_MAX","the upper bounds for the grid (initial range: the grid is extended when needed)");
406 99 : keys.add("optional","GRID_BIN","the number of bins for the grid");
407 99 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN). "
408 : "Without GRID_MIN and GRID_MAX it gives the grid spacing, and the grid is placed around the first Gaussian");
409 99 : keys.add("compulsory","GRID_SIZE_LIMIT","100000000","the maximum number of points of each grid, which is extended when a Gaussian does not fit in it");
410 99 : keys.addFlag("GRID_SPARSE",false,"use a sparse grid to store hills");
411 99 : keys.addFlag("GRID_NOSPLINE",false,"don't use spline interpolation with grids");
412 99 : keys.addFlag("STORE_GRIDS",false,"append every grid that is written to GRID_WFILES, instead of overwriting the previous one");
413 99 : keys.add("optional","GRID_WSTRIDE", "frequency for dumping the grid");
414 99 : keys.add("optional","GRID_WFILES", "dump grid for the bias (one file for each family; METAD: GRID_WFILE), default names are used if GRID_WSTRIDE is used without GRID_WFILES.");
415 99 : keys.add("optional","GRID_RFILES", "read grid for the bias (one file for each family; METAD: GRID_RFILE)");
416 99 : keys.add("optional","ADAPTIVE","use a geometric (=GEOM) or diffusion (=DIFF) based hills width scheme. Sigma is one number that has distance units or timestep dimensions");
417 99 : keys.add("optional","SIGMA_MAX","the upper bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
418 99 : keys.add("optional","SIGMA_MIN","the lower bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
419 99 : 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");
420 99 : 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 files");
421 198 : keys.addInputKeyword("numbered","PF", "scalar", "specify which CVs belong in a partitioned family. Once a PF is specified, all CVs in ARG must be placed in a PF even if there is one CV per PF");
422 99 : keys.add("optional","SELECTOR", "add forces and do update based on the value of SELECTOR");
423 99 : keys.add("optional","SELECTOR_ID", "value of SELECTOR");
424 99 : MetaDWalkers::registerKeywords(keys);
425 99 : keys.add("optional","INTERVAL_MIN","one dimensional lower limits (one for each family), outside the limits the system will not feel the biasing force (METAD: INTERVAL).");
426 99 : keys.add("optional","INTERVAL_MAX","one dimensional upper limits (one for each family), outside the limits the system will not feel the biasing force (METAD: INTERVAL).");
427 99 : keys.use("RESTART");
428 99 : keys.use("UPDATE_FROM");
429 99 : keys.use("UPDATE_UNTIL");
430 99 : keys.addDOI("10.1021/acs.jctc.5b00846");
431 99 : keys.addDOI("10.1103/PhysRevLett.100.020603");
432 99 : keys.addDOI("10.1021/ct3002464");
433 99 : keys.addDOI("10.2174/1877946811202010079");
434 99 : keys.addDOI("10.1021/acs.jctc.8b00448");
435 99 : keys.addDOI("10.1021/jp054359r");
436 99 : }
437 :
438 95 : PBMetaD::PBMetaD(const ActionOptions& ao):
439 95 : PLUMED_BIAS_INIT(ao) {
440 :
441 : // parse the flexible hills
442 : std::string adaptiveoption;
443 : adaptiveoption="NONE";
444 95 : parse("ADAPTIVE",adaptiveoption);
445 95 : if(adaptiveoption=="GEOM") {
446 1 : log.printf(" Uses Geometry-based hills width: sigma must be in distance units and only one sigma is needed\n");
447 1 : adaptive_=FlexibleBin::geometry;
448 94 : } else if(adaptiveoption=="DIFF") {
449 15 : log.printf(" Uses Diffusion-based hills width: sigma must be in time steps and only one sigma is needed\n");
450 15 : adaptive_=FlexibleBin::diffusion;
451 79 : } else if(adaptiveoption=="NONE") {
452 79 : adaptive_=FlexibleBin::none;
453 : } else {
454 0 : error("I do not know this type of adaptive scheme");
455 : }
456 :
457 95 : parse("FMT",fmt_);
458 :
459 : // Partitioned Families - fill with -1 to mark as invalid
460 98 : std::vector<int> pfs(getNumberOfArguments(), -1);
461 : // the first argument of each family (for HILLS & GRID files)
462 95 : std::vector<Value*> pfhold;
463 95 : std::vector<Value*> familyargs;
464 95 : for(int i = 0;; i++) {
465 244 : parseArgumentList("PF", i, familyargs);
466 122 : if (familyargs.empty()) {
467 : break;
468 : }
469 :
470 27 : do_pf_ = true;
471 27 : pfhold.push_back(nullptr);
472 27 : log << " Identified Partitioned Family " << i << ":";
473 63 : for (unsigned j = 0; j < familyargs.size(); j++) {
474 36 : log << " " << familyargs[j]->getName();
475 : // loop through the argument list to make sure it exists and assign it
476 : bool foundArg = false;
477 132 : for (unsigned argnum = 0; argnum < getNumberOfArguments(); argnum++) {
478 96 : if (familyargs[j]->getName() == getPntrToArgument(argnum)->getName()) {
479 : foundArg = true;
480 36 : if (pfs[argnum] != -1) {
481 0 : error(familyargs[j]->getName() + " already present in PF" + std::to_string(pfs[argnum]));
482 : }
483 36 : pfs[argnum] = i; // store the pf# for each cv
484 36 : if (pfhold[i] == nullptr) {
485 : // if this is the first argument in the family, store a pointer for it (this is for HILLS & GRID files)
486 27 : pfhold[i] = getPntrToArgument(argnum);
487 : }
488 : }
489 : }
490 36 : if (!foundArg) {
491 0 : error(familyargs[j]->getName() + " in PF" + std::to_string(i) + " not found in ARG");
492 : }
493 : }
494 27 : log << "\n";
495 27 : }
496 :
497 : // if PF were specified, every argument gets treated as its own PF
498 95 : if (!do_pf_) {
499 250 : for(unsigned i=0; i < getNumberOfArguments(); i++) {
500 169 : pfhold.push_back(getPntrToArgument(i));
501 169 : pfs[i] = i;
502 : }
503 : } else {
504 : // If we are doing PF, make sure each argument got assigned to a family.
505 50 : for (unsigned i = 0; i < getNumberOfArguments(); i++) {
506 36 : if (pfs[i] == -1) {
507 0 : error(getPntrToArgument(i)->getName() + " was not assigned a PF");
508 : }
509 : }
510 : }
511 95 : const unsigned pf_n=pfhold.size();
512 95 : pfs_.assign(pfs.begin(),pfs.end());
513 : // the arguments of a family share a one dimensional bias (and its hills file and grid)
514 299 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
515 : const Value* arg=getPntrToArgument(i);
516 205 : const Value* first=pfhold[pfs_[i]];
517 205 : bool same=(arg->isPeriodic()==first->isPeriodic());
518 205 : if(same && arg->isPeriodic()) {
519 : std::string amin,amax,fmin,fmax;
520 68 : arg->getDomain(amin,amax);
521 68 : first->getDomain(fmin,fmax);
522 68 : same=(amin==fmin && amax==fmax);
523 : }
524 205 : if(!same) {
525 4 : error("the arguments in a partitioned family must have the same periodicity: "+arg->getName()+" and "+first->getName()+" do not");
526 : }
527 : }
528 94 : families_.resize(pf_n);
529 297 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
530 203 : families_[pfs_[i]].args.push_back(i);
531 : }
532 289 : for(unsigned f=0; f<pf_n; f++) {
533 195 : families_[f].arg=pfhold[f];
534 195 : families_[f].kernel.setArguments(std::vector<Value*>(1,pfhold[f]));
535 : // stretched Gaussians (they go to zero at the cutoff), as in METAD and as written in the hills files
536 195 : if(!dp2cutoffNoStretch()) {
537 : families_[f].kernel.setStretch(dp2cutoffA,dp2cutoffB);
538 : }
539 : }
540 :
541 : // parse the sigma
542 188 : parseVector("SIGMA",sigma0_);
543 271 : for(const auto s : sigma0_) {
544 177 : if(!(s>0.0)) {
545 0 : error("SIGMA should be positive");
546 : }
547 : }
548 94 : if(adaptive_==FlexibleBin::none) {
549 : // if you use normal sigma you need one sigma per family
550 78 : if( sigma0_.size()!=pf_n ) {
551 0 : std::string fields = do_pf_ ? "PFs" : "arguments";
552 0 : error("number of " + fields + " does not match number of SIGMA parameters");
553 : }
554 78 : std::vector<double> smin,smax;
555 78 : parseVector("SIGMA_MIN",smin);
556 156 : parseVector("SIGMA_MAX",smax);
557 78 : if(smin.size()>0 || smax.size()>0) {
558 0 : error("SIGMA_MIN and SIGMA_MAX can only be used with ADAPTIVE");
559 : }
560 : } else {
561 : // if you use flexible hills you need one sigma
562 16 : if(sigma0_.size()!=1) {
563 0 : error("If you choose ADAPTIVE you need only one sigma according to your choice of type (GEOM/DIFF)");
564 : }
565 : // if adaptive then the number must be an integer
566 16 : if(adaptive_==FlexibleBin::diffusion) {
567 15 : if(int(sigma0_[0])-sigma0_[0]>1.e-9 || int(sigma0_[0])-sigma0_[0] <-1.e-9 || int(sigma0_[0])<1 ) {
568 0 : error("In case of adaptive hills with diffusion, the sigma must be an integer which is the number of time steps\n");
569 : }
570 : }
571 : // the limits of the widths (one value for each family)
572 32 : parseVector("SIGMA_MIN",sigma0min_);
573 16 : if(sigma0min_.size()>0 && sigma0min_.size()!=pf_n) {
574 0 : error("the number of SIGMA_MIN values should be the same as the number of arguments (PFs)");
575 16 : } else if(sigma0min_.size()==0) {
576 9 : sigma0min_.assign(pf_n,-1.);
577 : }
578 :
579 32 : parseVector("SIGMA_MAX",sigma0max_);
580 16 : if(sigma0max_.size()>0 && sigma0max_.size()!=pf_n) {
581 0 : error("the number of SIGMA_MAX values should be the same as the number of arguments (PFs)");
582 16 : } else if(sigma0max_.size()==0) {
583 15 : sigma0max_.assign(pf_n,-1.);
584 : }
585 :
586 : // one for each argument, with the limits of its family
587 54 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
588 38 : std::vector<double> tmp_smin(1,sigma0min_[pfs_[i]]);
589 38 : std::vector<double> tmp_smax(1,sigma0max_[pfs_[i]]);
590 76 : flexbin_.push_back(FlexibleBin(adaptive_,this,i,sigma0_[0],tmp_smin,tmp_smax));
591 : }
592 : }
593 94 : parseFlag("PRECAL_DIFF",precal_diff_);
594 94 : if(precal_diff_) {
595 1 : if(adaptive_!=FlexibleBin::diffusion) {
596 0 : error("PRECAL_DIFF can only be used with ADAPTIVE=DIFF");
597 : }
598 1 : log.printf(" The statistics for the hills width are accumulated at every step, also outside the update window (PRECAL_DIFF)\n");
599 : }
600 94 : bool diff_from_scratch=false;
601 94 : parseFlag("DIFF_FROM_SCRATCH",diff_from_scratch);
602 94 : if(diff_from_scratch && adaptive_!=FlexibleBin::diffusion) {
603 0 : error("DIFF_FROM_SCRATCH can only be used with ADAPTIVE=DIFF");
604 : }
605 : // the hills are kept only if needed (restart with ADAPTIVE=DIFF), see MetaDHillsFiles::readHills
606 289 : for(auto & family : families_) {
607 197 : family.diffrestart.pending=(adaptive_==FlexibleBin::diffusion && getRestart() && !diff_from_scratch);
608 : }
609 94 : if(adaptive_==FlexibleBin::diffusion && getRestart()) {
610 3 : if(diff_from_scratch) {
611 1 : log.printf(" The statistics for the hills width are accumulated from scratch after the restart (DIFF_FROM_SCRATCH)\n");
612 : } else {
613 2 : log.printf(" The statistics for the hills width are initialized from the last hills of the hills files, if available\n");
614 : }
615 : }
616 :
617 : // note: HEIGHT is not compulsory, since one could use the TAU keyword, see below
618 94 : parse("HEIGHT",height0_);
619 94 : if(height0_<0.0) {
620 0 : error("HEIGHT should not be negative");
621 : }
622 94 : parse("PACE",stride_);
623 94 : if(stride_<=0) {
624 3 : error("frequency for hill addition is nonsensical");
625 : }
626 :
627 94 : std::vector<std::string> hillsfname;
628 188 : parseVector("FILE",hillsfname);
629 94 : if(hillsfname.size()==0) {
630 36 : for(unsigned i=0; i< pf_n; i++) {
631 24 : std::string hillName = do_pf_ ? "HILLS.PF"+std::to_string(i) : "HILLS."+getPntrToArgument(i)->getName();
632 24 : hillsfname.push_back(hillName);
633 : }
634 : }
635 94 : if( hillsfname.size()!=pf_n ) {
636 0 : error("the number of FILE names does not match the number of arguments (PFs)");
637 : }
638 :
639 : // (-1: not given)
640 94 : double biasf=-1.0;
641 94 : parse("BIASFACTOR",biasf);
642 94 : if(biasf!=-1.0) {
643 85 : if( biasf<1.0 ) {
644 0 : error("well tempered bias factor is nonsensical");
645 : }
646 85 : biasf_=biasf;
647 85 : welltemp_=(biasf_>1.0);
648 85 : nobias_=(biasf_==1.0);
649 : }
650 94 : kbt_=getkBT();
651 : // the biases of the families are combined, and the heights of the hills are computed, with kBT
652 94 : if(kbt_==0.0) {
653 1 : error("Unless the MD engine passes the temperature to plumed, you must specify it using TEMP");
654 : }
655 93 : double tau=0.0;
656 93 : parse("TAU",tau);
657 93 : if(tau<0.0) {
658 1 : error("TAU should be positive");
659 : }
660 92 : if(tau==0.0) {
661 91 : if(height0_==std::numeric_limits<double>::max()) {
662 0 : error("At least one between HEIGHT and TAU should be specified");
663 : }
664 : // if tau is not set, we compute it here from the other input parameters
665 91 : if(welltemp_) {
666 82 : tau=(kbt_*(biasf_-1.0))/height0_*getTimeStep()*stride_;
667 : }
668 : } else {
669 1 : if(!welltemp_ && !nobias_) {
670 0 : error("TAU only makes sense in well-tempered metadynamics");
671 : }
672 1 : if(height0_!=std::numeric_limits<double>::max()) {
673 0 : error("At most one between HEIGHT and TAU should be specified");
674 : }
675 1 : if(welltemp_) {
676 0 : height0_=(kbt_*(biasf_-1.0))/tau*getTimeStep()*stride_;
677 : } else {
678 1 : height0_=kbt_/tau*getTimeStep()*stride_; // special case for BIASFACTOR=1, as in METAD
679 : }
680 : }
681 :
682 : // Multiple walkers
683 92 : mw_.parse(*this);
684 :
685 : // Grid file
686 92 : parse("GRID_WSTRIDE",wgridstride_);
687 92 : if(wgridstride_<0) {
688 0 : error("GRID_WSTRIDE should not be negative");
689 : }
690 92 : std::vector<std::string> gridfilenames;
691 92 : parseVector("GRID_WFILES",gridfilenames);
692 92 : if (wgridstride_ == 0 && gridfilenames.size() > 0) {
693 0 : error("frequency with which to output grid not specified use GRID_WSTRIDE");
694 : }
695 92 : if(gridfilenames.size()==0 && wgridstride_ > 0) {
696 12 : for(unsigned i=0; i<pf_n; i++) {
697 8 : std::string gridfname = do_pf_ ? "GRID.PF"+std::to_string(i) : "GRID."+getPntrToArgument(i)->getName();
698 8 : gridfilenames.push_back(gridfname);
699 : }
700 : }
701 92 : if(gridfilenames.size() > 0 && hillsfname.size() > 0 && gridfilenames.size() != hillsfname.size()) {
702 0 : error("number of GRID_WFILES arguments does not match number of HILLS files");
703 : }
704 :
705 : // Read grid
706 92 : std::vector<std::string> gridreadfilenames;
707 92 : parseVector("GRID_RFILES",gridreadfilenames);
708 :
709 : // Grid Stuff
710 : // (parseVector checks the number of values: the vectors are empty if the keywords are not given)
711 92 : std::vector<std::string> gmin(pf_n);
712 92 : parseVector("GRID_MIN",gmin);
713 92 : std::vector<std::string> gmax(pf_n);
714 92 : parseVector("GRID_MAX",gmax);
715 92 : std::vector<unsigned> gbin(pf_n);
716 92 : parseVector("GRID_BIN",gbin);
717 92 : std::vector<double> gspacing(pf_n);
718 92 : parseVector("GRID_SPACING",gspacing);
719 92 : const bool gridnobounds=MetaDBiasGrid::setupBins(*this,pfhold,gmin,gmax,gbin,gspacing,adaptive_!=FlexibleBin::none,sigma0_,sigma0min_);
720 92 : if(gbin.size()>0) {
721 40 : grid_=true;
722 : }
723 92 : double gridsizelimit=1.0e8;
724 92 : parse("GRID_SIZE_LIMIT",gridsizelimit);
725 92 : MetaDBiasGrid::checkSizeLimit(*this,gridsizelimit);
726 :
727 92 : bool sparsegrid=false;
728 92 : parseFlag("GRID_SPARSE",sparsegrid);
729 92 : parseFlag("STORE_GRIDS",storeOldGrids_);
730 92 : bool nospline=false;
731 92 : parseFlag("GRID_NOSPLINE",nospline);
732 92 : bool spline=!nospline;
733 92 : if(!grid_&&gridfilenames.size() > 0) {
734 0 : error("To write a grid you need first to define it!");
735 : }
736 92 : if(!grid_&&gridreadfilenames.size() > 0) {
737 0 : error("To read a grid you need first to define it!");
738 : }
739 92 : if(gridreadfilenames.size()>0 && gridreadfilenames.size()!=pf_n) {
740 0 : error("number of GRID_RFILES arguments does not match number of HILLS files");
741 : }
742 :
743 : // Interval keyword
744 92 : std::vector<double> lowI;
745 92 : std::vector<double> uppI;
746 92 : parseVector("INTERVAL_MIN",lowI);
747 184 : parseVector("INTERVAL_MAX",uppI);
748 : // various checks
749 92 : if(lowI.size()!=uppI.size()) {
750 0 : error("both the lower and the upper limits must be provided with INTERVAL_MIN and INTERVAL_MAX");
751 : }
752 92 : if(lowI.size()!=0 && lowI.size()!=pf_n) {
753 0 : error("the number of INTERVAL_MIN and INTERVAL_MAX values should be the same as the number of arguments (PFs)");
754 : }
755 106 : for(unsigned i=0; i<lowI.size(); ++i) {
756 14 : if(uppI[i]<lowI[i]) {
757 0 : error("The Upper limit must be greater than the Lower limit!");
758 : }
759 : // (a value must be given for every family: the values for the periodic ones are not used)
760 14 : if(pfhold[i]->isPeriodic()) {
761 4 : warning("INTERVAL is not used for the periodic variable "+pfhold[i]->getName());
762 4 : continue;
763 : }
764 10 : families_[i].interval=true;
765 10 : families_[i].lower=lowI[i];
766 10 : families_[i].upper=uppI[i];
767 10 : families_[i].kernel.setInterval(lowI[i],uppI[i]);
768 : }
769 :
770 : // parse selector stuff
771 184 : parse("SELECTOR", selector_);
772 92 : if(selector_.length()>0) {
773 3 : do_select_ = true;
774 3 : select_value_ = 0; // default value, used if SELECTOR_ID is not given
775 6 : parse("SELECTOR_ID", select_value_);
776 : }
777 :
778 92 : checkRead();
779 :
780 92 : log.printf(" Gaussian width ");
781 92 : if (adaptive_==FlexibleBin::diffusion) {
782 15 : log.printf(" (Note: The units of sigma are in timesteps) ");
783 : }
784 92 : if (adaptive_==FlexibleBin::geometry) {
785 1 : log.printf(" (Note: The units of sigma are in dist units) ");
786 : }
787 265 : for(unsigned i=0; i<sigma0_.size(); ++i) {
788 173 : log.printf(" %f",sigma0_[i]);
789 : }
790 92 : log.printf("\n Gaussian height %f\n",height0_);
791 92 : log.printf(" Gaussian deposition pace %d\n",stride_);
792 92 : log.printf(" Gaussian files ");
793 283 : for(unsigned i=0; i<hillsfname.size(); ++i) {
794 191 : log.printf("%s ",hillsfname[i].c_str());
795 : }
796 92 : log.printf("\n");
797 92 : if(nobias_) {
798 2 : log.printf(" Bias factor 1: no bias is applied, the hills files and grids contain the estimate of -F\n");
799 : }
800 92 : if(welltemp_) {
801 82 : log.printf(" Well-Tempered Bias Factor %f\n",biasf_);
802 82 : log.printf(" Hills relaxation time (tau) %f\n",tau);
803 82 : log.printf(" KbT %f\n",kbt_);
804 : }
805 :
806 92 : if(do_select_) {
807 3 : log.printf(" Add forces and update bias based on the value of SELECTOR %s\n",selector_.c_str());
808 3 : log.printf(" Id of the SELECTOR for this action %u\n", select_value_);
809 : }
810 :
811 92 : mw_.setup(*this);
812 :
813 283 : for(unsigned f=0; f<pf_n; f++) {
814 191 : if(families_[f].interval) {
815 10 : log.printf(" Upper and Lower limits boundaries for the bias of %s%s are activated at %f - %f\n",
816 20 : do_pf_ ? ("PF"+std::to_string(f)+": ").c_str() : "",pfhold[f]->getName().c_str(),families_[f].lower,families_[f].upper);
817 : }
818 : }
819 92 : if(grid_ && gridnobounds) {
820 8 : log.printf(" Grid spacing");
821 24 : for(unsigned i=0; i<gspacing.size(); ++i) {
822 32 : log.printf(" %s",exactString(gspacing[i]).c_str());
823 : }
824 8 : log.printf(" (periodic variables: grid bin");
825 24 : for(unsigned i=0; i<gbin.size(); ++i) {
826 16 : if(pfhold[i]->isPeriodic()) {
827 8 : log.printf(" %u",gbin[i]);
828 : } else {
829 8 : log.printf(" -");
830 : }
831 : }
832 8 : log.printf(")\n The grids are placed around the first Gaussian\n");
833 84 : } else if(grid_) {
834 32 : log.printf(" Grid min");
835 100 : for(unsigned i=0; i<gmin.size(); ++i) {
836 68 : log.printf(" %s",gmin[i].c_str() );
837 : }
838 32 : log.printf("\n");
839 32 : log.printf(" Grid max");
840 100 : for(unsigned i=0; i<gmax.size(); ++i) {
841 68 : log.printf(" %s",gmax[i].c_str() );
842 : }
843 32 : log.printf("\n");
844 32 : log.printf(" Grid bin");
845 100 : for(unsigned i=0; i<gbin.size(); ++i) {
846 68 : log.printf(" %u",gbin[i]);
847 : }
848 32 : log.printf("\n");
849 : }
850 92 : if(grid_) {
851 40 : if(spline) {
852 40 : log.printf(" Grid uses spline interpolation\n");
853 : }
854 40 : if(sparsegrid) {
855 1 : log.printf(" Grid uses sparse grid\n");
856 : }
857 40 : if(wgridstride_>0) {
858 25 : for(unsigned i=0; i<gridfilenames.size(); ++i) {
859 17 : log.printf(" Grid is written on file %s with stride %d\n",gridfilenames[i].c_str(),wgridstride_);
860 : }
861 : }
862 40 : if(gridreadfilenames.size()>0) {
863 15 : for(unsigned i=0; i<gridreadfilenames.size(); ++i) {
864 10 : log.printf(" Reading bias from grid in file %s \n",gridreadfilenames[i].c_str());
865 : }
866 : }
867 80 : log.printf(" The grids are extended when a Gaussian does not fit in them, up to %s points\n",exactString(gridsizelimit).c_str());
868 : }
869 :
870 : // restart from external grid
871 : bool restartedFromGrid=false;
872 : // hills files of the families, and the numbers of hills of the other walkers in the grids read from file
873 92 : std::vector<MetaDHillsInput> inputs(pf_n);
874 :
875 : // initializing and checking grid
876 92 : if(grid_) {
877 40 : std::string funcl=getLabel() + ".bias";
878 124 : for(unsigned f=0; f<pf_n; ++f) {
879 : // check for mesh and sigma size
880 84 : MetaDBiasGrid::checkSpacing(*this,std::vector<std::string>(1,gmin[f]),std::vector<std::string>(1,gmax[f]),std::vector<unsigned>(1,gbin[f]),
881 239 : adaptive_!=FlexibleBin::none,std::vector<double>(1,sigma0_[adaptive_!=FlexibleBin::none ? 0 : f]),
882 168 : std::vector<double>(1,adaptive_!=FlexibleBin::none ? sigma0min_[f] : -1.0));
883 84 : std::vector<Value*> args(1,pfhold[f]);
884 84 : std::vector<std::string> gmin_t(1,gmin[f]);
885 84 : std::vector<std::string> gmax_t(1,gmax[f]);
886 84 : std::vector<unsigned> gbin_t(1,gbin[f]);
887 84 : families_[f].grid=Tools::make_unique<MetaDBiasGrid>(*this);
888 : MetaDBiasGrid& grid=*families_[f].grid;
889 : // Read grid from file
890 84 : if(gridreadfilenames.size()>0) {
891 20 : inputs[f].gridWalkerHills=readBiasGrid(*this,mw_,grid,gridreadfilenames[f],funcl,args,gmin_t,gmax_t,gbin_t,sparsegrid,spline,gridnobounds);
892 10 : if(getRestart()) {
893 : restartedFromGrid=true;
894 : }
895 : } else {
896 74 : grid.create(funcl,args,gmin_t,gmax_t,gbin_t,sparsegrid,spline,gridnobounds);
897 : }
898 84 : grid.setSizeLimit(gridsizelimit);
899 : // with INTERVAL the bias beyond the limits is constant: the grid contains the interval and at least one point
900 : // beyond each limit, which holds this constant value
901 84 : if(families_[f].interval) {
902 9 : grid.setInterval(families_[f].lower,families_[f].upper,getTime());
903 : }
904 84 : }
905 : }
906 :
907 : // open all files at the beginning and read Gaussians if restarting
908 283 : for(unsigned f=0; f<pf_n; ++f) {
909 191 : Family& family=families_[f];
910 573 : family.hillsfiles=Tools::make_unique<MetaDHillsFiles>(*this,comm,std::vector<Value*>(1,pfhold[f]));
911 : MetaDHillsFiles::Bias bias;
912 191 : bias.add=[this,f](const Gaussian& hill) {
913 1436 : addGaussian(f,hill);
914 191 : };
915 191 : bias.grid=family.grid.get();
916 191 : bias.diffrestart=&family.diffrestart;
917 191 : if(welltemp_) {
918 171 : bias.heightscale=(biasf_-1.0)/biasf_;
919 : }
920 191 : bias.tolerance=0.5*stride_*getTimeStep();
921 : family.hillsfiles->setBias(bias);
922 : // ADAPTIVE=DIFF: hills of the last ~6 averaging times (older ones have weight < exp(-6)), for each argument
923 : // of the family; with MPI walkers the file contains the hills of all walkers
924 382 : family.diffrestart.maxhills=(static_cast<unsigned>(std::ceil(6.0*sigma0_[0]/stride_))+2)*
925 295 : std::max(1u,mw_.mpi_nw)*static_cast<unsigned>(family.args.size());
926 191 : inputs[f].files=family.hillsfiles.get();
927 191 : inputs[f].name=hillsfname[f];
928 : }
929 92 : readHillsAtStart(*this,mw_,inputs,restartedFromGrid,schedule_);
930 :
931 : // open hills files for writing
932 283 : for(unsigned f=0; f<pf_n; ++f) {
933 : auto ofile=Tools::make_unique<OFile>();
934 191 : mw_.openOutput(*this,*ofile,families_[f].hillsfiles->getInputName(mw_.id));
935 191 : if(fmt_.length()>0) {
936 131 : ofile->fmtField(fmt_);
937 : }
938 191 : ofile->addConstantField("multivariate");
939 382 : ofile->addConstantField("kerneltype");
940 191 : if(families_[f].interval) {
941 20 : ofile->addConstantField("lower_int").printField("lower_int",families_[f].lower);
942 20 : ofile->addConstantField("upper_int").printField("upper_int",families_[f].upper);
943 : }
944 191 : ofile->setHeavyFlush();
945 : // output periodicities of variables
946 191 : ofile->setupPrintValue( pfhold[f] ); //assuming cvs in the same family have the same periodicity and boundaries.
947 : families_[f].hillsOfile=std::move(ofile);
948 191 : }
949 :
950 : // Dump grid to files
951 92 : if(wgridstride_ > 0) {
952 25 : for(unsigned f = 0; f < gridfilenames.size(); ++f) {
953 : auto ofile=Tools::make_unique<OFile>();
954 17 : mw_.openOutput(*this,*ofile,gridfilenames[f]);
955 17 : ofile->setHeavyFlush();
956 : families_[f].gridOfile=std::move(ofile);
957 17 : }
958 : }
959 :
960 184 : log<<" Bibliography "<<plumed.cite("Pfaendtner and Bonomi. J. Chem. Theory Comput. 11, 5062 (2015)");
961 : bool interval=false;
962 283 : for(const auto & family : families_) {
963 191 : interval=interval || family.interval;
964 : }
965 92 : if(interval)
966 12 : log<<plumed.cite(
967 12 : "Baftizadeh, Cossio, Pietrucci, and Laio, Curr. Phys. Chem. 2, 79 (2012)");
968 92 : if(mw_.n>1||mw_.mpi)
969 120 : log<<plumed.cite(
970 120 : "Raiteri, Laio, Gervasio, Micheletti, and Parrinello, J. Phys. Chem. B 110, 3533 (2006)");
971 92 : if(adaptive_!=FlexibleBin::none)
972 32 : log<<plumed.cite(
973 32 : "Branduardi, Bussi, and Parrinello, J. Chem. Theory Comput. 8, 2247 (2012)");
974 92 : if (do_pf_) {
975 26 : log<<plumed.cite("Prakash, Fu, Bonomi, and Pfaendtner, J. Chem. Theory Comput. 14, 4985 (2018)");
976 : }
977 92 : log<<"\n";
978 296 : }
979 :
980 2807 : void PBMetaD::initializeDiffusionStatistics() {
981 2807 : const double dt=getTimeStep();
982 8967 : for(auto & family : families_) {
983 6160 : MetaDDiffusionRestart& diffrestart=family.diffrestart;
984 6160 : if(!diffrestart.pending) {
985 6156 : continue;
986 : }
987 4 : diffrestart.pending=false;
988 : // In the hills file of the family, each deposition is a block with a hill for each argument of the family, in
989 : // their order (and with MPI walkers for each walker, in the order of the replicas)
990 4 : const unsigned nf=family.args.size();
991 10 : for(unsigned p=0; p<nf; p++) {
992 6 : const unsigned i=family.args[p];
993 6 : std::vector<double> average;
994 6 : std::vector<double> variance;
995 6 : double lasttime=0.0;
996 : // times are written with a finite precision (FMT): times closer than half a PACE are the same time
997 12 : const unsigned n=diffrestart.compute(0.5*stride_*dt,getTime(),(mw_.mpi ? mw_.mpi_nw : 1)*nf,(mw_.mpi ? mw_.mpi_id : 0)*nf+p,
998 6 : 1.0/sigma0_[0],dt*getStride(),
999 357 : [&](unsigned,double a,double b) {
1000 357 : return difference(i,a,b);
1001 : },
1002 6 : [&](unsigned,double x) {
1003 6 : return bringBackInPbc(i,x);
1004 : },
1005 : average,variance,lasttime);
1006 6 : if(n==0) {
1007 0 : log.printf(" ADAPTIVE=DIFF: no hills available for %s, the statistics for the hills width start from scratch\n",
1008 0 : getPntrToArgument(i)->getName().c_str());
1009 : continue;
1010 : }
1011 6 : flexbin_[i].setStatistics(average,variance,2);
1012 6 : log.printf(" ADAPTIVE=DIFF: the statistics for the hills width of %s are initialized from %u hills of %s (last one at time %f)\n",
1013 6 : getPntrToArgument(i)->getName().c_str(),n,diffrestart.getFile().c_str(),lasttime);
1014 : }
1015 4 : diffrestart.clear();
1016 : }
1017 2807 : }
1018 :
1019 9084 : void PBMetaD::writeGaussian(unsigned f, const Gaussian& hill) {
1020 9084 : double height=hill.height;
1021 9084 : if(welltemp_) {
1022 8594 : height *= biasf_/(biasf_-1.0);
1023 : }
1024 9084 : families_[f].hillsfiles->write(*families_[f].hillsOfile,getTimeStep()*getStep(),hill,height,biasf_,mw_.n>1);
1025 9084 : }
1026 :
1027 10520 : void PBMetaD::addGaussian(unsigned f, const Gaussian& hill) {
1028 10520 : Family& family=families_[f];
1029 10520 : if(family.grid) {
1030 6323 : family.grid->addHill(hill,family.kernel,getTime());
1031 : } else {
1032 4197 : family.hills.push_back(hill);
1033 4197 : family.center.push_back(hill.center[0]);
1034 4197 : family.invsigma.push_back(hill.multivariate ? 0.0 : hill.invsigma[0]);
1035 4197 : family.height.push_back(hill.height);
1036 5721 : family.diagonal.push_back(!hill.multivariate && !family.interval);
1037 : }
1038 10520 : }
1039 :
1040 14962 : void PBMetaD::getBiases(const std::vector<double>& cv, std::vector<double>& bias, std::vector<double>* der) {
1041 14962 : const unsigned nargs=getNumberOfArguments();
1042 14962 : bias.assign(nargs,0.0);
1043 14962 : if(der) {
1044 11944 : der->assign(nargs,0.0);
1045 : }
1046 14962 : std::vector<double> cvi(1);
1047 14962 : std::vector<double> deri(1);
1048 14962 : if(grid_) {
1049 28191 : for(unsigned i=0; i<nargs; ++i) {
1050 19719 : const MetaDBiasGrid& grid=*families_[pfs_[i]].grid;
1051 19719 : cvi[0]=cv[i];
1052 19719 : if(der) {
1053 15580 : bias[i]=grid.getBiasAndDerivatives(cvi,deri);
1054 15580 : (*der)[i]=deri[0];
1055 : } else {
1056 4139 : bias[i]=grid.getBias(cvi);
1057 : }
1058 : }
1059 : return;
1060 : }
1061 : // the hills are split among the MPI processes and the OpenMP threads
1062 6490 : const unsigned stride=comm.Get_size();
1063 6490 : const unsigned rank=comm.Get_rank();
1064 6490 : const unsigned nt=OpenMP::getNumThreads();
1065 22329 : for(unsigned i=0; i<nargs; ++i) {
1066 15839 : const Family& family=families_[pfs_[i]];
1067 15839 : const unsigned nhills=family.hills.size();
1068 15839 : cvi[0]=cv[i];
1069 : double b=0.0;
1070 15839 : double d=0.0;
1071 15839 : if(nt==1 || nhills<2*nt*stride) {
1072 3480 : for(unsigned h=rank; h<nhills; h+=stride) {
1073 2404 : b+=family.evaluate(cvi,h,der ? &d : nullptr);
1074 : }
1075 : } else {
1076 14315 : #pragma omp parallel num_threads(nt)
1077 : {
1078 : double dt=0.0;
1079 : #pragma omp for reduction(+:b) nowait
1080 : for(unsigned h=rank; h<nhills; h+=stride) {
1081 : b+=family.evaluate(cvi,h,der ? &dt : nullptr);
1082 : }
1083 : #pragma omp atomic
1084 : d+=dt;
1085 : }
1086 : }
1087 15839 : bias[i]=b;
1088 15839 : if(der) {
1089 12996 : (*der)[i]=d;
1090 : }
1091 : }
1092 6490 : if(stride>1) {
1093 : // a single reduction for the biases and the derivatives
1094 1064 : std::vector<double> buffer(bias);
1095 1064 : if(der) {
1096 648 : buffer.insert(buffer.end(),der->begin(),der->end());
1097 : }
1098 1064 : comm.Sum(buffer);
1099 : std::copy(buffer.begin(),buffer.begin()+nargs,bias.begin());
1100 1064 : if(der) {
1101 : std::copy(buffer.begin()+nargs,buffer.end(),der->begin());
1102 : }
1103 : }
1104 : }
1105 :
1106 11944 : void PBMetaD::calculate() {
1107 11944 : schedule_.start(getStep());
1108 :
1109 11944 : const unsigned nargs=getNumberOfArguments();
1110 11944 : std::vector<double> cv(nargs);
1111 40520 : for(unsigned i=0; i<nargs; ++i) {
1112 28576 : cv[i]=getArgument(i);
1113 : }
1114 :
1115 : // during a replica exchange attempt the bias is computed on the coordinates of the partner.
1116 : // The MD code does not tell whether the exchange is accepted: with ADAPTIVE=DIFF this is
1117 : // detected at the next update, comparing the CVs with those of the partner (see updateFlexibleBins)
1118 11944 : if(adaptive_==FlexibleBin::diffusion && getExchangeStep()) {
1119 6 : exchange_.attempted=true;
1120 6 : exchange_.partner_cv=cv;
1121 : }
1122 11944 : std::vector<double> bias;
1123 11944 : std::vector<double> deriv;
1124 11944 : getBiases(cv,bias,&deriv);
1125 :
1126 : // the biases are combined as -kT log(1/N sum exp(-V_i/kT)), shifting the exponents by the smallest bias
1127 11944 : const double ncv = static_cast<double>(nargs);
1128 11944 : const double bmin = *std::min_element(bias.begin(),bias.end());
1129 : double ene = 0.;
1130 40520 : for(unsigned i=0; i<nargs; ++i) {
1131 28576 : ene += std::exp((-bias[i]+bmin)/kbt_);
1132 : }
1133 :
1134 : // set Forces - set them to zero if SELECTOR is active
1135 11944 : if(do_select_) {
1136 205 : const auto selector=plumed.passMap.find(selector_);
1137 205 : if(selector==plumed.passMap.end()) {
1138 0 : error("SELECTOR "+selector_+" has not been defined");
1139 : }
1140 205 : current_value_ = static_cast<unsigned>(selector->second);
1141 : }
1142 :
1143 : // with BIASFACTOR=1 no bias is applied
1144 11944 : const bool selected=(!do_select_ || select_value_==current_value_);
1145 40520 : for(unsigned i=0; i<nargs; ++i) {
1146 28576 : const double force = (selected && !nobias_ ? - std::exp((-bias[i]+bmin)/kbt_) / (ene) * deriv[i] : 0.0);
1147 28576 : setOutputForce(i, force);
1148 : }
1149 :
1150 : // set bias
1151 11944 : ene = -kbt_ * (std::log(ene) - std::log(ncv)) + bmin;
1152 23688 : setBias(nobias_ ? 0.0 : ene);
1153 11944 : }
1154 :
1155 11938 : void PBMetaD::beforeUpdate() {
1156 : // called at every step, also outside the update window and when the action is not active
1157 : // (e.g. with STRIDE), in which case the arguments have not been computed
1158 11938 : if(precal_diff_ && isActive()) {
1159 100 : updateFlexibleBins(false);
1160 : }
1161 11938 : }
1162 :
1163 12434 : bool PBMetaD::isHillStep() const {
1164 12434 : return schedule_.isHillStep(getStep(),stride_,getTime(),0.5*stride_*getTimeStep());
1165 : }
1166 :
1167 3353 : void PBMetaD::updateFlexibleBins(bool nowAddAHill) {
1168 3353 : const unsigned nargs=getNumberOfArguments();
1169 3353 : if(adaptive_==FlexibleBin::diffusion) {
1170 2807 : initializeDiffusionStatistics();
1171 : }
1172 3353 : if(adaptive_==FlexibleBin::diffusion) {
1173 2807 : std::vector<double> cv(nargs);
1174 9813 : for(unsigned i=0; i<nargs; ++i) {
1175 7006 : cv[i]=getArgument(i);
1176 : }
1177 2807 : if(exchange_.attempted && exchange_.last_cv.size()==nargs) {
1178 : // the exchange has been accepted if the CVs are now closer to those of the partner than
1179 : // to their own value at the previous update (distances in units of the current fluctuations)
1180 : double dpartner=0.0;
1181 : double down=0.0;
1182 18 : for(unsigned i=0; i<nargs; i++) {
1183 12 : const std::vector<double> var=flexbin_[i].getMatrix();
1184 12 : const double s2=(var.size()>0 && var[0]>0.0 ? var[0] : 1.0);
1185 12 : const double dp=difference(i,exchange_.partner_cv[i],cv[i]);
1186 12 : const double dl=difference(i,exchange_.last_cv[i],cv[i]);
1187 12 : dpartner+=dp*dp/s2;
1188 12 : down+=dl*dl/s2;
1189 : }
1190 6 : if(dpartner<down) {
1191 : // the CVs jumped: restart the running averages, but keep the variances
1192 12 : for(auto & fb : flexbin_) {
1193 8 : fb.restartAverage();
1194 : }
1195 : }
1196 : }
1197 2807 : exchange_.attempted=false;
1198 2807 : exchange_.last_cv=cv;
1199 : }
1200 11997 : for(unsigned i=0; i<nargs; i++) {
1201 8644 : flexbin_[i].update(nowAddAHill,i);
1202 : }
1203 3353 : }
1204 :
1205 11888 : void PBMetaD::update() {
1206 : // adding hills criteria
1207 11888 : bool nowAddAHill=isHillStep();
1208 :
1209 11888 : const unsigned nargs=getNumberOfArguments();
1210 : // if you use adaptive, call the FlexibleBin
1211 : // (with PRECAL_DIFF the statistics have already been updated in beforeUpdate)
1212 : bool multivariate=false;
1213 11888 : if(adaptive_!=FlexibleBin::none) {
1214 3328 : if(!precal_diff_) {
1215 3253 : updateFlexibleBins(nowAddAHill);
1216 : }
1217 : multivariate=true;
1218 : }
1219 :
1220 : // with ADAPTIVE=DIFF the width of a hill requires the fluctuations of at least two CV values
1221 : // (this happens at the first update, e.g. after UPDATE_FROM or when restarting)
1222 11888 : if(nowAddAHill && adaptive_==FlexibleBin::diffusion) {
1223 1654 : for(const auto & fb : flexbin_) {
1224 1195 : if(fb.getNumberOfSamples()<2) {
1225 : nowAddAHill=false;
1226 : }
1227 : }
1228 : }
1229 :
1230 : // with SELECTOR, hills are added only when the selector has the value SELECTOR_ID; with WALKERS_MPI all the walkers
1231 : // take part in the exchange of the hills (also those that do not add hills)
1232 11888 : const bool selected=(!do_select_ || select_value_==current_value_);
1233 11888 : if(nowAddAHill && (selected || mw_.mpi)) {
1234 : // get all biases and heights
1235 3027 : std::vector<double> cv(nargs);
1236 3027 : std::vector<double> thissigma(nargs);
1237 3027 : std::vector<double> height(nargs);
1238 10009 : for(unsigned i=0; i<nargs && selected; ++i) {
1239 : // the width of the hills of each argument
1240 6982 : if(adaptive_!=FlexibleBin::none) {
1241 1517 : thissigma[i]=flexbin_[i].getInverseMatrix(i)[0];
1242 : } else {
1243 5465 : thissigma[i]=sigma0_[pfs_[i]];
1244 : }
1245 6982 : cv[i] = getArgument(i);
1246 : }
1247 3027 : if(selected) {
1248 3018 : std::vector<double> bias;
1249 3018 : getBiases(cv,bias,nullptr);
1250 3018 : if(nobias_) {
1251 : // BIASFACTOR=1: the grids and hills contain -F and the applied bias is zero: the conditional term is 1/N,
1252 : // and the well-tempered factor is exp(-V/kT), the limit for a bias factor going to 1 (as in METAD)
1253 114 : for(unsigned i=0; i<nargs; ++i) {
1254 76 : height[i] = height0_/static_cast<double>(nargs)*std::exp(-bias[i]/kbt_);
1255 : }
1256 : } else {
1257 2980 : const double bmin=*std::min_element(bias.begin(),bias.end());
1258 : // calculate heights and norm
1259 : double norm = 0.0;
1260 9886 : for(unsigned i=0; i<nargs; ++i) {
1261 6906 : double h = std::exp((-bias[i]+bmin)/kbt_);
1262 6906 : norm += h;
1263 6906 : height[i] = h;
1264 : }
1265 : // normalize and apply welltemp correction
1266 9886 : for(unsigned i=0; i<nargs; ++i) {
1267 6906 : height[i] *= height0_ / norm;
1268 6906 : if(welltemp_) {
1269 6492 : height[i] *= std::exp(-bias[i]/(kbt_*(biasf_-1.0)));
1270 : }
1271 : }
1272 : }
1273 : }
1274 :
1275 3027 : std::vector<double> cv_tmp(1);
1276 3027 : std::vector<double> sigma_tmp(1);
1277 : // MPI Multiple walkers: share hills and add them all
1278 3027 : if(mw_.mpi) {
1279 : // the hills of all walkers (packed in a single array: center, sigma, height of each argument, and whether the
1280 : // walker adds hills, so that one collective is enough)
1281 1060 : const std::size_t recordsize=3*nargs+1;
1282 1060 : std::vector<double> record(recordsize,0.0);
1283 3180 : for(unsigned i=0; i<nargs; ++i) {
1284 2120 : record[i]=cv[i];
1285 2120 : record[nargs+i]=thissigma[i];
1286 2120 : record[2*nargs+i]=height[i];
1287 : }
1288 1069 : record[3*nargs]=(selected ? 1.0 : 0.0);
1289 1060 : std::vector<double> all_records(mw_.mpi_nw*recordsize, 0.0);
1290 1060 : if(comm.Get_rank()==0) {
1291 : // Communicate (only root)
1292 540 : multi_sim_comm.Allgather(record,all_records);
1293 : }
1294 : // Share info with group members
1295 1060 : comm.Bcast(all_records,0);
1296 : // now add hills one by one
1297 3180 : for(unsigned j=0; j<mw_.mpi_nw; ++j) {
1298 2120 : const double* rec=&all_records[j*recordsize];
1299 2120 : if(rec[3*nargs]==0.0) {
1300 18 : continue;
1301 : }
1302 6306 : for(unsigned i=0; i<nargs; ++i) {
1303 : // Add CVs of same family together and write to same file
1304 4204 : cv_tmp[0] = rec[i];
1305 4204 : sigma_tmp[0] = rec[nargs+i];
1306 4204 : Gaussian newhill(multivariate, rec[2*nargs+i], cv_tmp, sigma_tmp);
1307 4204 : addGaussian(pfs_[i], newhill);
1308 4204 : writeGaussian(pfs_[i], newhill);
1309 4204 : }
1310 : }
1311 : // just add your own hills
1312 : } else {
1313 6847 : for(unsigned i=0; i<nargs; ++i) {
1314 : // Add CVs of same family together and write to same file
1315 4880 : cv_tmp[0] = cv[i];
1316 4880 : sigma_tmp[0] = thissigma[i];
1317 4880 : Gaussian newhill(multivariate, height[i], cv_tmp, sigma_tmp);
1318 4880 : addGaussian(pfs_[i], newhill);
1319 4880 : writeGaussian(pfs_[i], newhill);
1320 4880 : }
1321 : }
1322 : }
1323 :
1324 : // with multiple walkers the hills must reach the files before the other walkers read them
1325 : // (this should be outside of the if block in case WALKERS_RSTRIDE is not a multiple of PACE)
1326 11888 : if(mw_.isReadStep(getStep())) {
1327 300 : for(auto & family : families_) {
1328 200 : family.hillsOfile->flush();
1329 : }
1330 : }
1331 :
1332 : // write grid files
1333 11888 : if(wgridstride_>0 && (getStep()%wgridstride_==0 || getCPT()) && mw_.writesFiles()) {
1334 53 : for(auto & family : families_) {
1335 36 : if(family.gridOfile) {
1336 : // the hills of the other walkers in the grid (see readHillsAtStart), in a comment after the grid
1337 108 : family.grid->write(*family.gridOfile,storeOldGrids_,mw_.n>1 ? family.hillsfiles->getGridComment(mw_.id) : std::string());
1338 : }
1339 : }
1340 : }
1341 :
1342 : // if multiple walkers and time to read Gaussians
1343 11888 : if(mw_.isReadStep(getStep())) {
1344 100 : std::vector<MetaDHillsFiles*> hillsfiles;
1345 300 : for(auto & family : families_) {
1346 200 : hillsfiles.push_back(family.hillsfiles.get());
1347 : }
1348 100 : readOtherWalkersHills(*this,mw_,hillsfiles);
1349 100 : }
1350 :
1351 11888 : }
1352 :
1353 11950 : bool PBMetaD::checkNeedsGradients()const {
1354 11950 : return adaptive_==FlexibleBin::geometry && isHillStep();
1355 : }
1356 :
1357 : }
1358 : }
|
