Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2018 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "Bias.h"
23 : #include "ActionRegister.h"
24 : #include "core/ActionSet.h"
25 : #include "tools/Grid.h"
26 : #include "core/PlumedMain.h"
27 : #include "core/Atoms.h"
28 : #include "tools/Exception.h"
29 : #include "core/FlexibleBin.h"
30 : #include "tools/Matrix.h"
31 : #include "tools/Random.h"
32 : #include <string>
33 : #include <cstring>
34 : #include "tools/File.h"
35 : #include <iostream>
36 : #include <limits>
37 : #include <ctime>
38 :
39 : #define DP2CUTOFF 6.25
40 :
41 : using namespace std;
42 :
43 :
44 : namespace PLMD {
45 : namespace bias {
46 :
47 : //+PLUMEDOC BIAS METAD
48 : /*
49 : Used to performed MetaDynamics on one or more collective variables.
50 :
51 : In a metadynamics simulations a history dependent bias composed of
52 : intermittently added Gaussian functions is added to the potential \cite metad.
53 :
54 : \f[
55 : V(\vec{s},t) = \sum_{ k \tau < t} W(k \tau)
56 : \exp\left(
57 : -\sum_{i=1}^{d} \frac{(s_i-s_i^{(0)}(k \tau))^2}{2\sigma_i^2}
58 : \right).
59 : \f]
60 :
61 : This potential forces the system away from the kinetic traps in the potential energy surface
62 : and out into the unexplored parts of the energy landscape. Information on the Gaussian
63 : functions from which this potential is composed is output to a file called HILLS, which
64 : is used both the restart the calculation and to reconstruct the free energy as a function of the CVs.
65 : The free energy can be reconstructed from a metadynamics calculation because the final bias is given
66 : by:
67 :
68 : \f[
69 : V(\vec{s}) = -F(\vec(s))
70 : \f]
71 :
72 : During post processing the free energy can be calculated in this way using the \ref sum_hills
73 : utility.
74 :
75 : In the simplest possible implementation of a metadynamics calculation the expense of a metadynamics
76 : calculation increases with the length of the simulation as one has to, at every step, evaluate
77 : the values of a larger and larger number of Gaussians. To avoid this issue you can
78 : store the bias on a grid. This approach is similar to that proposed in \cite babi08jcp but has the
79 : advantage that the grid spacing is independent on the Gaussian width.
80 : Notice that you should
81 : provide either the number of bins for every collective variable (GRID_BIN) or
82 : the desired grid spacing (GRID_SPACING). In case you provide both PLUMED will use
83 : the most conservative choice (highest number of bins) for each dimension.
84 : In case you do not provide any information about bin size (neither GRID_BIN nor GRID_SPACING)
85 : and if Gaussian width is fixed PLUMED will use 1/5 of the Gaussian width as grid spacing.
86 : This default choice should be reasonable for most applications.
87 :
88 : Metadynamics can be restarted either from a HILLS file as well as from a GRID, in this second
89 : case one can first save a GRID using GRID_WFILE (and GRID_WSTRIDE) and at a later stage read
90 : it using GRID_RFILE.
91 :
92 : Another option that is available in plumed is well-tempered metadynamics \cite Barducci:2008. In this
93 : varient of metadynamics the heights of the Gaussian hills are rescaled at each step so the bias is now
94 : given by:
95 :
96 : \f[
97 : V({s},t)= \sum_{t'=0,\tau_G,2\tau_G,\dots}^{t'<t} W e^{-V({s}({q}(t'),t')/\Delta T} \exp\left(
98 : -\sum_{i=1}^{d} \frac{(s_i({q})-s_i({q}(t'))^2}{2\sigma_i^2}
99 : \right),
100 : \f]
101 :
102 : This method ensures that the bias converges more smoothly. It should be noted that, in the case of well-tempered metadynamics, in
103 : the output printed the Gaussian height is re-scaled using the bias factor.
104 : Also notice that with well-tempered metadynamics the HILLS file does not contain the bias,
105 : but the negative of the free-energy estimate. This choice has the advantage that
106 : one can restart a simulation using a different value for the \f$\Delta T\f$. The applied bias will be scaled accordingly.
107 :
108 : Note that you can use here also the flexible gaussian approach \cite Branduardi:2012dl
109 : in which you can adapt the gaussian to the extent of Cartesian space covered by a variable or
110 : to the space in collective variable covered in a given time. In this case the width of the deposited
111 : gaussian potential is denoted by one value only that is a Cartesian space (ADAPTIVE=GEOM) or a time
112 : (ADAPTIVE=DIFF). Note that a specific integration technique for the deposited gaussians
113 : should be used in this case. Check the documentation for utility sum_hills.
114 :
115 : With the keyword INTERVAL one changes the metadynamics algorithm setting the bias force equal to zero
116 : outside boundary \cite baftizadeh2012protein. If, for example, metadynamics is performed on a CV s and one is interested only
117 : to the free energy for s > sw, the history dependent potential is still updated according to the above
118 : equations but the metadynamics force is set to zero for s < sw. Notice that Gaussians are added also
119 : if s < sw, as the tails of these Gaussians influence VG in the relevant region s > sw. In this way, the
120 : force on the system in the region s > sw comes from both metadynamics and the force field, in the region
121 : s < sw only from the latter. This approach allows obtaining a history-dependent bias potential VG that
122 : fluctuates around a stable estimator, equal to the negative of the free energy far enough from the
123 : boundaries. Note that:
124 : - It works only for one-dimensional biases;
125 : - It works both with and without GRID;
126 : - The interval limit sw in a region where the free energy derivative is not large;
127 : - If in the region outside the limit sw the system has a free energy minimum, the INTERVAL keyword should
128 : be used together with a \ref UPPER_WALLS or \ref LOWER_WALLS at sw.
129 :
130 : As a final note, since version 2.0.2 when the system is outside of the selected interval the force
131 : is set to zero and the bias value to the value at the corresponding boundary. This allows acceptances
132 : for replica exchange methods to be computed correctly.
133 :
134 : Multiple walkers \cite multiplewalkers can also be used. See below the examples.
135 :
136 :
137 : The c(t) reweighting factor can also be calculated on the fly using the equations
138 : presented in \cite Tiwary_jp504920s.
139 : The expression used to calculate c(t) follows directly from using Eq. 12 in
140 : Eq. 3 in \cite Tiwary_jp504920s and gives smoother results than equivalent Eqs. 13
141 : and Eqs. 14 in that paper. The c(t) is given by the rct component while the bias
142 : normalized by c(t) is given by the rbias component (rbias=bias-ct) which can be used
143 : to obtain a reweighted histogram.
144 : The calculation of c(t) is enabled by using the keyword REWEIGHTING_NGRID where the grid used for the
145 : calculation is specified. This grid should have a size that is equal or larger than the grid given in GRID_BIN./
146 : By default c(t) is updated every 50 Gaussian hills but this
147 : can be changed by using the REWEIGHTING_NHILLS keyword.
148 : This option can only be employed together with Well-Tempered Metadynamics and requires that
149 : a grid is used.
150 :
151 : Additional material and examples can be also found in the tutorials:
152 :
153 : - \ref belfast-6
154 : - \ref belfast-7
155 : - \ref belfast-8
156 :
157 : Notice that at variance with PLUMED 1.3 it is now straightforward to apply concurrent metadynamics
158 : as done e.g. in Ref. \cite gil2015enhanced . This indeed can be obtained by using the METAD
159 : action multiple times in the same input file.
160 :
161 : \par Examples
162 :
163 : The following input is for a standard metadynamics calculation using as
164 : collective variables the distance between atoms 3 and 5
165 : and the distance between atoms 2 and 4. The value of the CVs and
166 : the metadynamics bias potential are written to the COLVAR file every 100 steps.
167 : \verbatim
168 : DISTANCE ATOMS=3,5 LABEL=d1
169 : DISTANCE ATOMS=2,4 LABEL=d2
170 : METAD ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3 PACE=500 LABEL=restraint
171 : PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
172 : \endverbatim
173 : (See also \ref DISTANCE \ref PRINT).
174 :
175 : \par
176 : If you use adaptive Gaussians, with diffusion scheme where you use
177 : a Gaussian that should cover the space of 20 timesteps in collective variables.
178 : Note that in this case the histogram correction is needed when summing up hills.
179 : \verbatim
180 : DISTANCE ATOMS=3,5 LABEL=d1
181 : DISTANCE ATOMS=2,4 LABEL=d2
182 : METAD ARG=d1,d2 SIGMA=20 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=DIFF
183 : PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
184 : \endverbatim
185 :
186 : \par
187 : If you use adaptive Gaussians, with geometrical scheme where you use
188 : a Gaussian that should cover the space of 0.05 nm in Cartesian space.
189 : Note that in this case the histogram correction is needed when summing up hills.
190 : \verbatim
191 : DISTANCE ATOMS=3,5 LABEL=d1
192 : DISTANCE ATOMS=2,4 LABEL=d2
193 : METAD ARG=d1,d2 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=GEOM
194 : PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
195 : \endverbatim
196 :
197 : \par
198 : When using adaptive Gaussians you might want to limit how the hills width can change.
199 : You can use SIGMA_MIN and SIGMA_MAX keywords.
200 : The sigmas should specified in terms of CV so you should use the CV units.
201 : Note that if you use a negative number, this means that the limit is not set.
202 : Note also that in this case the histogram correction is needed when summing up hills.
203 : \verbatim
204 : DISTANCE ATOMS=3,5 LABEL=d1
205 : DISTANCE ATOMS=2,4 LABEL=d2
206 : METAD ...
207 : ARG=d1,d2 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=GEOM
208 : SIGMA_MIN=0.2,0.1 SIGMA_MAX=0.5,1.0
209 : ... METAD
210 : PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
211 : \endverbatim
212 :
213 : \par
214 : Multiple walkers can be also use as in \cite multiplewalkers
215 : These are enabled by setting the number of walker used, the id of the
216 : current walker which interprets the input file, the directory where the
217 : hills containing files resides, and the frequency to read the other walkers.
218 : Here is an example
219 : \verbatim
220 : DISTANCE ATOMS=3,5 LABEL=d1
221 : METAD ...
222 : ARG=d1 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint
223 : WALKERS_N=10
224 : WALKERS_ID=3
225 : WALKERS_DIR=../
226 : WALKERS_RSTRIDE=100
227 : ... METAD
228 : \endverbatim
229 : where WALKERS_N is the total number of walkers, WALKERS_ID is the
230 : id of the present walker (starting from 0 ) and the WALKERS_DIR is the directory
231 : where all the walkers are located. WALKERS_RSTRIDE is the number of step between
232 : one update and the other. Since version 2.2.5, hills files are automatically
233 : flushed every WALKERS_RSTRIDE steps.
234 :
235 : \par
236 : The c(t) reweighting factor can be calculated on the fly using the equations
237 : presented in \cite Tiwary_jp504920s as described above.
238 : This is enabled by using the keyword REWEIGHTING_NGRID where the grid used for
239 : the calculation is set. The number of grid points given in REWEIGHTING_NGRID
240 : should be equal or larger than the number of grid points given in GRID_BIN.
241 : \verbatim
242 : METAD ...
243 : LABEL=metad
244 : ARG=phi,psi SIGMA=0.20,0.20 HEIGHT=1.20 BIASFACTOR=5 TEMP=300.0 PACE=500
245 : GRID_MIN=-pi,-pi GRID_MAX=pi,pi GRID_BIN=150,150
246 : REWEIGHTING_NGRID=150,150
247 : REWEIGHTING_NHILLS=20
248 : ... METAD
249 : \endverbatim
250 : Here we have asked that the calculation is performed every 20 hills by using
251 : REWEIGHTING_NHILLS keyword. If this keyword is not given the calculation will
252 : by default be performed every 50 hills. The c(t) reweighting factor will be given
253 : in the rct component while the instantaneous value of the bias potential
254 : normalized using the c(t) reweighting factor is given in the rbias component
255 : [rbias=bias-c(t)] which can be used to obtain a reweighted histogram or
256 : free energy surface using the \ref HISTOGRAM analysis.
257 :
258 : \par
259 : The kinetics of the transitions between basins can also be analysed on the fly as
260 : in \cite PRL230602. The flag ACCELERATION turn on accumulation of the acceleration
261 : factor that can then be used to determine the rate. This method can be used together
262 : with \ref COMMITTOR analysis to stop the simulation when the system get to the target basin.
263 : It must be used together with Well-Tempered Metadynamics.
264 :
265 : \par
266 : You can also provide a target distribution using the keyword TARGET
267 : \cite white2015designing
268 : \cite marinelli2015ensemble
269 : \cite gil2016empirical
270 : The TARGET should be a grid containing a free-energy (i.e. the -kbT*log of the desired target distribution).
271 : Gaussians will then be scaled by a factor
272 : \f[
273 : e^{\beta(\tilde{F}(s)-\tilde{F}_{max})}
274 : \f]
275 : Here \f$\tilde{F}(s)\f$ is the free energy defined on the grid and \f$\tilde{F}_{max}\f$ its maximum value.
276 : Notice that we here used the maximum value as in ref \cite gil2016empirical
277 : This choice allows to avoid exceedingly large Gaussians to be added. However,
278 : it could make the Gaussian too small. You should always choose carefully the HEIGHT parameter
279 : in this case.
280 : The grid file should be similar to other PLUMED grid files in that it should contain
281 : both the target free-energy and its derivatives.
282 :
283 : Notice that if you wish your simulation to converge to the target free energy you should use
284 : the DAMPFACTOR command to provide a global tempering \cite dama2014well
285 : Alternatively, if you use a BIASFACTOR yout simulation will converge to a free
286 : energy that is a linear combination of the target free energy and of the intrinsic free energy
287 : determined by the original force field.
288 :
289 : \verbatim
290 : DISTANCE ATOMS=3,5 LABEL=d1
291 : METAD ...
292 : LABEL=t1
293 : ARG=d1 SIGMA=0.05 TAU=200 DAMPFACTOR=100 PACE=250
294 : GRID_MIN=0 GRID_MAX=2 GRID_BIN=200
295 : TARGET=dist.dat
296 : ... METAD
297 :
298 : PRINT ARG=d1,t1.bias STRIDE=100 FILE=COLVAR
299 : \endverbatim
300 :
301 : The header in the file dist.dat for this calculation would read:
302 :
303 : \verbatim
304 : #! FIELDS d1 t1.target der_d1
305 : #! SET min_d1 0
306 : #! SET max_d1 2
307 : #! SET nbins_d1 200
308 : #! SET periodic_d1 false
309 : \endverbatim
310 :
311 : */
312 : //+ENDPLUMEDOC
313 :
314 : class MetaD : public Bias {
315 :
316 : private:
317 29960 : struct Gaussian {
318 : vector<double> center;
319 : vector<double> sigma;
320 : double height;
321 : bool multivariate; // this is required to discriminate the one dimensional case
322 : vector<double> invsigma;
323 5282 : Gaussian(const vector<double> & center,const vector<double> & sigma,double height, bool multivariate ):
324 5282 : center(center),sigma(sigma),height(height),multivariate(multivariate),invsigma(sigma) {
325 : // to avoid troubles from zero element in flexible hills
326 5282 : for(unsigned i=0; i<invsigma.size(); ++i) abs(invsigma[i])>1.e-20?invsigma[i]=1.0/invsigma[i]:0.;
327 5282 : }
328 : };
329 : vector<double> sigma0_;
330 : vector<double> sigma0min_;
331 : vector<double> sigma0max_;
332 : vector<Gaussian> hills_;
333 : OFile hillsOfile_;
334 : OFile gridfile_;
335 : Grid* BiasGrid_;
336 : bool storeOldGrids_;
337 : int wgridstride_;
338 : bool grid_;
339 : double height0_;
340 : double biasf_;
341 : double dampfactor_;
342 : std::string targetfilename_;
343 : Grid* TargetGrid_;
344 : double kbt_;
345 : int stride_;
346 : bool welltemp_;
347 : double* dp_;
348 : int adaptive_;
349 : FlexibleBin *flexbin;
350 : int mw_n_;
351 : string mw_dir_;
352 : int mw_id_;
353 : int mw_rstride_;
354 : bool walkers_mpi;
355 : unsigned mpi_nw_;
356 : bool acceleration;
357 : double acc;
358 : vector<IFile*> ifiles;
359 : vector<string> ifilesnames;
360 : double uppI_;
361 : double lowI_;
362 : bool doInt_;
363 : bool isFirstStep;
364 : double reweight_factor;
365 : vector<unsigned> rewf_grid_;
366 : unsigned rewf_ustride_;
367 : double work_;
368 : long int last_step_warn_grid;
369 :
370 : void readGaussians(IFile*);
371 : bool readChunkOfGaussians(IFile *ifile, unsigned n);
372 : void writeGaussian(const Gaussian&,OFile&);
373 : void addGaussian(const Gaussian&);
374 : double getHeight(const vector<double>&);
375 : double getBiasAndDerivatives(const vector<double>&,double* der=NULL);
376 : double evaluateGaussian(const vector<double>&, const Gaussian&,double* der=NULL);
377 : void finiteDifferenceGaussian(const vector<double>&, const Gaussian&);
378 : double getGaussianNormalization( const Gaussian& );
379 : vector<unsigned> getGaussianSupport(const Gaussian&);
380 : bool scanOneHill(IFile *ifile, vector<Value> &v, vector<double> ¢er, vector<double> &sigma, double &height, bool &multivariate);
381 : void computeReweightingFactor();
382 : string fmt;
383 :
384 : public:
385 : explicit MetaD(const ActionOptions&);
386 : ~MetaD();
387 : void calculate();
388 : void update();
389 : static void registerKeywords(Keywords& keys);
390 5743 : bool checkNeedsGradients()const {if(adaptive_==FlexibleBin::geometry) {return true;} else {return false;}}
391 : };
392 :
393 2617 : PLUMED_REGISTER_ACTION(MetaD,"METAD")
394 :
395 95 : void MetaD::registerKeywords(Keywords& keys) {
396 95 : Bias::registerKeywords(keys);
397 : keys.addOutputComponent("rbias","REWEIGHTING_NGRID","the instantaneous value of the bias normalized using the \\f$c(t)\\f$ reweighting factor [rbias=bias-c(t)]."
398 95 : "This component can be used to obtain a reweighted histogram.");
399 95 : keys.addOutputComponent("rct","REWEIGHTING_NGRID","the reweighting factor \\f$c(t)\\f$.");
400 95 : keys.addOutputComponent("work","default","accumulator for work");
401 95 : keys.addOutputComponent("acc","ACCELERATION","the metadynamics acceleration factor");
402 95 : keys.use("ARG");
403 95 : keys.add("compulsory","SIGMA","the widths of the Gaussian hills");
404 95 : keys.add("compulsory","PACE","the frequency for hill addition");
405 95 : keys.add("compulsory","FILE","HILLS","a file in which the list of added hills is stored");
406 95 : keys.add("optional","HEIGHT","the heights of the Gaussian hills. Compulsory unless TAU and either BIASFACTOR or DAMPFACTOR are given");
407 95 : keys.add("optional","FMT","specify format for HILLS files (useful for decrease the number of digits in regtests)");
408 95 : keys.add("optional","BIASFACTOR","use well tempered metadynamics and use this biasfactor. Please note you must also specify temp");
409 95 : keys.add("optional","DAMPFACTOR","damp hills with exp(-max(V)/(kbT*DAMPFACTOR)");
410 95 : keys.add("optional","TARGET","target to a predefined distribution");
411 95 : keys.add("optional","TEMP","the system temperature - this is only needed if you are doing well-tempered metadynamics");
412 95 : keys.add("optional","TAU","in well tempered metadynamics, sets height to (kb*DeltaT*pace*timestep)/tau");
413 95 : keys.add("optional","GRID_MIN","the lower bounds for the grid");
414 95 : keys.add("optional","GRID_MAX","the upper bounds for the grid");
415 95 : keys.add("optional","GRID_BIN","the number of bins for the grid");
416 95 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
417 : keys.add("optional","REWEIGHTING_NGRID","calculate the c(t) reweighting factor and use that to obtain the normalized bias [rbias=bias-c(t)]."
418 : "Here you should specify the number of grid points required in each dimension."
419 : "The number of grid points should be equal or larger to the number of grid points given in GRID_BIN."
420 95 : "This method is not compatible with metadynamics not on a grid.");
421 : keys.add("optional","REWEIGHTING_NHILLS","how many Gaussian hills should be deposited between calculating the c(t) reweighting factor."
422 95 : "The default is to do this every 50 hills.");
423 95 : keys.addFlag("GRID_SPARSE",false,"use a sparse grid to store hills");
424 95 : keys.addFlag("GRID_NOSPLINE",false,"don't use spline interpolation with grids");
425 95 : keys.add("optional","GRID_WSTRIDE","write the grid to a file every N steps");
426 95 : keys.add("optional","GRID_WFILE","the file on which to write the grid");
427 95 : keys.add("optional","GRID_RFILE","a grid file from which the bias should be read at the initial step of the simulation");
428 95 : keys.addFlag("STORE_GRIDS",false,"store all the grid files the calculation generates. They will be deleted if this keyword is not present");
429 95 : 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");
430 95 : keys.add("optional","WALKERS_ID", "walker id");
431 95 : keys.add("optional","WALKERS_N", "number of walkers");
432 95 : keys.add("optional","WALKERS_DIR", "shared directory with the hills files from all the walkers");
433 95 : keys.add("optional","WALKERS_RSTRIDE","stride for reading hills files");
434 95 : keys.add("optional","INTERVAL","monodimensional lower and upper limits, outside the limits the system will not feel the biasing force.");
435 95 : keys.add("optional","SIGMA_MAX","the upper bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
436 95 : keys.add("optional","SIGMA_MIN","the lower bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
437 95 : keys.addFlag("WALKERS_MPI",false,"Switch on MPI version of multiple walkers - not compatible with other WALKERS_* options");
438 95 : keys.addFlag("ACCELERATION",false,"Set to TRUE if you want to compute the metadynamics acceleration factor.");
439 95 : keys.use("RESTART");
440 95 : keys.use("UPDATE_FROM");
441 95 : keys.use("UPDATE_UNTIL");
442 95 : }
443 :
444 376 : MetaD::~MetaD() {
445 94 : if(flexbin) delete flexbin;
446 94 : if(BiasGrid_) delete BiasGrid_;
447 94 : if(TargetGrid_) delete TargetGrid_;
448 94 : hillsOfile_.close();
449 94 : if(wgridstride_>0) gridfile_.close();
450 94 : delete [] dp_;
451 : // close files
452 200 : for(int i=0; i<mw_n_; ++i) {
453 106 : if(ifiles[i]->isOpen()) ifiles[i]->close();
454 106 : delete ifiles[i];
455 : }
456 282 : }
457 :
458 94 : MetaD::MetaD(const ActionOptions& ao):
459 : PLUMED_BIAS_INIT(ao),
460 : // Grid stuff initialization
461 : BiasGrid_(NULL), wgridstride_(0), grid_(false),
462 : // Metadynamics basic parameters
463 94 : height0_(std::numeric_limits<double>::max()), biasf_(1.0), dampfactor_(0.0), TargetGrid_(NULL),
464 : kbt_(0.0),
465 : stride_(0), welltemp_(false),
466 : // Other stuff
467 : dp_(NULL), adaptive_(FlexibleBin::none),
468 : flexbin(NULL),
469 : // Multiple walkers initialization
470 : mw_n_(1), mw_dir_("./"), mw_id_(0), mw_rstride_(1),
471 : walkers_mpi(false), mpi_nw_(0),
472 : acceleration(false), acc(0.0),
473 : // Interval initialization
474 : uppI_(-1), lowI_(-1), doInt_(false),
475 : isFirstStep(true),
476 : reweight_factor(0.0),
477 : rewf_ustride_(1),
478 : work_(0),
479 188 : last_step_warn_grid(0)
480 : {
481 : // parse the flexible hills
482 94 : string adaptiveoption;
483 94 : adaptiveoption="NONE";
484 94 : parse("ADAPTIVE",adaptiveoption);
485 94 : if(adaptiveoption=="GEOM") {
486 19 : log.printf(" Uses Geometry-based hills width: sigma must be in distance units and only one sigma is needed\n");
487 19 : adaptive_=FlexibleBin::geometry;
488 75 : } else if(adaptiveoption=="DIFF") {
489 2 : log.printf(" Uses Diffusion-based hills width: sigma must be in timesteps and only one sigma is needed\n");
490 2 : adaptive_=FlexibleBin::diffusion;
491 73 : } else if(adaptiveoption=="NONE") {
492 73 : adaptive_=FlexibleBin::none;
493 : } else {
494 0 : error("I do not know this type of adaptive scheme");
495 : }
496 :
497 94 : parse("FMT",fmt);
498 :
499 : // parse the sigma
500 94 : parseVector("SIGMA",sigma0_);
501 94 : if(adaptive_==FlexibleBin::none) {
502 : // if you use normal sigma you need one sigma per argument
503 73 : if( sigma0_.size()!=getNumberOfArguments() ) error("number of arguments does not match number of SIGMA parameters");
504 : } else {
505 : // if you use flexible hills you need one sigma
506 21 : if(sigma0_.size()!=1) {
507 0 : error("If you choose ADAPTIVE you need only one sigma according to your choice of type (GEOM/DIFF)");
508 : }
509 : // if adaptive then the number must be an integer
510 21 : if(adaptive_==FlexibleBin::diffusion) {
511 2 : if(int(sigma0_[0])-sigma0_[0]>1.e-9 || int(sigma0_[0])-sigma0_[0] <-1.e-9 || int(sigma0_[0])<1 ) {
512 0 : error("In case of adaptive hills with diffusion, the sigma must be an integer which is the number of timesteps\n");
513 : }
514 : }
515 : // here evtl parse the sigma min and max values
516 21 : parseVector("SIGMA_MIN",sigma0min_);
517 21 : if(sigma0min_.size()>0 && sigma0min_.size()!=getNumberOfArguments()) {
518 0 : error("the number of SIGMA_MIN values be the same of the number of the arguments");
519 21 : } else if(sigma0min_.size()==0) {
520 21 : sigma0min_.resize(getNumberOfArguments());
521 21 : for(unsigned i=0; i<getNumberOfArguments(); i++) {sigma0min_[i]=-1.;}
522 : }
523 :
524 21 : parseVector("SIGMA_MAX",sigma0max_);
525 21 : if(sigma0max_.size()>0 && sigma0max_.size()!=getNumberOfArguments()) {
526 0 : error("the number of SIGMA_MAX values be the same of the number of the arguments");
527 21 : } else if(sigma0max_.size()==0) {
528 21 : sigma0max_.resize(getNumberOfArguments());
529 21 : for(unsigned i=0; i<getNumberOfArguments(); i++) {sigma0max_[i]=-1.;}
530 : }
531 :
532 21 : flexbin=new FlexibleBin(adaptive_,this,sigma0_[0],sigma0min_,sigma0max_);
533 : }
534 : // note: HEIGHT is not compulsory, since one could use the TAU keyword, see below
535 94 : parse("HEIGHT",height0_);
536 94 : parse("PACE",stride_);
537 94 : if(stride_<=0 ) error("frequency for hill addition is nonsensical");
538 188 : string hillsfname="HILLS";
539 94 : parse("FILE",hillsfname);
540 94 : parse("BIASFACTOR",biasf_);
541 94 : if( biasf_<1.0 ) error("well tempered bias factor is nonsensical");
542 94 : parse("DAMPFACTOR",dampfactor_);
543 94 : double temp=0.0;
544 94 : parse("TEMP",temp);
545 94 : if(temp>0.0) kbt_=plumed.getAtoms().getKBoltzmann()*temp;
546 80 : else kbt_=plumed.getAtoms().getKbT();
547 94 : if(biasf_>1.0) {
548 12 : if(kbt_==0.0) error("Unless the MD engine passes the temperature to plumed, with well-tempered metad you must specify it using TEMP");
549 12 : welltemp_=true;
550 : }
551 94 : if(dampfactor_>0.0) {
552 1 : if(kbt_==0.0) error("Unless the MD engine passes the temperature to plumed, with damped metad you must specify it using TEMP");
553 : }
554 94 : parse("TARGET",targetfilename_);
555 94 : if(targetfilename_.length()>0 && kbt_==0.0) error("with TARGET temperature must be specified");
556 94 : double tau=0.0;
557 94 : parse("TAU",tau);
558 94 : if(tau==0.0) {
559 92 : if(height0_==std::numeric_limits<double>::max()) error("At least one between HEIGHT and TAU should be specified");
560 : // if tau is not set, we compute it here from the other input parameters
561 92 : if(welltemp_) tau=(kbt_*(biasf_-1.0))/height0_*getTimeStep()*stride_;
562 81 : else if(dampfactor_>0.0) tau=(kbt_*dampfactor_)/height0_*getTimeStep()*stride_;
563 : } else {
564 2 : if(height0_!=std::numeric_limits<double>::max()) error("At most one between HEIGHT and TAU should be specified");
565 2 : if(welltemp_) height0_=(kbt_*(biasf_-1.0))/tau*getTimeStep()*stride_;
566 1 : else if(dampfactor_>0.0) height0_=(kbt_*dampfactor_)/tau*getTimeStep()*stride_;
567 0 : else error("TAU only makes sense in well-tempered or damped metadynamics");
568 : }
569 :
570 : // Grid Stuff
571 188 : vector<std::string> gmin(getNumberOfArguments());
572 94 : parseVector("GRID_MIN",gmin);
573 94 : if(gmin.size()!=getNumberOfArguments() && gmin.size()!=0) error("not enough values for GRID_MIN");
574 188 : vector<std::string> gmax(getNumberOfArguments());
575 94 : parseVector("GRID_MAX",gmax);
576 94 : if(gmax.size()!=getNumberOfArguments() && gmax.size()!=0) error("not enough values for GRID_MAX");
577 188 : vector<unsigned> gbin(getNumberOfArguments());
578 188 : vector<double> gspacing;
579 94 : parseVector("GRID_BIN",gbin);
580 94 : if(gbin.size()!=getNumberOfArguments() && gbin.size()!=0) error("not enough values for GRID_BIN");
581 94 : parseVector("GRID_SPACING",gspacing);
582 94 : if(gspacing.size()!=getNumberOfArguments() && gspacing.size()!=0) error("not enough values for GRID_SPACING");
583 94 : if(gmin.size()!=gmax.size()) error("GRID_MAX and GRID_MIN should be either present or absent");
584 94 : if(gspacing.size()!=0 && gmin.size()==0) error("If GRID_SPACING is present also GRID_MIN should be present");
585 94 : if(gbin.size()!=0 && gmin.size()==0) error("If GRID_SPACING is present also GRID_MIN should be present");
586 94 : if(gmin.size()!=0) {
587 34 : if(gbin.size()==0 && gspacing.size()==0) {
588 1 : if(adaptive_==FlexibleBin::none) {
589 1 : log<<" Binsize not specified, 1/5 of sigma will be be used\n";
590 1 : plumed_assert(sigma0_.size()==getNumberOfArguments());
591 1 : gspacing.resize(getNumberOfArguments());
592 1 : for(unsigned i=0; i<gspacing.size(); i++) gspacing[i]=0.2*sigma0_[i];
593 : } else {
594 : // with adaptive hills and grid a sigma min must be specified
595 0 : for(unsigned i=0; i<sigma0min_.size(); i++) if(sigma0min_[i]<=0) error("When using Adaptive Gaussians on a grid SIGMA_MIN must be specified");
596 0 : log<<" Binsize not specified, 1/5 of sigma_min will be be used\n";
597 0 : gspacing.resize(getNumberOfArguments());
598 0 : for(unsigned i=0; i<gspacing.size(); i++) gspacing[i]=0.2*sigma0min_[i];
599 : }
600 33 : } else if(gspacing.size()!=0 && gbin.size()==0) {
601 1 : log<<" The number of bins will be estimated from GRID_SPACING\n";
602 32 : } else if(gspacing.size()!=0 && gbin.size()!=0) {
603 1 : log<<" You specified both GRID_BIN and GRID_SPACING\n";
604 1 : log<<" The more conservative (highest) number of bins will be used for each variable\n";
605 : }
606 34 : if(gbin.size()==0) gbin.assign(getNumberOfArguments(),1);
607 40 : if(gspacing.size()!=0) for(unsigned i=0; i<getNumberOfArguments(); i++) {
608 : double a,b;
609 6 : Tools::convert(gmin[i],a);
610 6 : Tools::convert(gmax[i],b);
611 6 : unsigned n=((b-a)/gspacing[i])+1;
612 6 : if(gbin[i]<n) gbin[i]=n;
613 : }
614 : }
615 94 : bool sparsegrid=false;
616 94 : parseFlag("GRID_SPARSE",sparsegrid);
617 94 : bool nospline=false;
618 94 : parseFlag("GRID_NOSPLINE",nospline);
619 94 : bool spline=!nospline;
620 94 : if(gbin.size()>0) {grid_=true;}
621 94 : parse("GRID_WSTRIDE",wgridstride_);
622 188 : string gridfilename_;
623 94 : parse("GRID_WFILE",gridfilename_);
624 94 : parseFlag("STORE_GRIDS",storeOldGrids_);
625 94 : if(grid_ && gridfilename_.length()>0) {
626 16 : if(wgridstride_==0 ) error("frequency with which to output grid not specified use GRID_WSTRIDE");
627 : }
628 :
629 94 : if(grid_ && wgridstride_>0) {
630 16 : if(gridfilename_.length()==0) error("grid filename not specified use GRID_WFILE");
631 : }
632 188 : string gridreadfilename_;
633 94 : parse("GRID_RFILE",gridreadfilename_);
634 :
635 94 : if(!grid_&&gridfilename_.length()> 0) error("To write a grid you need first to define it!");
636 94 : if(!grid_&&gridreadfilename_.length()>0) error("To read a grid you need first to define it!");
637 :
638 94 : if(grid_) {
639 34 : parseVector("REWEIGHTING_NGRID",rewf_grid_);
640 34 : if(rewf_grid_.size()>0 && rewf_grid_.size()!=getNumberOfArguments()) {
641 0 : error("size mismatch for REWEIGHTING_NGRID keyword");
642 34 : } else if(rewf_grid_.size()==getNumberOfArguments()) {
643 15 : for(unsigned j=0; j<getNumberOfArguments(); ++j) {
644 10 : if( !getPntrToArgument(j)->isPeriodic() ) rewf_grid_[j] += 1;
645 : }
646 : }
647 34 : if(adaptive_==FlexibleBin::diffusion || adaptive_==FlexibleBin::geometry) warning("reweighting has not been proven to work with adaptive Gaussians");
648 34 : rewf_ustride_=50; parse("REWEIGHTING_NHILLS",rewf_ustride_);
649 : }
650 94 : if(dampfactor_>0.0) {
651 1 : if(!grid_) error("With DAMPFACTOR you should use grids");
652 : }
653 :
654 : // Multiple walkers
655 94 : parse("WALKERS_N",mw_n_);
656 94 : parse("WALKERS_ID",mw_id_);
657 94 : if(mw_n_<=mw_id_) error("walker ID should be a numerical value less than the total number of walkers");
658 94 : parse("WALKERS_DIR",mw_dir_);
659 94 : parse("WALKERS_RSTRIDE",mw_rstride_);
660 :
661 : // MPI version
662 94 : parseFlag("WALKERS_MPI",walkers_mpi);
663 :
664 : // Inteval keyword
665 188 : vector<double> tmpI(2);
666 94 : parseVector("INTERVAL",tmpI);
667 94 : if(tmpI.size()!=2&&tmpI.size()!=0) error("both a lower and an upper limits must be provided with INTERVAL");
668 94 : else if(tmpI.size()==2) {
669 2 : lowI_=tmpI.at(0);
670 2 : uppI_=tmpI.at(1);
671 2 : if(getNumberOfArguments()!=1) error("INTERVAL limits correction works only for monodimensional metadynamics!");
672 2 : if(uppI_<lowI_) error("The Upper limit must be greater than the Lower limit!");
673 2 : if(getPntrToArgument(0)->isPeriodic()) error("INTERVAL cannot be used with periodic variables!");
674 2 : doInt_=true;
675 : }
676 :
677 94 : acceleration=false;
678 94 : parseFlag("ACCELERATION",acceleration);
679 :
680 94 : checkRead();
681 :
682 94 : log.printf(" Gaussian width ");
683 94 : if (adaptive_==FlexibleBin::diffusion)log.printf(" (Note: The units of sigma are in timesteps) ");
684 94 : if (adaptive_==FlexibleBin::geometry)log.printf(" (Note: The units of sigma are in dist units) ");
685 94 : for(unsigned i=0; i<sigma0_.size(); ++i) log.printf(" %f",sigma0_[i]);
686 94 : log.printf(" Gaussian height %f\n",height0_);
687 94 : log.printf(" Gaussian deposition pace %d\n",stride_);
688 94 : log.printf(" Gaussian file %s\n",hillsfname.c_str());
689 94 : if(welltemp_) {
690 12 : log.printf(" Well-Tempered Bias Factor %f\n",biasf_);
691 12 : log.printf(" Hills relaxation time (tau) %f\n",tau);
692 12 : log.printf(" KbT %f\n",kbt_);
693 : }
694 94 : if(doInt_) log.printf(" Upper and Lower limits boundaries for the bias are activated at %f - %f\n", lowI_, uppI_);
695 94 : if(grid_) {
696 34 : log.printf(" Grid min");
697 34 : for(unsigned i=0; i<gmin.size(); ++i) log.printf(" %s",gmin[i].c_str() );
698 34 : log.printf("\n");
699 34 : log.printf(" Grid max");
700 34 : for(unsigned i=0; i<gmax.size(); ++i) log.printf(" %s",gmax[i].c_str() );
701 34 : log.printf("\n");
702 34 : log.printf(" Grid bin");
703 34 : for(unsigned i=0; i<gbin.size(); ++i) log.printf(" %u",gbin[i]);
704 34 : log.printf("\n");
705 34 : if(spline) {log.printf(" Grid uses spline interpolation\n");}
706 34 : if(sparsegrid) {log.printf(" Grid uses sparse grid\n");}
707 34 : if(wgridstride_>0) {log.printf(" Grid is written on file %s with stride %d\n",gridfilename_.c_str(),wgridstride_);}
708 : }
709 :
710 94 : if(mw_n_>1) {
711 6 : if(walkers_mpi) error("MPI version of multiple walkers is not compatible with filesystem version of multiple walkers");
712 6 : log.printf(" %d multiple walkers active\n",mw_n_);
713 6 : log.printf(" walker id %d\n",mw_id_);
714 6 : log.printf(" reading stride %d\n",mw_rstride_);
715 6 : log.printf(" directory with hills files %s\n",mw_dir_.c_str());
716 : } else {
717 88 : if(walkers_mpi) {
718 18 : log.printf(" Multiple walkers active using MPI communnication\n");
719 18 : if(comm.Get_rank()==0) {
720 : // Only root of group can communicate with other walkers
721 12 : mpi_nw_=multi_sim_comm.Get_size();
722 : }
723 : // Communicate to the other members of the same group
724 : // info abount number of walkers and walker index
725 18 : comm.Bcast(mpi_nw_,0);
726 : }
727 : }
728 :
729 94 : if( rewf_grid_.size()>0 ) {
730 5 : addComponent("rbias"); componentIsNotPeriodic("rbias");
731 5 : addComponent("rct"); componentIsNotPeriodic("rct");
732 5 : log.printf(" the c(t) reweighting factor will be calculated every %u hills\n",rewf_ustride_);
733 5 : getPntrToComponent("rct")->set(reweight_factor);
734 : }
735 94 : addComponent("work"); componentIsNotPeriodic("work");
736 :
737 94 : if(acceleration) {
738 0 : if(!welltemp_) error("The calculation of the acceleration works only if Well-Tempered Metadynamics is on");
739 0 : log.printf(" calculation on the fly of the acceleration factor");
740 0 : addComponent("acc"); componentIsNotPeriodic("acc");
741 : }
742 :
743 : // for performance
744 94 : dp_ = new double[getNumberOfArguments()];
745 :
746 : // initializing and checking grid
747 94 : if(grid_) {
748 : // check for mesh and sigma size
749 99 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
750 : double a,b;
751 65 : Tools::convert(gmin[i],a);
752 65 : Tools::convert(gmax[i],b);
753 65 : double mesh=(b-a)/((double)gbin[i]);
754 65 : if(adaptive_==FlexibleBin::none) {
755 65 : if(mesh>0.5*sigma0_[i]) log<<" WARNING: Using a METAD with a Grid Spacing larger than half of the Gaussians width can produce artifacts\n";
756 : } else {
757 0 : if(mesh>0.5*sigma0min_[i]||sigma0min_[i]<0.) log<<" WARNING: to use a METAD with a GRID and ADAPTIVE you need to set a Grid Spacing larger than half of the Gaussians \n";
758 : }
759 : }
760 34 : std::string funcl=getLabel() + ".bias";
761 34 : if(!sparsegrid) {BiasGrid_=new Grid(funcl,getArguments(),gmin,gmax,gbin,spline,true);}
762 5 : else {BiasGrid_=new SparseGrid(funcl,getArguments(),gmin,gmax,gbin,spline,true);}
763 68 : std::vector<std::string> actualmin=BiasGrid_->getMin();
764 68 : std::vector<std::string> actualmax=BiasGrid_->getMax();
765 99 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
766 65 : if(gmin[i]!=actualmin[i]) log<<" WARNING: GRID_MIN["<<i<<"] has been adjusted to "<<actualmin[i]<<" to fit periodicity\n";
767 65 : if(gmax[i]!=actualmax[i]) log<<" WARNING: GRID_MAX["<<i<<"] has been adjusted to "<<actualmax[i]<<" to fit periodicity\n";
768 34 : }
769 : }
770 :
771 : // restart from external grid
772 94 : bool restartedFromGrid=false;
773 94 : if(gridreadfilename_.length()>0) {
774 : // read the grid in input, find the keys
775 3 : IFile gridfile;
776 3 : gridfile.link(*this);
777 3 : if(gridfile.FileExist(gridreadfilename_)) {
778 3 : gridfile.open(gridreadfilename_);
779 : } else {
780 0 : error("The GRID file you want to read: " + gridreadfilename_ + ", cannot be found!");
781 : }
782 6 : std::string funcl=getLabel() + ".bias";
783 3 : BiasGrid_=Grid::create(funcl, getArguments(), gridfile, gmin, gmax, gbin, sparsegrid, spline, true);
784 3 : gridfile.close();
785 3 : if(BiasGrid_->getDimension()!=getNumberOfArguments()) error("mismatch between dimensionality of input grid and number of arguments");
786 9 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
787 6 : if( getPntrToArgument(i)->isPeriodic()!=BiasGrid_->getIsPeriodic()[i] ) error("periodicity mismatch between arguments and input bias");
788 : double a, b;
789 6 : Tools::convert(gmin[i],a);
790 6 : Tools::convert(gmax[i],b);
791 6 : double mesh=(b-a)/((double)gbin[i]);
792 6 : if(mesh>0.5*sigma0_[i]) log<<" WARNING: Using a METAD with a Grid Spacing larger than half of the Gaussians width can produce artifacts\n";
793 : }
794 3 : log.printf(" Restarting from %s:",gridreadfilename_.c_str());
795 6 : if(getRestart()) restartedFromGrid=true;
796 : }
797 :
798 : // initializing and checking grid
799 94 : if(grid_&&!(gridreadfilename_.length()>0)) {
800 : // check for adaptive and sigma_min
801 31 : if(sigma0min_.size()==0&&adaptive_!=FlexibleBin::none) error("When using Adaptive Gaussians on a grid SIGMA_MIN must be specified");
802 : // check for mesh and sigma size
803 90 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
804 : double a,b;
805 59 : Tools::convert(gmin[i],a);
806 59 : Tools::convert(gmax[i],b);
807 59 : double mesh=(b-a)/((double)gbin[i]);
808 59 : if(mesh>0.5*sigma0_[i]) log<<" WARNING: Using a METAD with a Grid Spacing larger than half of the Gaussians width can produce artifacts\n";
809 : }
810 31 : std::string funcl=getLabel() + ".bias";
811 31 : if(!sparsegrid) {BiasGrid_=new Grid(funcl,getArguments(),gmin,gmax,gbin,spline,true);}
812 5 : else {BiasGrid_=new SparseGrid(funcl,getArguments(),gmin,gmax,gbin,spline,true);}
813 62 : std::vector<std::string> actualmin=BiasGrid_->getMin();
814 62 : std::vector<std::string> actualmax=BiasGrid_->getMax();
815 90 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
816 59 : if(gmin[i]!=actualmin[i]) log<<" WARNING: GRID_MIN["<<i<<"] has been adjusted to "<<actualmin[i]<<" to fit periodicity\n";
817 59 : if(gmax[i]!=actualmax[i]) log<<" WARNING: GRID_MAX["<<i<<"] has been adjusted to "<<actualmax[i]<<" to fit periodicity\n";
818 31 : }
819 : }
820 :
821 : // creating vector of ifile* for hills reading
822 : // open all files at the beginning and read Gaussians if restarting
823 200 : for(int i=0; i<mw_n_; ++i) {
824 106 : string fname;
825 106 : if(mw_n_>1) {
826 18 : stringstream out; out << i;
827 18 : fname = mw_dir_+"/"+hillsfname+"."+out.str();
828 : } else {
829 88 : fname = hillsfname;
830 : }
831 106 : IFile *ifile = new IFile();
832 106 : ifile->link(*this);
833 106 : ifiles.push_back(ifile);
834 106 : ifilesnames.push_back(fname);
835 106 : if(ifile->FileExist(fname)) {
836 34 : ifile->open(fname);
837 34 : if(getRestart()&&!restartedFromGrid) {
838 17 : log.printf(" Restarting from %s:",ifilesnames[i].c_str());
839 17 : readGaussians(ifiles[i]);
840 : }
841 34 : ifiles[i]->reset(false);
842 : // close only the walker own hills file for later writing
843 34 : if(i==mw_id_) ifiles[i]->close();
844 : }
845 106 : }
846 :
847 94 : comm.Barrier();
848 : // this barrier is needed when using walkers_mpi
849 : // to be sure that all files have been read before
850 : // backing them up
851 : // it should not be used when walkers_mpi is false otherwise
852 : // it would introduce troubles when using replicas without METAD
853 : // (e.g. in bias exchange with a neutral replica)
854 : // see issue #168 on github
855 94 : if(comm.Get_rank()==0 && walkers_mpi) multi_sim_comm.Barrier();
856 94 : if(targetfilename_.length()>0) {
857 1 : IFile gridfile; gridfile.open(targetfilename_);
858 2 : std::string funcl=getLabel() + ".target";
859 1 : TargetGrid_=Grid::create(funcl,getArguments(),gridfile,false,false,true);
860 1 : gridfile.close();
861 1 : if(TargetGrid_->getDimension()!=getNumberOfArguments()) error("mismatch between dimensionality of input grid and number of arguments");
862 2 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
863 1 : if( getPntrToArgument(i)->isPeriodic()!=TargetGrid_->getIsPeriodic()[i] ) error("periodicity mismatch between arguments and input bias");
864 1 : }
865 : }
866 :
867 : // Calculate the Tiwary-Parrinello reweighting factor if we are restarting from previous hills
868 94 : if(getRestart() && rewf_grid_.size()>0 ) computeReweightingFactor();
869 :
870 : // open grid file for writing
871 94 : if(wgridstride_>0) {
872 16 : gridfile_.link(*this);
873 16 : if(walkers_mpi) {
874 0 : int r=0;
875 0 : if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
876 0 : comm.Bcast(r,0);
877 0 : if(r>0) gridfilename_="/dev/null";
878 0 : gridfile_.enforceSuffix("");
879 : }
880 16 : if(mw_n_>1) gridfile_.enforceSuffix("");
881 16 : gridfile_.open(gridfilename_);
882 : }
883 :
884 : // open hills file for writing
885 94 : hillsOfile_.link(*this);
886 94 : if(walkers_mpi) {
887 18 : int r=0;
888 18 : if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
889 18 : comm.Bcast(r,0);
890 18 : if(r>0) ifilesnames[mw_id_]="/dev/null";
891 18 : hillsOfile_.enforceSuffix("");
892 : }
893 94 : if(mw_n_>1) hillsOfile_.enforceSuffix("");
894 94 : hillsOfile_.open(ifilesnames[mw_id_]);
895 94 : if(fmt.length()>0) hillsOfile_.fmtField(fmt);
896 94 : hillsOfile_.addConstantField("multivariate");
897 94 : if(doInt_) {
898 2 : hillsOfile_.addConstantField("lower_int").printField("lower_int",lowI_);
899 2 : hillsOfile_.addConstantField("upper_int").printField("upper_int",uppI_);
900 : }
901 94 : hillsOfile_.setHeavyFlush();
902 : // output periodicities of variables
903 94 : for(unsigned i=0; i<getNumberOfArguments(); ++i) hillsOfile_.setupPrintValue( getPntrToArgument(i) );
904 :
905 94 : bool concurrent=false;
906 94 : const ActionSet&actionSet(plumed.getActionSet());
907 94 : for(ActionSet::const_iterator p=actionSet.begin(); p!=actionSet.end(); ++p) if(dynamic_cast<MetaD*>(*p)) { concurrent=true; break; }
908 94 : if(concurrent) log<<" You are using concurrent metadynamics\n";
909 :
910 94 : log<<" Bibliography "<<plumed.cite("Laio and Parrinello, PNAS 99, 12562 (2002)");
911 94 : if(welltemp_) log<<plumed.cite(
912 12 : "Barducci, Bussi, and Parrinello, Phys. Rev. Lett. 100, 020603 (2008)");
913 94 : if(mw_n_>1||walkers_mpi) log<<plumed.cite(
914 24 : "Raiteri, Laio, Gervasio, Micheletti, and Parrinello, J. Phys. Chem. B 110, 3533 (2006)");
915 94 : if(adaptive_!=FlexibleBin::none) log<<plumed.cite(
916 21 : "Branduardi, Bussi, and Parrinello, J. Chem. Theory Comput. 8, 2247 (2012)");
917 94 : if(doInt_) log<<plumed.cite(
918 2 : "Baftizadeh, Cossio, Pietrucci, and Laio, Curr. Phys. Chem. 2, 79 (2012)");
919 94 : if(acceleration) log<<plumed.cite(
920 0 : "Pratyush and Parrinello, Phys. Rev. Lett. 111, 230602 (2013)");
921 94 : if(rewf_grid_.size()>0) log<<plumed.cite(
922 5 : "Pratyush and Parrinello, J. Phys. Chem. B, 119, 736 (2015)");
923 94 : if(concurrent) log<<plumed.cite(
924 47 : "Gil-Ley and Bussi, J. Chem. Theory Comput. 11, 1077 (2015)");
925 94 : if(targetfilename_.length()>0) {
926 1 : log<<plumed.cite("White, Dama, and Voth, J. Chem. Theory Comput. 11, 2451 (2015)");
927 1 : log<<plumed.cite("Marinelli and Faraldo-Gómez, Biophys. J. 108, 2779 (2015)");
928 1 : log<<plumed.cite("Gil-Ley, Bottaro, and Bussi, J. Chem. Theory Comput. 12, 2790 (2016)");
929 : }
930 188 : log<<"\n";
931 94 : }
932 :
933 6035 : void MetaD::readGaussians(IFile *ifile)
934 : {
935 6035 : unsigned ncv=getNumberOfArguments();
936 6035 : vector<double> center(ncv);
937 12070 : vector<double> sigma(ncv);
938 : double height;
939 6035 : int nhills=0;
940 6035 : bool multivariate=false;
941 :
942 12070 : std::vector<Value> tmpvalues;
943 6035 : for(unsigned j=0; j<getNumberOfArguments(); ++j) tmpvalues.push_back( Value( this, getPntrToArgument(j)->getName(), false ) );
944 :
945 14989 : while(scanOneHill(ifile,tmpvalues,center,sigma,height,multivariate)) {
946 : ;
947 2919 : nhills++;
948 2919 : if(welltemp_) {height*=(biasf_-1.0)/biasf_;}
949 2919 : addGaussian(Gaussian(center,sigma,height,multivariate));
950 : }
951 12070 : log.printf(" %d Gaussians read\n",nhills);
952 6035 : }
953 :
954 0 : bool MetaD::readChunkOfGaussians(IFile *ifile, unsigned n)
955 : {
956 0 : unsigned ncv=getNumberOfArguments();
957 0 : vector<double> center(ncv);
958 0 : vector<double> sigma(ncv);
959 : double height;
960 0 : unsigned nhills=0;
961 0 : bool multivariate=false;
962 0 : std::vector<Value> tmpvalues;
963 0 : for(unsigned j=0; j<getNumberOfArguments(); ++j) tmpvalues.push_back( Value( this, getPntrToArgument(j)->getName(), false ) );
964 :
965 0 : while(scanOneHill(ifile,tmpvalues,center,sigma,height,multivariate)) {
966 : ;
967 0 : if(welltemp_) height*=(biasf_-1.0)/biasf_;
968 0 : addGaussian(Gaussian(center,sigma,height,multivariate));
969 0 : if(nhills==n) {
970 0 : log.printf(" %u Gaussians read\n",nhills);
971 0 : return true;
972 : }
973 0 : nhills++;
974 : }
975 0 : log.printf(" %u Gaussians read\n",nhills);
976 0 : return false;
977 : }
978 :
979 2363 : void MetaD::writeGaussian(const Gaussian& hill, OFile&file)
980 : {
981 2363 : unsigned ncv=getNumberOfArguments();
982 2363 : file.printField("time",getTimeStep()*getStep());
983 6772 : for(unsigned i=0; i<ncv; ++i) {
984 4409 : file.printField(getPntrToArgument(i),hill.center[i]);
985 : }
986 2363 : if(hill.multivariate) {
987 446 : hillsOfile_.printField("multivariate","true");
988 446 : Matrix<double> mymatrix(ncv,ncv);
989 446 : unsigned k=0;
990 1047 : for(unsigned i=0; i<ncv; i++) {
991 1357 : for(unsigned j=i; j<ncv; j++) {
992 : // recompose the full inverse matrix
993 756 : mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k];
994 756 : k++;
995 : }
996 : }
997 : // invert it
998 892 : Matrix<double> invmatrix(ncv,ncv);
999 446 : Invert(mymatrix,invmatrix);
1000 : // enforce symmetry
1001 1047 : for(unsigned i=0; i<ncv; i++) {
1002 1357 : for(unsigned j=i; j<ncv; j++) {
1003 756 : invmatrix(i,j)=invmatrix(j,i);
1004 : }
1005 : }
1006 :
1007 : // do cholesky so to have a "sigma like" number
1008 892 : Matrix<double> lower(ncv,ncv);
1009 446 : cholesky(invmatrix,lower);
1010 : // loop in band form
1011 1047 : for(unsigned i=0; i<ncv; i++) {
1012 1357 : for(unsigned j=0; j<ncv-i; j++) {
1013 756 : file.printField("sigma_"+getPntrToArgument(j+i)->getName()+"_"+getPntrToArgument(j)->getName(),lower(j+i,j));
1014 : }
1015 446 : }
1016 : } else {
1017 1917 : hillsOfile_.printField("multivariate","false");
1018 5725 : for(unsigned i=0; i<ncv; ++i)
1019 3808 : file.printField("sigma_"+getPntrToArgument(i)->getName(),hill.sigma[i]);
1020 : }
1021 2363 : double height=hill.height;
1022 2363 : if(welltemp_) height*=biasf_/(biasf_-1.0);
1023 2363 : file.printField("height",height).printField("biasf",biasf_);
1024 2363 : if(mw_n_>1) file.printField("clock",int(std::time(0)));
1025 2363 : file.printField();
1026 2363 : }
1027 :
1028 5282 : void MetaD::addGaussian(const Gaussian& hill)
1029 : {
1030 5282 : if(!grid_) hills_.push_back(hill);
1031 : else {
1032 217 : unsigned ncv=getNumberOfArguments();
1033 217 : vector<unsigned> nneighb=getGaussianSupport(hill);
1034 434 : vector<Grid::index_t> neighbors=BiasGrid_->getNeighbors(hill.center,nneighb);
1035 434 : vector<double> der(ncv);
1036 434 : vector<double> xx(ncv);
1037 217 : if(comm.Get_size()==1) {
1038 45580 : for(unsigned i=0; i<neighbors.size(); ++i) {
1039 45451 : Grid::index_t ineigh=neighbors[i];
1040 45451 : for(unsigned j=0; j<ncv; ++j) der[j]=0.0;
1041 45451 : BiasGrid_->getPoint(ineigh,xx);
1042 45451 : double bias=evaluateGaussian(xx,hill,&der[0]);
1043 45451 : BiasGrid_->addValueAndDerivatives(ineigh,bias,der);
1044 : }
1045 : } else {
1046 88 : unsigned stride=comm.Get_size();
1047 88 : unsigned rank=comm.Get_rank();
1048 88 : vector<double> allder(ncv*neighbors.size(),0.0);
1049 176 : vector<double> allbias(neighbors.size(),0.0);
1050 27104 : for(unsigned i=rank; i<neighbors.size(); i+=stride) {
1051 27016 : Grid::index_t ineigh=neighbors[i];
1052 27016 : BiasGrid_->getPoint(ineigh,xx);
1053 27016 : allbias[i]=evaluateGaussian(xx,hill,&allder[ncv*i]);
1054 : }
1055 88 : comm.Sum(allbias);
1056 88 : comm.Sum(allder);
1057 103120 : for(unsigned i=0; i<neighbors.size(); ++i) {
1058 103032 : Grid::index_t ineigh=neighbors[i];
1059 103032 : for(unsigned j=0; j<ncv; ++j) {der[j]=allder[ncv*i+j];}
1060 103032 : BiasGrid_->addValueAndDerivatives(ineigh,allbias[i],der);
1061 88 : }
1062 217 : }
1063 : }
1064 5282 : }
1065 :
1066 217 : vector<unsigned> MetaD::getGaussianSupport(const Gaussian& hill)
1067 : {
1068 217 : vector<unsigned> nneigh;
1069 434 : vector<double> cutoff;
1070 217 : unsigned ncv=getNumberOfArguments();
1071 :
1072 : // traditional or flexible hill?
1073 217 : if(hill.multivariate) {
1074 0 : unsigned k=0;
1075 0 : Matrix<double> mymatrix(ncv,ncv);
1076 0 : for(unsigned i=0; i<ncv; i++) {
1077 0 : for(unsigned j=i; j<ncv; j++) {
1078 : // recompose the full inverse matrix
1079 0 : mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k];
1080 0 : k++;
1081 : }
1082 : }
1083 : // Reinvert so to have the ellipses
1084 0 : Matrix<double> myinv(ncv,ncv);
1085 0 : Invert(mymatrix,myinv);
1086 0 : Matrix<double> myautovec(ncv,ncv);
1087 0 : vector<double> myautoval(ncv); //should I take this or their square root?
1088 0 : diagMat(myinv,myautoval,myautovec);
1089 0 : double maxautoval=0.;
1090 0 : unsigned ind_maxautoval; ind_maxautoval=ncv;
1091 0 : for(unsigned i=0; i<ncv; i++) {
1092 0 : if(myautoval[i]>maxautoval) {maxautoval=myautoval[i]; ind_maxautoval=i;}
1093 : }
1094 0 : for(unsigned i=0; i<ncv; i++) {
1095 0 : cutoff.push_back(sqrt(2.0*DP2CUTOFF)*abs(sqrt(maxautoval)*myautovec(i,ind_maxautoval)));
1096 0 : }
1097 : } else {
1098 638 : for(unsigned i=0; i<ncv; ++i) {
1099 421 : cutoff.push_back(sqrt(2.0*DP2CUTOFF)*hill.sigma[i]);
1100 : }
1101 : }
1102 :
1103 217 : if(doInt_) {
1104 2 : if(hill.center[0]+cutoff[0] > uppI_ || hill.center[0]-cutoff[0] < lowI_) {
1105 : // in this case, we updated the entire grid to avoid problems
1106 2 : return BiasGrid_->getNbin();
1107 : } else {
1108 0 : nneigh.push_back( static_cast<unsigned>(ceil(cutoff[0]/BiasGrid_->getDx()[0])) );
1109 0 : return nneigh;
1110 : }
1111 : } else {
1112 634 : for(unsigned i=0; i<ncv; i++) {
1113 419 : nneigh.push_back( static_cast<unsigned>(ceil(cutoff[i]/BiasGrid_->getDx()[i])) );
1114 : }
1115 : }
1116 :
1117 432 : return nneigh;
1118 : }
1119 :
1120 916305 : double MetaD::getBiasAndDerivatives(const vector<double>& cv, double* der)
1121 : {
1122 916305 : double bias=0.0;
1123 916305 : if(!grid_) {
1124 13957 : if(hills_.size()>10000 && (getStep()-last_step_warn_grid)>10000) {
1125 0 : std::string msg;
1126 0 : Tools::convert(hills_.size(),msg);
1127 0 : msg="You have accumulated "+msg+" hills, you should enable GRIDs to avoid serious performance hits";
1128 0 : warning(msg);
1129 0 : last_step_warn_grid=getStep();
1130 : }
1131 13957 : unsigned stride=comm.Get_size();
1132 13957 : unsigned rank=comm.Get_rank();
1133 6972246 : for(unsigned i=rank; i<hills_.size(); i+=stride) {
1134 6958289 : bias+=evaluateGaussian(cv,hills_[i],der);
1135 : //finite difference test
1136 : //finiteDifferenceGaussian(cv,hills_[i]);
1137 : }
1138 13957 : comm.Sum(bias);
1139 13957 : if(der) comm.Sum(der,getNumberOfArguments());
1140 : } else {
1141 902348 : if(der) {
1142 900740 : vector<double> vder(getNumberOfArguments());
1143 900740 : bias=BiasGrid_->getValueAndDerivatives(cv,vder);
1144 900740 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {der[i]=vder[i];}
1145 : } else {
1146 1608 : bias = BiasGrid_->getValue(cv);
1147 : }
1148 : }
1149 :
1150 916305 : return bias;
1151 : }
1152 :
1153 0 : double MetaD::getGaussianNormalization( const Gaussian& hill )
1154 : {
1155 0 : double norm=1;
1156 0 : unsigned ncv=hill.center.size();
1157 :
1158 0 : if(hill.multivariate) {
1159 : // recompose the full sigma from the upper diag cholesky
1160 0 : unsigned k=0;
1161 0 : Matrix<double> mymatrix(ncv,ncv);
1162 0 : for(unsigned i=0; i<ncv; i++) {
1163 0 : for(unsigned j=i; j<ncv; j++) {
1164 0 : mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k]; // recompose the full inverse matrix
1165 0 : k++;
1166 : }
1167 0 : double ldet; logdet( mymatrix, ldet );
1168 0 : norm = exp( ldet ); // Not sure here if mymatrix is sigma or inverse
1169 0 : }
1170 : } else {
1171 0 : for(unsigned i=0; i<hill.sigma.size(); ++i) norm*=hill.sigma[i];
1172 : }
1173 :
1174 0 : return norm*pow(2*pi,static_cast<double>(ncv)/2.0);
1175 : }
1176 :
1177 7030756 : double MetaD::evaluateGaussian(const vector<double>& cv, const Gaussian& hill, double* der)
1178 : {
1179 7030756 : double dp2=0.0;
1180 7030756 : double bias=0.0;
1181 : // I use a pointer here because cv is const (and should be const)
1182 : // but when using doInt it is easier to locally replace cv[0] with
1183 : // the upper/lower limit in case it is out of range
1184 7030756 : const double *pcv=NULL; // pointer to cv
1185 : double tmpcv[1]; // tmp array with cv (to be used with doInt_)
1186 7030756 : if(cv.size()>0) pcv=&cv[0];
1187 7030756 : if(doInt_) {
1188 1402 : plumed_assert(cv.size()==1);
1189 1402 : tmpcv[0]=cv[0];
1190 1402 : if(cv[0]<lowI_) tmpcv[0]=lowI_;
1191 1402 : if(cv[0]>uppI_) tmpcv[0]=uppI_;
1192 1402 : pcv=&(tmpcv[0]);
1193 : }
1194 7030756 : if(hill.multivariate) {
1195 230564 : unsigned k=0;
1196 230564 : unsigned ncv=cv.size();
1197 : // recompose the full sigma from the upper diag cholesky
1198 230564 : Matrix<double> mymatrix(ncv,ncv);
1199 462552 : for(unsigned i=0; i<ncv; i++) {
1200 465400 : for(unsigned j=i; j<ncv; j++) {
1201 233412 : mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k]; // recompose the full inverse matrix
1202 233412 : k++;
1203 : }
1204 : }
1205 462552 : for(unsigned i=0; i<cv.size(); ++i) {
1206 231988 : double dp_i=difference(i,hill.center[i],pcv[i]);
1207 231988 : dp_[i]=dp_i;
1208 465400 : for(unsigned j=i; j<cv.size(); ++j) {
1209 233412 : if(i==j) {
1210 231988 : dp2+=dp_i*dp_i*mymatrix(i,j)*0.5;
1211 : } else {
1212 1424 : double dp_j=difference(j,hill.center[j],pcv[j]);
1213 1424 : dp2+=dp_i*dp_j*mymatrix(i,j);
1214 : }
1215 : }
1216 : }
1217 230564 : if(dp2<DP2CUTOFF) {
1218 221813 : bias=hill.height*exp(-dp2);
1219 221813 : if(der) {
1220 155673 : for(unsigned i=0; i<cv.size(); ++i) {
1221 77990 : double tmp=0.0;
1222 156594 : for(unsigned j=0; j<cv.size(); ++j) {
1223 78604 : tmp += dp_[j]*mymatrix(i,j)*bias;
1224 : }
1225 77990 : der[i]-=tmp;
1226 : }
1227 : }
1228 230564 : }
1229 : } else {
1230 20397449 : for(unsigned i=0; i<cv.size(); ++i) {
1231 13597257 : double dp=difference(i,hill.center[i],pcv[i])*hill.invsigma[i];
1232 13597257 : dp2+=dp*dp;
1233 13597257 : dp_[i]=dp;
1234 : }
1235 6800192 : dp2*=0.5;
1236 6800192 : if(dp2<DP2CUTOFF) {
1237 3939855 : bias=hill.height*exp(-dp2);
1238 3939855 : if(der) {
1239 1347065 : for(unsigned i=0; i<cv.size(); ++i) {der[i]+=-bias*dp_[i]*hill.invsigma[i];}
1240 : }
1241 : }
1242 : }
1243 :
1244 7030756 : if(doInt_ && der) {
1245 602 : if(cv[0]<lowI_ || cv[0]>uppI_) for(unsigned i=0; i<cv.size(); ++i) der[i]=0;
1246 : }
1247 :
1248 7030756 : return bias;
1249 : }
1250 :
1251 2219 : double MetaD::getHeight(const vector<double>& cv)
1252 : {
1253 2219 : double height=height0_;
1254 2219 : if(welltemp_) {
1255 128 : double vbias = getBiasAndDerivatives(cv);
1256 128 : height = height0_*exp(-vbias/(kbt_*(biasf_-1.0)));
1257 : }
1258 2219 : if(dampfactor_>0.0) {
1259 9 : plumed_assert(BiasGrid_);
1260 9 : double m=BiasGrid_->getMaxValue();
1261 9 : height*=exp(-m/(kbt_*(dampfactor_)));
1262 : }
1263 2219 : if(TargetGrid_) {
1264 9 : double f=TargetGrid_->getValue(cv)-TargetGrid_->getMaxValue();
1265 9 : height*=exp(f/kbt_);
1266 : }
1267 2219 : return height;
1268 : }
1269 :
1270 5743 : void MetaD::calculate()
1271 : {
1272 : // this is because presently there is no way to properly pass information
1273 : // on adaptive hills (diff) after exchanges:
1274 5743 : if(adaptive_==FlexibleBin::diffusion && getExchangeStep()) error("ADAPTIVE=DIFF is not compatible with replica exchange");
1275 :
1276 5743 : const unsigned ncv=getNumberOfArguments();
1277 5743 : vector<double> cv(ncv);
1278 5743 : double* der = new double[ncv];
1279 15297 : for(unsigned i=0; i<ncv; ++i) {
1280 9554 : cv[i]=getArgument(i);
1281 9554 : der[i]=0.;
1282 : }
1283 5743 : const double ene = getBiasAndDerivatives(cv,der);
1284 5743 : setBias(ene);
1285 5743 : if( rewf_grid_.size()>0 ) getPntrToComponent("rbias")->set(ene - reweight_factor);
1286 : // calculate the acceleration factor
1287 5743 : if(acceleration&&!isFirstStep) {
1288 0 : acc += exp(ene/(kbt_));
1289 0 : const double mean_acc = acc/((double) getStep());
1290 0 : getPntrToComponent("acc")->set(mean_acc);
1291 : }
1292 5743 : getPntrToComponent("work")->set(work_);
1293 : // set Forces
1294 15297 : for(unsigned i=0; i<ncv; ++i) {
1295 9554 : setOutputForce(i,-der[i]);
1296 : }
1297 5743 : delete [] der;
1298 5743 : }
1299 :
1300 5217 : void MetaD::update() {
1301 5217 : vector<double> cv(getNumberOfArguments());
1302 10434 : vector<double> thissigma;
1303 : bool multivariate;
1304 :
1305 : // adding hills criteria (could be more complex though)
1306 : bool nowAddAHill;
1307 5217 : if(getStep()%stride_==0 && !isFirstStep )nowAddAHill=true;
1308 : else {
1309 2998 : nowAddAHill=false;
1310 2998 : isFirstStep=false;
1311 : }
1312 :
1313 5217 : for(unsigned i=0; i<cv.size(); ++i) cv[i] = getArgument(i);
1314 :
1315 5217 : double vbias=getBiasAndDerivatives(cv);
1316 :
1317 : // if you use adaptive, call the FlexibleBin
1318 5217 : if(adaptive_!=FlexibleBin::none) {
1319 778 : flexbin->update(nowAddAHill);
1320 778 : multivariate=true;
1321 : } else {
1322 4439 : multivariate=false;
1323 : }
1324 :
1325 5217 : if(nowAddAHill) {
1326 : // add a Gaussian
1327 2219 : double height=getHeight(cv);
1328 : // returns upper diagonal inverse
1329 2219 : if(adaptive_!=FlexibleBin::none) thissigma=flexbin->getInverseMatrix();
1330 : // returns normal sigma
1331 1845 : else thissigma=sigma0_;
1332 :
1333 : // In case we use walkers_mpi, it is now necessary to communicate with other replicas.
1334 2219 : if(walkers_mpi) {
1335 : // Allocate arrays to store all walkers hills
1336 72 : std::vector<double> all_cv(mpi_nw_*cv.size(),0.0);
1337 144 : std::vector<double> all_sigma(mpi_nw_*thissigma.size(),0.0);
1338 144 : std::vector<double> all_height(mpi_nw_,0.0);
1339 144 : std::vector<int> all_multivariate(mpi_nw_,0);
1340 72 : if(comm.Get_rank()==0) {
1341 : // Communicate (only root)
1342 48 : multi_sim_comm.Allgather(cv,all_cv);
1343 48 : multi_sim_comm.Allgather(thissigma,all_sigma);
1344 48 : multi_sim_comm.Allgather(height,all_height);
1345 48 : multi_sim_comm.Allgather(int(multivariate),all_multivariate);
1346 : }
1347 : // Share info with group members
1348 72 : comm.Bcast(all_cv,0);
1349 72 : comm.Bcast(all_sigma,0);
1350 72 : comm.Bcast(all_height,0);
1351 72 : comm.Bcast(all_multivariate,0);
1352 288 : for(unsigned i=0; i<mpi_nw_; i++) {
1353 : // actually add hills one by one
1354 216 : std::vector<double> cv_now(cv.size());
1355 432 : std::vector<double> sigma_now(thissigma.size());
1356 216 : for(unsigned j=0; j<cv.size(); j++) cv_now[j]=all_cv[i*cv.size()+j];
1357 216 : for(unsigned j=0; j<thissigma.size(); j++) sigma_now[j]=all_sigma[i*thissigma.size()+j];
1358 432 : Gaussian newhill=Gaussian(cv_now,sigma_now,all_height[i],all_multivariate[i]);
1359 216 : addGaussian(newhill);
1360 216 : writeGaussian(newhill,hillsOfile_);
1361 288 : }
1362 : } else {
1363 2147 : Gaussian newhill=Gaussian(cv,thissigma,height,multivariate);
1364 2147 : addGaussian(newhill);
1365 : // print on HILLS file
1366 2147 : writeGaussian(newhill,hillsOfile_);
1367 : }
1368 : }
1369 :
1370 : // this should be outside of the if block in case
1371 : // mw_rstride_ is not a multiple of stride_
1372 5217 : if(mw_n_>1 && getStep()%mw_rstride_==0) {
1373 3012 : hillsOfile_.flush();
1374 : }
1375 :
1376 5217 : double vbias1=getBiasAndDerivatives(cv);
1377 5217 : work_+=vbias1-vbias;
1378 :
1379 : // dump grid on file
1380 5217 : if(wgridstride_>0&&(getStep()%wgridstride_==0||getCPT())) {
1381 : // in case old grids are stored, a sequence of grids should appear
1382 : // this call results in a repetition of the header:
1383 80 : if(storeOldGrids_) gridfile_.clearFields();
1384 : // in case only latest grid is stored, file should be rewound
1385 : // this will overwrite previously written grids
1386 : else {
1387 40 : int r = 0;
1388 40 : if(walkers_mpi) {
1389 0 : if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
1390 0 : comm.Bcast(r,0);
1391 : }
1392 40 : if(r==0) gridfile_.rewind();
1393 : }
1394 80 : BiasGrid_->writeToFile(gridfile_);
1395 : // if a single grid is stored, it is necessary to flush it, otherwise
1396 : // the file might stay empty forever (when a single grid is not large enough to
1397 : // trigger flushing from the operating system).
1398 : // on the other hand, if grids are stored one after the other this is
1399 : // no necessary, and we leave the flushing control to the user as usual
1400 : // (with FLUSH keyword)
1401 80 : if(!storeOldGrids_) gridfile_.flush();
1402 : }
1403 :
1404 : // if multiple walkers and time to read Gaussians
1405 5217 : if(mw_n_>1 && getStep()%mw_rstride_==0) {
1406 12048 : for(int i=0; i<mw_n_; ++i) {
1407 : // don't read your own Gaussians
1408 9036 : if(i==mw_id_) continue;
1409 : // if the file is not open yet
1410 6024 : if(!(ifiles[i]->isOpen())) {
1411 : // check if it exists now and open it!
1412 6 : if(ifiles[i]->FileExist(ifilesnames[i])) {
1413 6 : ifiles[i]->open(ifilesnames[i]);
1414 6 : ifiles[i]->reset(false);
1415 : }
1416 : // otherwise read the new Gaussians
1417 : } else {
1418 6018 : log.printf(" Reading hills from %s:",ifilesnames[i].c_str());
1419 6018 : readGaussians(ifiles[i]);
1420 6018 : ifiles[i]->reset(false);
1421 : }
1422 : }
1423 : }
1424 10434 : if(getStep()%(stride_*rewf_ustride_)==0 && nowAddAHill && rewf_grid_.size()>0 ) computeReweightingFactor();
1425 5217 : }
1426 :
1427 0 : void MetaD::finiteDifferenceGaussian(const vector<double>& cv, const Gaussian& hill)
1428 : {
1429 0 : log<<"--------- finiteDifferenceGaussian: size "<<cv.size() <<"------------\n";
1430 : // for each cv
1431 : // first get the bias and the derivative
1432 0 : vector<double> oldder(cv.size());
1433 0 : vector<double> der(cv.size());
1434 0 : vector<double> mycv(cv.size());
1435 0 : mycv=cv;
1436 0 : double step=1.e-6;
1437 0 : Random random;
1438 : // just displace a tiny bit
1439 0 : for(unsigned i=0; i<cv.size(); i++) log<<"CV "<<i<<" V "<<mycv[i]<<"\n";
1440 0 : for(unsigned i=0; i<cv.size(); i++) mycv[i]+=1.e-2*2*(random.RandU01()-0.5);
1441 0 : for(unsigned i=0; i<cv.size(); i++) log<<"NENEWWCV "<<i<<" V "<<mycv[i]<<"\n";
1442 0 : double oldbias=evaluateGaussian(mycv,hill,&oldder[0]);
1443 0 : for(unsigned i=0; i<mycv.size(); i++) {
1444 0 : double delta=step*2*(random.RandU01()-0.5);
1445 0 : mycv[i]+=delta;
1446 0 : double newbias=evaluateGaussian(mycv,hill,&der[0]);
1447 0 : log<<"CV "<<i;
1448 0 : log<<" ANAL "<<oldder[i]<<" NUM "<<(newbias-oldbias)/delta<<" DIFF "<<(oldder[i]-(newbias-oldbias)/delta)<<"\n";
1449 0 : mycv[i]-=delta;
1450 : }
1451 0 : log<<"--------- END finiteDifferenceGaussian ------------\n";
1452 0 : }
1453 :
1454 : /// takes a pointer to the file and a template string with values v and gives back the next center, sigma and height
1455 8954 : bool MetaD::scanOneHill(IFile *ifile, vector<Value> &tmpvalues, vector<double> ¢er, vector<double> &sigma, double &height, bool &multivariate)
1456 : {
1457 : double dummy;
1458 8954 : multivariate=false;
1459 8954 : if(ifile->scanField("time",dummy)) {
1460 2919 : unsigned ncv; ncv=tmpvalues.size();
1461 8752 : for(unsigned i=0; i<ncv; ++i) {
1462 5833 : ifile->scanField( &tmpvalues[i] );
1463 5833 : if( tmpvalues[i].isPeriodic() && ! getPntrToArgument(i)->isPeriodic() ) {
1464 0 : error("in hills file periodicity for variable " + tmpvalues[i].getName() + " does not match periodicity in input");
1465 5833 : } else if( tmpvalues[i].isPeriodic() ) {
1466 0 : std::string imin, imax; tmpvalues[i].getDomain( imin, imax );
1467 0 : std::string rmin, rmax; getPntrToArgument(i)->getDomain( rmin, rmax );
1468 0 : if( imin!=rmin || imax!=rmax ) {
1469 0 : error("in hills file periodicity for variable " + tmpvalues[i].getName() + " does not match periodicity in input");
1470 0 : }
1471 : }
1472 5833 : center[i]=tmpvalues[i].get();
1473 : }
1474 : // scan for multivariate label: record the actual file position so to eventually rewind
1475 2919 : std::string sss;
1476 2919 : ifile->scanField("multivariate",sss);
1477 2919 : if(sss=="true") multivariate=true;
1478 2919 : else if(sss=="false") multivariate=false;
1479 0 : else plumed_merror("cannot parse multivariate = "+ sss);
1480 2919 : if(multivariate) {
1481 0 : sigma.resize(ncv*(ncv+1)/2);
1482 0 : Matrix<double> upper(ncv,ncv);
1483 0 : Matrix<double> lower(ncv,ncv);
1484 0 : for(unsigned i=0; i<ncv; i++) {
1485 0 : for(unsigned j=0; j<ncv-i; j++) {
1486 0 : ifile->scanField("sigma_"+getPntrToArgument(j+i)->getName()+"_"+getPntrToArgument(j)->getName(),lower(j+i,j));
1487 0 : upper(j,j+i)=lower(j+i,j);
1488 : }
1489 : }
1490 0 : Matrix<double> mymult(ncv,ncv);
1491 0 : Matrix<double> invmatrix(ncv,ncv);
1492 0 : mult(lower,upper,mymult);
1493 : // now invert and get the sigmas
1494 0 : Invert(mymult,invmatrix);
1495 : // put the sigmas in the usual order: upper diagonal (this time in normal form and not in band form)
1496 0 : unsigned k=0;
1497 0 : for(unsigned i=0; i<ncv; i++) {
1498 0 : for(unsigned j=i; j<ncv; j++) {
1499 0 : sigma[k]=invmatrix(i,j);
1500 0 : k++;
1501 : }
1502 0 : }
1503 : } else {
1504 8752 : for(unsigned i=0; i<ncv; ++i) {
1505 5833 : ifile->scanField("sigma_"+getPntrToArgument(i)->getName(),sigma[i]);
1506 : }
1507 : }
1508 :
1509 2919 : ifile->scanField("height",height);
1510 2919 : ifile->scanField("biasf",dummy);
1511 2919 : if(ifile->FieldExist("clock")) ifile->scanField("clock",dummy);
1512 2919 : if(ifile->FieldExist("lower_int")) ifile->scanField("lower_int",dummy);
1513 2919 : if(ifile->FieldExist("upper_int")) ifile->scanField("upper_int",dummy);
1514 2919 : ifile->scanField();
1515 2919 : return true;
1516 : } else {
1517 6035 : return false;
1518 : }
1519 : }
1520 :
1521 100 : void MetaD::computeReweightingFactor()
1522 : {
1523 100 : if( !welltemp_ ) error("cannot compute the c(t) reweighting factors for non well-tempered metadynamics");
1524 :
1525 : // Recover the minimum values for the grid
1526 100 : unsigned ncv=getNumberOfArguments();
1527 100 : unsigned ntotgrid=1;
1528 200 : std::vector<double> dmin( ncv ),dmax( ncv ), grid_spacing( ncv ), vals( ncv );
1529 300 : for(unsigned j=0; j<ncv; ++j) {
1530 200 : Tools::convert( BiasGrid_->getMin()[j], dmin[j] );
1531 200 : Tools::convert( BiasGrid_->getMax()[j], dmax[j] );
1532 200 : grid_spacing[j] = ( dmax[j] - dmin[j] ) / static_cast<double>( rewf_grid_[j] );
1533 200 : if( !getPntrToArgument(j)->isPeriodic() ) dmax[j] += grid_spacing[j];
1534 200 : ntotgrid *= rewf_grid_[j];
1535 : }
1536 :
1537 : // Now sum over whole grid
1538 200 : reweight_factor=0.0; double* der=new double[ncv]; std::vector<unsigned> t_index( ncv );
1539 100 : double sum1=0.0; double sum2=0.0;
1540 100 : double afactor = biasf_ / (kbt_*(biasf_-1.0)); double afactor2 = 1.0 / (kbt_*(biasf_-1.0));
1541 100 : unsigned rank=comm.Get_rank(), stride=comm.Get_size();
1542 900100 : for(unsigned i=rank; i<ntotgrid; i+=stride) {
1543 900000 : t_index[0]=(i%rewf_grid_[0]);
1544 900000 : unsigned kk=i;
1545 900000 : for(unsigned j=1; j<ncv-1; ++j) { kk=(kk-t_index[j-1])/rewf_grid_[j-1]; t_index[j]=(kk%rewf_grid_[j]); }
1546 900000 : if( ncv>=2 ) t_index[ncv-1]=((kk-t_index[ncv-2])/rewf_grid_[ncv-2]);
1547 :
1548 900000 : for(unsigned j=0; j<ncv; ++j) vals[j]=dmin[j] + t_index[j]*grid_spacing[j];
1549 :
1550 900000 : double currentb=getBiasAndDerivatives(vals,der);
1551 900000 : sum1 += exp( afactor*currentb );
1552 900000 : sum2 += exp( afactor2*currentb );
1553 : }
1554 100 : delete [] der;
1555 100 : comm.Sum( sum1 ); comm.Sum( sum2 );
1556 100 : reweight_factor = kbt_ * std::log( sum1/sum2 );
1557 200 : getPntrToComponent("rct")->set(reweight_factor);
1558 100 : }
1559 :
1560 : }
1561 2523 : }
|
