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 "ActionRegister.h"
24 : #include "core/ActionSet.h"
25 : #include "core/PlumedMain.h"
26 : #include "core/Atoms.h"
27 : #include "core/FlexibleBin.h"
28 : #include "tools/Exception.h"
29 : #include "tools/Grid.h"
30 : #include "tools/Matrix.h"
31 : #include "tools/OpenMP.h"
32 : #include "tools/Random.h"
33 : #include "tools/File.h"
34 : #include <ctime>
35 : #include <numeric>
36 : #if defined(__PLUMED_HAS_GETCWD)
37 : #include <unistd.h>
38 : #endif
39 :
40 : namespace PLMD {
41 : namespace bias {
42 :
43 : //+PLUMEDOC BIAS PBMETAD
44 : /*
45 : Used to performed Parallel Bias metadynamics.
46 :
47 : This action activate Parallel Bias Metadynamics (PBMetaD) \cite pbmetad, a version of metadynamics \cite metad in which
48 : multiple low-dimensional bias potentials are applied in parallel.
49 : In the current implementation, these have the form of mono-dimensional metadynamics bias
50 : potentials:
51 :
52 : \f[
53 : {V(s_1,t), ..., V(s_N,t)}
54 : \f]
55 :
56 : where:
57 :
58 : \f[
59 : V(s_i,t) = \sum_{ k \tau < t} W_i(k \tau)
60 : \exp\left(
61 : - \frac{(s_i-s_i^{(0)}(k \tau))^2}{2\sigma_i^2}
62 : \right).
63 : \f]
64 :
65 : To ensure the convergence of each mono-dimensional bias potential to the corresponding free energy,
66 : at each deposition step the Gaussian heights are multiplied by the so-called conditional term:
67 :
68 : \f[
69 : W_i(k \tau)=W_0 \frac{\exp\left(
70 : - \frac{V(s_i,k \tau)}{k_B T}
71 : \right)}{\sum_{i=1}^N
72 : \exp\left(
73 : - \frac{V(s_i,k \tau)}{k_B T}
74 : \right)}
75 : \f]
76 :
77 : where \f$W_0\f$ is the initial Gaussian height.
78 :
79 : The PBMetaD bias potential is defined by:
80 :
81 : \f[
82 : V_{PB}(\vec{s},t) = -k_B T \log{\sum_{i=1}^N
83 : \exp\left(
84 : - \frac{V(s_i,t)}{k_B T}
85 : \right)}.
86 : \f]
87 :
88 :
89 : Information on the Gaussian functions that build each bias potential are printed to
90 : multiple HILLS files, which
91 : are used both to restart the calculation and to reconstruct the mono-dimensional
92 : free energies as a function of the corresponding CVs.
93 : These can be reconstructed using the \ref sum_hills utility because the final bias is given
94 : by:
95 :
96 : \f[
97 : V(s_i) = -F(s_i)
98 : \f]
99 :
100 : Currently, only a subset of the \ref METAD options are available in PBMetaD.
101 :
102 : The bias potentials can be stored on a grid to increase performances of long PBMetaD simulations.
103 : You should
104 : provide either the number of bins for every collective variable (GRID_BIN) or
105 : the desired grid spacing (GRID_SPACING). In case you provide both PLUMED will use
106 : the most conservative choice (highest number of bins) for each dimension.
107 : In case you do not provide any information about bin size (neither GRID_BIN nor GRID_SPACING)
108 : and if Gaussian width is fixed PLUMED will use 1/5 of the Gaussian width as grid spacing.
109 : This default choice should be reasonable for most applications.
110 :
111 : Another option that is available is well-tempered metadynamics \cite Barducci:2008. In this
112 : variant of PBMetaD the heights of the Gaussian hills are scaled at each step by the
113 : additional well-tempered metadynamics term.
114 : This ensures that each bias converges more smoothly. It should be noted that, in the case of well-tempered metadynamics, in
115 : the output printed the Gaussian height is re-scaled using the bias factor.
116 : Also notice that with well-tempered metadynamics the HILLS files do not contain the bias,
117 : but the negative of the free-energy estimate. This choice has the advantage that
118 : one can restart a simulation using a different value for the \f$\Delta T\f$. The applied bias will be scaled accordingly.
119 :
120 : Note that you can use here also the flexible Gaussian approach \cite Branduardi:2012dl
121 : in which you can adapt the Gaussian to the extent of Cartesian space covered by a variable or
122 : to the space in collective variable covered in a given time. In this case the width of the deposited
123 : Gaussian potential is denoted by one value only that is a Cartesian space (ADAPTIVE=GEOM) or a time
124 : (ADAPTIVE=DIFF). Note that a specific integration technique for the deposited Gaussian kernels
125 : should be used in this case. Check the documentation for utility sum_hills.
126 :
127 : With the keyword INTERVAL one changes the metadynamics algorithm setting the bias force equal to zero
128 : outside boundary \cite baftizadeh2012protein. If, for example, metadynamics is performed on a CV s and one is interested only
129 : to the free energy for s > boundary, the history dependent potential is still updated according to the above
130 : equations but the metadynamics force is set to zero for s < boundary. Notice that Gaussian kernels are added also
131 : if s < boundary, as the tails of these Gaussian kernels influence VG in the relevant region s > boundary. In this way, the
132 : force on the system in the region s > boundary comes from both metadynamics and the force field, in the region
133 : s < boundary only from the latter. This approach allows obtaining a history-dependent bias potential VG that
134 : fluctuates around a stable estimator, equal to the negative of the free energy far enough from the
135 : boundaries. Note that:
136 : - It works only for one-dimensional biases;
137 : - It works both with and without GRID;
138 : - The interval limit boundary in a region where the free energy derivative is not large;
139 : - If in the region outside the limit boundary the system has a free energy minimum, the INTERVAL keyword should
140 : be used together with a \ref UPPER_WALLS or \ref LOWER_WALLS at boundary.
141 :
142 : Multiple walkers \cite multiplewalkers can also be used. See below the examples.
143 :
144 : \par Examples
145 :
146 : The following input is for PBMetaD calculation using as
147 : collective variables the distance between atoms 3 and 5
148 : and the distance between atoms 2 and 4. The value of the CVs and
149 : the PBMetaD bias potential are written to the COLVAR file every 100 steps.
150 : \plumedfile
151 : DISTANCE ATOMS=3,5 LABEL=d1
152 : DISTANCE ATOMS=2,4 LABEL=d2
153 : PBMETAD ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3 PACE=500 LABEL=pb FILE=HILLS_d1,HILLS_d2
154 : PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
155 : \endplumedfile
156 : (See also \ref DISTANCE and \ref PRINT).
157 :
158 : \par
159 : If you use well-tempered metadynamics, you should specify a single bias factor and initial
160 : Gaussian height.
161 : \plumedfile
162 : DISTANCE ATOMS=3,5 LABEL=d1
163 : DISTANCE ATOMS=2,4 LABEL=d2
164 : PBMETAD ...
165 : ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3
166 : PACE=500 BIASFACTOR=8 LABEL=pb
167 : FILE=HILLS_d1,HILLS_d2
168 : ... PBMETAD
169 : PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
170 : \endplumedfile
171 :
172 : \par
173 : The following input enables the MPI version of multiple-walkers.
174 : \plumedfile
175 : DISTANCE ATOMS=3,5 LABEL=d1
176 : DISTANCE ATOMS=2,4 LABEL=d2
177 : PBMETAD ...
178 : ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3
179 : PACE=500 BIASFACTOR=8 LABEL=pb
180 : FILE=HILLS_d1,HILLS_d2
181 : WALKERS_MPI
182 : ... PBMETAD
183 : PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
184 : \endplumedfile
185 :
186 : \par
187 : The disk version of multiple-walkers can be
188 : enabled by setting the number of walker used, the id of the
189 : current walker which interprets the input file, the directory where the
190 : hills containing files resides, and the frequency to read the other walkers.
191 : Here is an example
192 : \plumedfile
193 : DISTANCE ATOMS=3,5 LABEL=d1
194 : DISTANCE ATOMS=2,4 LABEL=d2
195 : PBMETAD ...
196 : ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3
197 : PACE=500 BIASFACTOR=8 LABEL=pb
198 : FILE=HILLS_d1,HILLS_d2
199 : WALKERS_N=10
200 : WALKERS_ID=3
201 : WALKERS_DIR=../
202 : WALKERS_RSTRIDE=100
203 : ... PBMETAD
204 : PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
205 : \endplumedfile
206 : where WALKERS_N is the total number of walkers, WALKERS_ID is the
207 : id of the present walker (starting from 0 ) and the WALKERS_DIR is the directory
208 : where all the walkers are located. WALKERS_RSTRIDE is the number of step between
209 : one update and the other.
210 :
211 : */
212 : //+ENDPLUMEDOC
213 :
214 : class PBMetaD : public Bias {
215 :
216 : private:
217 : struct Gaussian {
218 : std::vector<double> center;
219 : std::vector<double> sigma;
220 : double height;
221 : bool multivariate; // this is required to discriminate the one dimensional case
222 : std::vector<double> invsigma;
223 1072 : Gaussian(const std::vector<double> & center,const std::vector<double> & sigma, double height, bool multivariate):
224 1072 : center(center),sigma(sigma),height(height),multivariate(multivariate),invsigma(sigma) {
225 : // to avoid troubles from zero element in flexible hills
226 2144 : for(unsigned i=0; i<invsigma.size(); ++i) if(std::abs(invsigma[i])>1.e-20) invsigma[i]=1.0/invsigma[i] ; else invsigma[i]=0.0;
227 1072 : }
228 : };
229 : // general setup
230 : double kbt_;
231 : int stride_;
232 : // well-tempered MetaD
233 : bool welltemp_;
234 : double biasf_;
235 : // output files format
236 : std::string fmt_;
237 : // first step?
238 : bool isFirstStep_;
239 : // Gaussian starting parameters
240 : double height0_;
241 : std::vector<double> sigma0_;
242 : std::vector<double> sigma0min_;
243 : std::vector<double> sigma0max_;
244 : // Gaussians
245 : std::vector<std::vector<Gaussian> > hills_;
246 : std::vector<FlexibleBin> flexbin_;
247 : int adaptive_;
248 : std::vector<std::unique_ptr<OFile>> hillsOfiles_;
249 : std::vector<std::unique_ptr<IFile>> ifiles_;
250 : std::vector<std::string> ifilesnames_;
251 : // Grids
252 : bool grid_;
253 : std::vector<std::unique_ptr<GridBase>> BiasGrids_;
254 : std::vector<std::unique_ptr<OFile>> gridfiles_;
255 : int wgridstride_;
256 : // multiple walkers
257 : int mw_n_;
258 : std::string mw_dir_;
259 : int mw_id_;
260 : int mw_rstride_;
261 : bool walkers_mpi_;
262 : size_t mpi_nw_;
263 : unsigned mpi_id_;
264 : std::vector<std::string> hillsfname_;
265 : // intervals
266 : std::vector<double> uppI_;
267 : std::vector<double> lowI_;
268 : std::vector<bool> doInt_;
269 : // variable for selector
270 : std::string selector_;
271 : bool do_select_;
272 : unsigned select_value_;
273 : unsigned current_value_;
274 :
275 : double stretchA=1.0;
276 : double stretchB=0.0;
277 :
278 : bool noStretchWarningDone=false;
279 :
280 0 : void noStretchWarning() {
281 0 : if(!noStretchWarningDone) {
282 0 : log<<"\nWARNING: you are using a HILLS file with Gaussian kernels, PLUMED 2.8 uses stretched Gaussians by default\n";
283 : }
284 0 : noStretchWarningDone=true;
285 0 : }
286 :
287 : void readGaussians(unsigned iarg, IFile*);
288 : void writeGaussian(unsigned iarg, const Gaussian&, OFile*);
289 : void addGaussian(unsigned iarg, const Gaussian&);
290 : double getBiasAndDerivatives(unsigned iarg, const std::vector<double>&, double* der=NULL);
291 : double evaluateGaussian(unsigned iarg, const std::vector<double>&, const Gaussian&,double* der=NULL);
292 : std::vector<unsigned> getGaussianSupport(unsigned iarg, const Gaussian&);
293 : bool scanOneHill(unsigned iarg, IFile *ifile, std::vector<Value> &v, std::vector<double> ¢er, std::vector<double> &sigma, double &height, bool &multivariate);
294 :
295 : public:
296 : explicit PBMetaD(const ActionOptions&);
297 : void calculate() override;
298 : void update() override;
299 : static void registerKeywords(Keywords& keys);
300 : bool checkNeedsGradients()const override;
301 : };
302 :
303 10495 : PLUMED_REGISTER_ACTION(PBMetaD,"PBMETAD")
304 :
305 39 : void PBMetaD::registerKeywords(Keywords& keys) {
306 39 : Bias::registerKeywords(keys);
307 39 : keys.use("ARG");
308 78 : keys.add("compulsory","SIGMA","the widths of the Gaussian hills");
309 78 : keys.add("compulsory","PACE","the frequency for hill addition, one for all biases");
310 78 : 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");
311 78 : keys.add("optional","HEIGHT","the height of the Gaussian hills, one for all biases. Compulsory unless TAU, TEMP and BIASFACTOR are given");
312 78 : keys.add("optional","FMT","specify format for HILLS files (useful for decrease the number of digits in regtests)");
313 78 : keys.add("optional","BIASFACTOR","use well tempered metadynamics with this bias factor, one for all biases. Please note you must also specify temp");
314 78 : keys.add("optional","TEMP","the system temperature - this is only needed if you are doing well-tempered metadynamics");
315 78 : keys.add("optional","TAU","in well tempered metadynamics, sets height to (\\f$k_B \\Delta T\\f$*pace*timestep)/tau");
316 78 : keys.add("optional","GRID_MIN","the lower bounds for the grid");
317 78 : keys.add("optional","GRID_MAX","the upper bounds for the grid");
318 78 : keys.add("optional","GRID_BIN","the number of bins for the grid");
319 78 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
320 78 : keys.addFlag("GRID_SPARSE",false,"use a sparse grid to store hills");
321 78 : keys.addFlag("GRID_NOSPLINE",false,"don't use spline interpolation with grids");
322 78 : keys.add("optional","GRID_WSTRIDE", "frequency for dumping the grid");
323 78 : keys.add("optional","GRID_WFILES", "dump grid for the bias, default names are used if GRID_WSTRIDE is used without GRID_WFILES.");
324 78 : keys.add("optional","GRID_RFILES", "read grid for the bias");
325 78 : 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");
326 78 : keys.add("optional","SIGMA_MAX","the upper bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
327 78 : keys.add("optional","SIGMA_MIN","the lower bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
328 78 : keys.add("optional","SELECTOR", "add forces and do update based on the value of SELECTOR");
329 78 : keys.add("optional","SELECTOR_ID", "value of SELECTOR");
330 78 : keys.add("optional","WALKERS_ID", "walker id");
331 78 : keys.add("optional","WALKERS_N", "number of walkers");
332 78 : keys.add("optional","WALKERS_DIR", "shared directory with the hills files from all the walkers");
333 78 : keys.add("optional","WALKERS_RSTRIDE","stride for reading hills files");
334 78 : keys.addFlag("WALKERS_MPI",false,"Switch on MPI version of multiple walkers - not compatible with WALKERS_* options other than WALKERS_DIR");
335 78 : keys.add("optional","INTERVAL_MIN","one dimensional lower limits, outside the limits the system will not feel the biasing force.");
336 78 : keys.add("optional","INTERVAL_MAX","one dimensional upper limits, outside the limits the system will not feel the biasing force.");
337 39 : keys.use("RESTART");
338 39 : keys.use("UPDATE_FROM");
339 39 : keys.use("UPDATE_UNTIL");
340 39 : }
341 :
342 38 : PBMetaD::PBMetaD(const ActionOptions& ao):
343 : PLUMED_BIAS_INIT(ao),
344 38 : kbt_(0.0),
345 38 : stride_(0),
346 38 : welltemp_(false),
347 38 : biasf_(1.0),
348 38 : isFirstStep_(true),
349 38 : height0_(std::numeric_limits<double>::max()),
350 38 : adaptive_(FlexibleBin::none),
351 38 : grid_(false),
352 38 : wgridstride_(0),
353 38 : mw_n_(1), mw_dir_(""), mw_id_(0), mw_rstride_(1),
354 38 : walkers_mpi_(false), mpi_nw_(0),
355 38 : do_select_(false)
356 : {
357 : // parse the flexible hills
358 : std::string adaptiveoption;
359 : adaptiveoption="NONE";
360 76 : parse("ADAPTIVE",adaptiveoption);
361 38 : if(adaptiveoption=="GEOM") {
362 0 : log.printf(" Uses Geometry-based hills width: sigma must be in distance units and only one sigma is needed\n");
363 0 : adaptive_=FlexibleBin::geometry;
364 38 : } else if(adaptiveoption=="DIFF") {
365 4 : log.printf(" Uses Diffusion-based hills width: sigma must be in time steps and only one sigma is needed\n");
366 4 : adaptive_=FlexibleBin::diffusion;
367 34 : } else if(adaptiveoption=="NONE") {
368 34 : adaptive_=FlexibleBin::none;
369 : } else {
370 0 : error("I do not know this type of adaptive scheme");
371 : }
372 :
373 38 : parse("FMT",fmt_);
374 :
375 : // parse the sigma
376 38 : parseVector("SIGMA",sigma0_);
377 38 : if(adaptive_==FlexibleBin::none) {
378 : // if you use normal sigma you need one sigma per argument
379 34 : if( sigma0_.size()!=getNumberOfArguments() ) error("number of arguments does not match number of SIGMA parameters");
380 : } else {
381 : // if you use flexible hills you need one sigma
382 4 : if(sigma0_.size()!=1) {
383 0 : error("If you choose ADAPTIVE you need only one sigma according to your choice of type (GEOM/DIFF)");
384 : }
385 : // if adaptive then the number must be an integer
386 4 : if(adaptive_==FlexibleBin::diffusion) {
387 4 : if(int(sigma0_[0])-sigma0_[0]>1.e-9 || int(sigma0_[0])-sigma0_[0] <-1.e-9 || int(sigma0_[0])<1 ) {
388 0 : error("In case of adaptive hills with diffusion, the sigma must be an integer which is the number of time steps\n");
389 : }
390 : }
391 : // here evtl parse the sigma min and max values
392 8 : parseVector("SIGMA_MIN",sigma0min_);
393 4 : if(sigma0min_.size()>0 && sigma0min_.size()!=getNumberOfArguments()) {
394 0 : error("the number of SIGMA_MIN values be the same of the number of the arguments");
395 4 : } else if(sigma0min_.size()==0) {
396 0 : sigma0min_.resize(getNumberOfArguments());
397 0 : for(unsigned i=0; i<getNumberOfArguments(); i++) {sigma0min_[i]=-1.;}
398 : }
399 :
400 8 : parseVector("SIGMA_MAX",sigma0max_);
401 4 : if(sigma0max_.size()>0 && sigma0max_.size()!=getNumberOfArguments()) {
402 0 : error("the number of SIGMA_MAX values be the same of the number of the arguments");
403 4 : } else if(sigma0max_.size()==0) {
404 4 : sigma0max_.resize(getNumberOfArguments());
405 12 : for(unsigned i=0; i<getNumberOfArguments(); i++) {sigma0max_[i]=-1.;}
406 : }
407 :
408 12 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
409 : std::vector<double> tmp_smin, tmp_smax;
410 8 : tmp_smin.resize(1,sigma0min_[i]);
411 8 : tmp_smax.resize(1,sigma0max_[i]);
412 16 : flexbin_.push_back(FlexibleBin(adaptive_,this,i,sigma0_[0],tmp_smin,tmp_smax));
413 : }
414 : }
415 :
416 : // note: HEIGHT is not compulsory, since one could use the TAU keyword, see below
417 38 : parse("HEIGHT",height0_);
418 38 : parse("PACE",stride_);
419 38 : if(stride_<=0) error("frequency for hill addition is nonsensical");
420 :
421 76 : parseVector("FILE",hillsfname_);
422 38 : if(hillsfname_.size()==0) {
423 18 : for(unsigned i=0; i<getNumberOfArguments(); i++) hillsfname_.push_back("HILLS."+getPntrToArgument(i)->getName());
424 : }
425 38 : if( hillsfname_.size()!=getNumberOfArguments() ) {
426 0 : error("number of FILE arguments does not match number of HILLS files");
427 : }
428 :
429 38 : parse("BIASFACTOR",biasf_);
430 38 : if( biasf_<1.0 ) error("well tempered bias factor is nonsensical");
431 38 : double temp=0.0;
432 38 : parse("TEMP",temp);
433 38 : if(temp>0.0) kbt_=plumed.getAtoms().getKBoltzmann()*temp;
434 0 : else kbt_=plumed.getAtoms().getKbT();
435 38 : if(biasf_>1.0) {
436 37 : if(kbt_==0.0) error("Unless the MD engine passes the temperature to plumed, with well-tempered metad you must specify it using TEMP");
437 37 : welltemp_=true;
438 : }
439 38 : double tau=0.0;
440 38 : parse("TAU",tau);
441 38 : if(tau==0.0) {
442 38 : if(height0_==std::numeric_limits<double>::max()) error("At least one between HEIGHT and TAU should be specified");
443 : // if tau is not set, we compute it here from the other input parameters
444 38 : if(welltemp_) tau=(kbt_*(biasf_-1.0))/height0_*getTimeStep()*stride_;
445 : } else {
446 0 : if(!welltemp_)error("TAU only makes sense in well-tempered metadynamics");
447 0 : if(height0_!=std::numeric_limits<double>::max()) error("At most one between HEIGHT and TAU should be specified");
448 0 : height0_=(kbt_*(biasf_-1.0))/tau*getTimeStep()*stride_;
449 : }
450 :
451 : // Multiple walkers
452 38 : parse("WALKERS_N",mw_n_);
453 38 : parse("WALKERS_ID",mw_id_);
454 38 : if(mw_n_<=mw_id_) error("walker ID should be a numerical value less than the total number of walkers");
455 38 : parse("WALKERS_DIR",mw_dir_);
456 38 : parse("WALKERS_RSTRIDE",mw_rstride_);
457 :
458 : // MPI version
459 38 : parseFlag("WALKERS_MPI",walkers_mpi_);
460 :
461 : // Grid file
462 76 : parse("GRID_WSTRIDE",wgridstride_);
463 : std::vector<std::string> gridfilenames_;
464 38 : parseVector("GRID_WFILES",gridfilenames_);
465 38 : if (wgridstride_ == 0 && gridfilenames_.size() > 0) {
466 0 : error("frequency with which to output grid not specified use GRID_WSTRIDE");
467 : }
468 38 : if(gridfilenames_.size()==0 && wgridstride_ > 0) {
469 12 : for(unsigned i=0; i<getNumberOfArguments(); i++) gridfilenames_.push_back("GRID."+getPntrToArgument(i)->getName());
470 : }
471 38 : if(gridfilenames_.size() > 0 && hillsfname_.size() > 0 && gridfilenames_.size() != hillsfname_.size())
472 0 : error("number of GRID_WFILES arguments does not match number of HILLS files");
473 :
474 : // Read grid
475 : std::vector<std::string> gridreadfilenames_;
476 76 : parseVector("GRID_RFILES",gridreadfilenames_);
477 :
478 : // Grid Stuff
479 38 : std::vector<std::string> gmin(getNumberOfArguments());
480 76 : parseVector("GRID_MIN",gmin);
481 38 : if(gmin.size()!=getNumberOfArguments() && gmin.size()!=0) error("not enough values for GRID_MIN");
482 38 : std::vector<std::string> gmax(getNumberOfArguments());
483 76 : parseVector("GRID_MAX",gmax);
484 38 : if(gmax.size()!=getNumberOfArguments() && gmax.size()!=0) error("not enough values for GRID_MAX");
485 38 : std::vector<unsigned> gbin(getNumberOfArguments());
486 : std::vector<double> gspacing;
487 76 : parseVector("GRID_BIN",gbin);
488 38 : if(gbin.size()!=getNumberOfArguments() && gbin.size()!=0) error("not enough values for GRID_BIN");
489 76 : parseVector("GRID_SPACING",gspacing);
490 38 : if(gspacing.size()!=getNumberOfArguments() && gspacing.size()!=0) error("not enough values for GRID_SPACING");
491 38 : if(gmin.size()!=gmax.size()) error("GRID_MAX and GRID_MIN should be either present or absent");
492 38 : if(gspacing.size()!=0 && gmin.size()==0) error("If GRID_SPACING is present also GRID_MIN and GRID_MAX should be present");
493 38 : if(gbin.size()!=0 && gmin.size()==0) error("If GRID_BIN is present also GRID_MIN and GRID_MAX should be present");
494 38 : if(gmin.size()!=0) {
495 6 : if(gbin.size()==0 && gspacing.size()==0) {
496 6 : if(adaptive_==FlexibleBin::none) {
497 2 : log<<" Binsize not specified, 1/5 of sigma will be be used\n";
498 2 : plumed_assert(sigma0_.size()==getNumberOfArguments());
499 2 : gspacing.resize(getNumberOfArguments());
500 6 : for(unsigned i=0; i<gspacing.size(); i++) gspacing[i]=0.2*sigma0_[i];
501 : } else {
502 : // with adaptive hills and grid a sigma min must be specified
503 12 : 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");
504 4 : log<<" Binsize not specified, 1/5 of sigma_min will be be used\n";
505 4 : gspacing.resize(getNumberOfArguments());
506 12 : for(unsigned i=0; i<gspacing.size(); i++) gspacing[i]=0.2*sigma0min_[i];
507 : }
508 0 : } else if(gspacing.size()!=0 && gbin.size()==0) {
509 0 : log<<" The number of bins will be estimated from GRID_SPACING\n";
510 0 : } else if(gspacing.size()!=0 && gbin.size()!=0) {
511 0 : log<<" You specified both GRID_BIN and GRID_SPACING\n";
512 0 : log<<" The more conservative (highest) number of bins will be used for each variable\n";
513 : }
514 6 : if(gbin.size()==0) gbin.assign(getNumberOfArguments(),1);
515 18 : if(gspacing.size()!=0) for(unsigned i=0; i<getNumberOfArguments(); i++) {
516 : double a,b;
517 12 : Tools::convert(gmin[i],a);
518 12 : Tools::convert(gmax[i],b);
519 12 : unsigned n=((b-a)/gspacing[i])+1;
520 12 : if(gbin[i]<n) gbin[i]=n;
521 : }
522 : }
523 38 : if(gbin.size()>0) grid_=true;
524 :
525 38 : bool sparsegrid=false;
526 38 : parseFlag("GRID_SPARSE",sparsegrid);
527 38 : bool nospline=false;
528 38 : parseFlag("GRID_NOSPLINE",nospline);
529 38 : bool spline=!nospline;
530 38 : if(!grid_&&gridfilenames_.size() > 0) error("To write a grid you need first to define it!");
531 38 : if(!grid_&&gridreadfilenames_.size() > 0) error("To read a grid you need first to define it!");
532 :
533 38 : doInt_.resize(getNumberOfArguments(),false);
534 : // Interval keyword
535 38 : parseVector("INTERVAL_MIN",lowI_);
536 76 : parseVector("INTERVAL_MAX",uppI_);
537 : // various checks
538 38 : if(lowI_.size()!=uppI_.size()) error("both a lower and an upper limits must be provided with INTERVAL");
539 38 : if(lowI_.size()!=0 && lowI_.size()!=getNumberOfArguments()) error("check number of argument of INTERVAL");
540 46 : for(unsigned i=0; i<lowI_.size(); ++i) {
541 8 : if(uppI_[i]<lowI_[i]) error("The Upper limit must be greater than the Lower limit!");
542 8 : if(getPntrToArgument(i)->isPeriodic()) warning("INTERVAL is not used for periodic variables");
543 : else doInt_[i]=true;
544 : }
545 :
546 : // parse selector stuff
547 76 : parse("SELECTOR", selector_);
548 38 : if(selector_.length()>0) {
549 0 : do_select_ = true;
550 0 : parse("SELECTOR_ID", select_value_);
551 : }
552 :
553 38 : checkRead();
554 :
555 38 : log.printf(" Gaussian width ");
556 38 : if (adaptive_==FlexibleBin::diffusion)log.printf(" (Note: The units of sigma are in timesteps) ");
557 38 : if (adaptive_==FlexibleBin::geometry)log.printf(" (Note: The units of sigma are in dist units) ");
558 110 : for(unsigned i=0; i<sigma0_.size(); ++i) log.printf(" %f",sigma0_[i]);
559 38 : log.printf(" Gaussian height %f\n",height0_);
560 38 : log.printf(" Gaussian deposition pace %d\n",stride_);
561 38 : log.printf(" Gaussian files ");
562 114 : for(unsigned i=0; i<hillsfname_.size(); ++i) log.printf("%s ",hillsfname_[i].c_str());
563 38 : log.printf("\n");
564 38 : if(welltemp_) {
565 37 : log.printf(" Well-Tempered Bias Factor %f\n",biasf_);
566 37 : log.printf(" Hills relaxation time (tau) %f\n",tau);
567 37 : log.printf(" KbT %f\n",kbt_);
568 : }
569 :
570 38 : if(do_select_) {
571 0 : log.printf(" Add forces and update bias based on the value of SELECTOR %s\n",selector_.c_str());
572 0 : log.printf(" Id of the SELECTOR for this action %u\n", select_value_);
573 : }
574 :
575 38 : if(mw_n_>1) {
576 0 : if(walkers_mpi_) error("MPI version of multiple walkers is not compatible with filesystem version of multiple walkers");
577 0 : log.printf(" %d multiple walkers active\n",mw_n_);
578 0 : log.printf(" walker id %d\n",mw_id_);
579 0 : log.printf(" reading stride %d\n",mw_rstride_);
580 0 : if(mw_dir_!="")log.printf(" directory with hills files %s\n",mw_dir_.c_str());
581 : } else {
582 38 : if(walkers_mpi_) {
583 34 : log.printf(" Multiple walkers active using MPI communnication\n");
584 34 : if(mw_dir_!="")log.printf(" directory with hills files %s\n",mw_dir_.c_str());
585 34 : if(comm.Get_rank()==0) {
586 : // Only root of group can communicate with other walkers
587 18 : mpi_nw_ = multi_sim_comm.Get_size();
588 18 : mpi_id_ = multi_sim_comm.Get_rank();
589 : }
590 : // Communicate to the other members of the same group
591 : // info abount number of walkers and walker index
592 34 : comm.Bcast(mpi_nw_,0);
593 34 : comm.Bcast(mpi_id_,0);
594 : }
595 : }
596 :
597 114 : for(unsigned i=0; i<doInt_.size(); i++) {
598 76 : if(doInt_[i]) log.printf(" Upper and Lower limits boundaries for the bias of CV %u are activated\n", i);
599 : }
600 38 : if(grid_) {
601 6 : log.printf(" Grid min");
602 18 : for(unsigned i=0; i<gmin.size(); ++i) log.printf(" %s",gmin[i].c_str() );
603 6 : log.printf("\n");
604 6 : log.printf(" Grid max");
605 18 : for(unsigned i=0; i<gmax.size(); ++i) log.printf(" %s",gmax[i].c_str() );
606 6 : log.printf("\n");
607 6 : log.printf(" Grid bin");
608 18 : for(unsigned i=0; i<gbin.size(); ++i) log.printf(" %u",gbin[i]);
609 6 : log.printf("\n");
610 6 : if(spline) {log.printf(" Grid uses spline interpolation\n");}
611 6 : if(sparsegrid) {log.printf(" Grid uses sparse grid\n");}
612 6 : if(wgridstride_>0) {
613 18 : for(unsigned i=0; i<gridfilenames_.size(); ++i) {
614 12 : log.printf(" Grid is written on file %s with stride %d\n",gridfilenames_[i].c_str(),wgridstride_);
615 : }
616 : }
617 6 : if(gridreadfilenames_.size()>0) {
618 3 : for(unsigned i=0; i<gridreadfilenames_.size(); ++i) {
619 2 : log.printf(" Reading bias from grid in file %s \n",gridreadfilenames_[i].c_str());
620 : }
621 : }
622 : }
623 :
624 : // initializing vector of hills
625 38 : hills_.resize(getNumberOfArguments());
626 :
627 : // restart from external grid
628 : bool restartedFromGrid=false;
629 :
630 : // initializing and checking grid
631 38 : if(grid_) {
632 : // check for mesh and sigma size
633 18 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
634 : double a,b;
635 12 : Tools::convert(gmin[i],a);
636 12 : Tools::convert(gmax[i],b);
637 12 : double mesh=(b-a)/((double)gbin[i]);
638 12 : if(adaptive_==FlexibleBin::none) {
639 4 : if(mesh>0.5*sigma0_[i]) log<<" WARNING: Using a PBMETAD with a Grid Spacing larger than half of the Gaussians width can produce artifacts\n";
640 : } else {
641 8 : if(mesh>0.5*sigma0min_[i]||sigma0min_[i]<0.) log<<" WARNING: to use a PBMETAD with a GRID and ADAPTIVE you need to set a Grid Spacing larger than half of the Gaussians \n";
642 : }
643 : }
644 6 : std::string funcl=getLabel() + ".bias";
645 18 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
646 12 : std::vector<Value*> args(1);
647 12 : args[0] = getPntrToArgument(i);
648 12 : std::vector<std::string> gmin_t(1);
649 12 : std::vector<std::string> gmax_t(1);
650 12 : std::vector<unsigned> gbin_t(1);
651 : gmin_t[0] = gmin[i];
652 : gmax_t[0] = gmax[i];
653 12 : gbin_t[0] = gbin[i];
654 12 : std::unique_ptr<GridBase> BiasGrid_;
655 : // Read grid from file
656 12 : if(gridreadfilenames_.size()>0) {
657 2 : IFile gridfile;
658 2 : gridfile.link(*this);
659 2 : if(gridfile.FileExist(gridreadfilenames_[i])) {
660 2 : gridfile.open(gridreadfilenames_[i]);
661 : } else {
662 0 : error("The GRID file you want to read: " + gridreadfilenames_[i] + ", cannot be found!");
663 : }
664 2 : std::string funcl = getLabel() + ".bias";
665 4 : BiasGrid_=GridBase::create(funcl, args, gridfile, gmin_t, gmax_t, gbin_t, sparsegrid, spline, true);
666 2 : if(BiasGrid_->getDimension() != args.size()) {
667 0 : error("mismatch between dimensionality of input grid and number of arguments");
668 : }
669 4 : if(getPntrToArgument(i)->isPeriodic() != BiasGrid_->getIsPeriodic()[0]) {
670 0 : error("periodicity mismatch between arguments and input bias");
671 : }
672 2 : log.printf(" Restarting from %s:\n",gridreadfilenames_[i].c_str());
673 2 : if(getRestart()) restartedFromGrid=true;
674 2 : } else {
675 10 : if(!sparsegrid) {BiasGrid_=Tools::make_unique<Grid>(funcl,args,gmin_t,gmax_t,gbin_t,spline,true);}
676 0 : else {BiasGrid_=Tools::make_unique<SparseGrid>(funcl,args,gmin_t,gmax_t,gbin_t,spline,true);}
677 10 : std::vector<std::string> actualmin=BiasGrid_->getMin();
678 10 : std::vector<std::string> actualmax=BiasGrid_->getMax();
679 : std::string is;
680 10 : Tools::convert(i,is);
681 10 : if(gmin_t[0]!=actualmin[0]) error("GRID_MIN["+is+"] must be adjusted to "+actualmin[0]+" to fit periodicity");
682 10 : if(gmax_t[0]!=actualmax[0]) error("GRID_MAX["+is+"] must be adjusted to "+actualmax[0]+" to fit periodicity");
683 10 : }
684 12 : BiasGrids_.emplace_back(std::move(BiasGrid_));
685 24 : }
686 : }
687 :
688 :
689 : // creating vector of ifile* for hills reading
690 : // open all files at the beginning and read Gaussians if restarting
691 76 : for(int j=0; j<mw_n_; ++j) {
692 114 : for(unsigned i=0; i<hillsfname_.size(); ++i) {
693 76 : unsigned k=j*hillsfname_.size()+i;
694 : std::string fname;
695 76 : if(mw_dir_!="") {
696 0 : if(mw_n_>1) {
697 0 : std::stringstream out; out << j;
698 0 : fname = mw_dir_+"/"+hillsfname_[i]+"."+out.str();
699 0 : } else if(walkers_mpi_) {
700 0 : fname = mw_dir_+"/"+hillsfname_[i];
701 : } else {
702 : fname = hillsfname_[i];
703 : }
704 : } else {
705 76 : if(mw_n_>1) {
706 0 : std::stringstream out; out << j;
707 0 : fname = hillsfname_[i]+"."+out.str();
708 0 : } else {
709 : fname = hillsfname_[i];
710 : }
711 : }
712 76 : ifiles_.emplace_back(Tools::make_unique<IFile>());
713 : // this is just a shortcut pointer to the last element:
714 : IFile *ifile = ifiles_.back().get();
715 76 : ifile->link(*this);
716 76 : ifilesnames_.push_back(fname);
717 76 : if(ifile->FileExist(fname)) {
718 18 : ifile->open(fname);
719 18 : if(getRestart()&&!restartedFromGrid) {
720 2 : log.printf(" Restarting from %s:",ifilesnames_[k].c_str());
721 2 : readGaussians(i,ifiles_[k].get());
722 : }
723 18 : ifiles_[k]->reset(false);
724 : // close only the walker own hills file for later writing
725 18 : if(j==mw_id_) ifiles_[k]->close();
726 : } else {
727 : // in case a file does not exist and we are restarting, complain that the file was not found
728 58 : if(getRestart()) log<<" WARNING: restart file "<<fname<<" not found\n";
729 : }
730 : }
731 : }
732 :
733 38 : comm.Barrier();
734 38 : if(comm.Get_rank()==0 && walkers_mpi_) multi_sim_comm.Barrier();
735 :
736 : // open hills files for writing
737 114 : for(unsigned i=0; i<hillsfname_.size(); ++i) {
738 76 : auto ofile=Tools::make_unique<OFile>();
739 76 : ofile->link(*this);
740 : // if MPI multiple walkers, only rank 0 will write to file
741 76 : if(walkers_mpi_) {
742 68 : int r=0;
743 68 : if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
744 68 : comm.Bcast(r,0);
745 68 : if(r>0) ifilesnames_[mw_id_*hillsfname_.size()+i]="/dev/null";
746 136 : ofile->enforceSuffix("");
747 : }
748 76 : if(mw_n_>1) ofile->enforceSuffix("");
749 76 : ofile->open(ifilesnames_[mw_id_*hillsfname_.size()+i]);
750 76 : if(fmt_.length()>0) ofile->fmtField(fmt_);
751 152 : ofile->addConstantField("multivariate");
752 152 : ofile->addConstantField("kerneltype");
753 76 : if(doInt_[i]) {
754 16 : ofile->addConstantField("lower_int").printField("lower_int",lowI_[i]);
755 16 : ofile->addConstantField("upper_int").printField("upper_int",uppI_[i]);
756 : }
757 76 : ofile->setHeavyFlush();
758 : // output periodicities of variables
759 76 : ofile->setupPrintValue( getPntrToArgument(i) );
760 : // push back
761 76 : hillsOfiles_.emplace_back(std::move(ofile));
762 76 : }
763 :
764 : // Dump grid to files
765 38 : if(wgridstride_ > 0) {
766 18 : for(unsigned i = 0; i < gridfilenames_.size(); ++i) {
767 12 : auto ofile=Tools::make_unique<OFile>();
768 12 : ofile->link(*this);
769 12 : std::string gridfname_tmp = gridfilenames_[i];
770 12 : if(walkers_mpi_) {
771 8 : int r = 0;
772 8 : if(comm.Get_rank() == 0) {
773 4 : r = multi_sim_comm.Get_rank();
774 : }
775 8 : comm.Bcast(r, 0);
776 8 : if(r>0) {
777 : gridfname_tmp = "/dev/null";
778 : }
779 16 : ofile->enforceSuffix("");
780 : }
781 12 : if(mw_n_>1) ofile->enforceSuffix("");
782 12 : ofile->open(gridfname_tmp);
783 12 : ofile->setHeavyFlush();
784 12 : gridfiles_.emplace_back(std::move(ofile));
785 12 : }
786 : }
787 :
788 114 : log<<" Bibliography "<<plumed.cite("Pfaendtner and Bonomi. J. Chem. Theory Comput. 11, 5062 (2015)");
789 50 : if(doInt_[0]) log<<plumed.cite(
790 8 : "Baftizadeh, Cossio, Pietrucci, and Laio, Curr. Phys. Chem. 2, 79 (2012)");
791 140 : if(mw_n_>1||walkers_mpi_) log<<plumed.cite(
792 68 : "Raiteri, Laio, Gervasio, Micheletti, and Parrinello, J. Phys. Chem. B 110, 3533 (2006)");
793 50 : if(adaptive_!=FlexibleBin::none) log<<plumed.cite(
794 8 : "Branduardi, Bussi, and Parrinello, J. Chem. Theory Comput. 8, 2247 (2012)");
795 38 : log<<"\n";
796 76 : }
797 :
798 2 : void PBMetaD::readGaussians(unsigned iarg, IFile *ifile)
799 : {
800 2 : std::vector<double> center(1);
801 2 : std::vector<double> sigma(1);
802 : double height;
803 : int nhills=0;
804 2 : bool multivariate=false;
805 :
806 : std::vector<Value> tmpvalues;
807 4 : tmpvalues.push_back( Value( this, getPntrToArgument(iarg)->getName(), false ) );
808 :
809 10 : while(scanOneHill(iarg,ifile,tmpvalues,center,sigma,height,multivariate)) {
810 : ;
811 8 : nhills++;
812 8 : if(welltemp_) {height*=(biasf_-1.0)/biasf_;}
813 8 : addGaussian(iarg, Gaussian(center,sigma,height,multivariate));
814 : }
815 2 : log.printf(" %d Gaussians read\n",nhills);
816 4 : }
817 :
818 1064 : void PBMetaD::writeGaussian(unsigned iarg, const Gaussian& hill, OFile *ofile)
819 : {
820 2128 : ofile->printField("time",getTimeStep()*getStep());
821 1064 : ofile->printField(getPntrToArgument(iarg),hill.center[0]);
822 :
823 2128 : ofile->printField("kerneltype","stretched-gaussian");
824 1064 : if(hill.multivariate) {
825 288 : ofile->printField("multivariate","true");
826 144 : double lower = std::sqrt(1./hill.sigma[0]);
827 288 : ofile->printField("sigma_"+getPntrToArgument(iarg)->getName()+"_"+
828 : getPntrToArgument(iarg)->getName(),lower);
829 : } else {
830 1840 : ofile->printField("multivariate","false");
831 1840 : ofile->printField("sigma_"+getPntrToArgument(iarg)->getName(),hill.sigma[0]);
832 : }
833 1064 : double height=hill.height;
834 1064 : if(welltemp_) height *= biasf_/(biasf_-1.0);
835 1064 : ofile->printField("height",height);
836 1064 : ofile->printField("biasf",biasf_);
837 1064 : if(mw_n_>1) ofile->printField("clock",int(std::time(0)));
838 1064 : ofile->printField();
839 1064 : }
840 :
841 1072 : void PBMetaD::addGaussian(unsigned iarg, const Gaussian& hill)
842 : {
843 1072 : if(!grid_) {hills_[iarg].push_back(hill);}
844 : else {
845 160 : std::vector<unsigned> nneighb=getGaussianSupport(iarg, hill);
846 160 : std::vector<Grid::index_t> neighbors=BiasGrids_[iarg]->getNeighbors(hill.center,nneighb);
847 160 : std::vector<double> der(1);
848 160 : std::vector<double> xx(1);
849 160 : if(comm.Get_size()==1) {
850 608 : for(unsigned i=0; i<neighbors.size(); ++i) {
851 592 : Grid::index_t ineigh=neighbors[i];
852 592 : der[0]=0.0;
853 592 : BiasGrids_[iarg]->getPoint(ineigh,xx);
854 592 : double bias=evaluateGaussian(iarg,xx,hill,&der[0]);
855 592 : BiasGrids_[iarg]->addValueAndDerivatives(ineigh,bias,der);
856 : }
857 : } else {
858 144 : unsigned stride=comm.Get_size();
859 144 : unsigned rank=comm.Get_rank();
860 144 : std::vector<double> allder(neighbors.size(),0.0);
861 144 : std::vector<double> allbias(neighbors.size(),0.0);
862 2808 : for(unsigned i=rank; i<neighbors.size(); i+=stride) {
863 2664 : Grid::index_t ineigh=neighbors[i];
864 2664 : BiasGrids_[iarg]->getPoint(ineigh,xx);
865 2664 : allbias[i]=evaluateGaussian(iarg,xx,hill,&allder[i]);
866 : }
867 144 : comm.Sum(allbias);
868 144 : comm.Sum(allder);
869 5472 : for(unsigned i=0; i<neighbors.size(); ++i) {
870 5328 : Grid::index_t ineigh=neighbors[i];
871 5328 : der[0]=allder[i];
872 5328 : BiasGrids_[iarg]->addValueAndDerivatives(ineigh,allbias[i],der);
873 : }
874 : }
875 : }
876 1072 : }
877 :
878 160 : std::vector<unsigned> PBMetaD::getGaussianSupport(unsigned iarg, const Gaussian& hill)
879 : {
880 : std::vector<unsigned> nneigh;
881 : double cutoff;
882 160 : if(hill.multivariate) {
883 144 : double maxautoval=1./hill.sigma[0];
884 144 : cutoff=std::sqrt(2.0*dp2cutoff*maxautoval);
885 : } else {
886 16 : cutoff=std::sqrt(2.0*dp2cutoff)*hill.sigma[0];
887 : }
888 :
889 160 : if(doInt_[iarg]) {
890 144 : if(hill.center[0]+cutoff > uppI_[iarg] || hill.center[0]-cutoff < lowI_[iarg]) {
891 : // in this case, we updated the entire grid to avoid problems
892 0 : return BiasGrids_[iarg]->getNbin();
893 : } else {
894 288 : nneigh.push_back( static_cast<unsigned>(ceil(cutoff/BiasGrids_[iarg]->getDx()[0])));
895 : return nneigh;
896 : }
897 : }
898 :
899 32 : nneigh.push_back( static_cast<unsigned>(ceil(cutoff/BiasGrids_[iarg]->getDx()[0])) );
900 :
901 : return nneigh;
902 : }
903 :
904 1188 : double PBMetaD::getBiasAndDerivatives(unsigned iarg, const std::vector<double>& cv, double* der)
905 : {
906 1188 : double bias=0.0;
907 1188 : if(!grid_) {
908 1000 : unsigned stride=comm.Get_size();
909 1000 : unsigned rank=comm.Get_rank();
910 5520 : for(unsigned i=rank; i<hills_[iarg].size(); i+=stride) {
911 4520 : bias += evaluateGaussian(iarg,cv,hills_[iarg][i],der);
912 : }
913 1000 : comm.Sum(bias);
914 1000 : if(der) comm.Sum(der,1);
915 : } else {
916 188 : if(der) {
917 100 : std::vector<double> vder(1);
918 100 : bias = BiasGrids_[iarg]->getValueAndDerivatives(cv,vder);
919 100 : der[0] = vder[0];
920 : } else {
921 88 : bias = BiasGrids_[iarg]->getValue(cv);
922 : }
923 : }
924 :
925 1188 : return bias;
926 : }
927 :
928 7776 : double PBMetaD::evaluateGaussian(unsigned iarg, const std::vector<double>& cv, const Gaussian& hill, double* der)
929 : {
930 : double bias=0.0;
931 : // I use a pointer here because cv is const (and should be const)
932 : // but when using doInt it is easier to locally replace cv[0] with
933 : // the upper/lower limit in case it is out of range
934 : const double *pcv=NULL;
935 : double tmpcv[1]; // tmp array with cv (to be used with doInt_)
936 7776 : tmpcv[0]=cv[0];
937 : bool isOutOfInt = false;
938 7776 : if(doInt_[iarg]) {
939 2664 : if(cv[0]<lowI_[iarg]) { tmpcv[0]=lowI_[iarg]; isOutOfInt = true; }
940 2664 : else if(cv[0]>uppI_[iarg]) { tmpcv[0]=uppI_[iarg]; isOutOfInt = true; }
941 : }
942 : pcv=&(tmpcv[0]);
943 :
944 7776 : if(hill.multivariate) {
945 2664 : double dp = difference(iarg, hill.center[0], pcv[0]);
946 2664 : double dp2 = 0.5 * dp * dp * hill.sigma[0];
947 2664 : if(dp2<dp2cutoff) {
948 2534 : bias = hill.height*std::exp(-dp2);
949 2534 : if(der && !isOutOfInt) {
950 2534 : der[0] += -bias * dp * hill.sigma[0] * stretchA;
951 : }
952 2534 : bias=stretchA*bias+hill.height*stretchB;
953 : }
954 : } else {
955 5112 : double dp = difference(iarg, hill.center[0], pcv[0]) * hill.invsigma[0];
956 5112 : double dp2 = 0.5 * dp * dp;
957 5112 : if(dp2<dp2cutoff) {
958 5088 : bias = hill.height*std::exp(-dp2);
959 5088 : if(der && !isOutOfInt) {
960 2832 : der[0] += -bias * dp * hill.invsigma[0] * stretchA;
961 : }
962 5088 : bias=stretchA*bias+hill.height*stretchB;
963 : }
964 : }
965 :
966 7776 : return bias;
967 : }
968 :
969 320 : void PBMetaD::calculate()
970 : {
971 : // this is because presently there is no way to properly pass information
972 : // on adaptive hills (diff) after exchanges:
973 320 : if(adaptive_==FlexibleBin::diffusion && getExchangeStep()) error("ADAPTIVE=DIFF is not compatible with replica exchange");
974 :
975 320 : std::vector<double> cv(1);
976 : double der[1];
977 320 : std::vector<double> bias(getNumberOfArguments());
978 320 : std::vector<double> deriv(getNumberOfArguments());
979 :
980 320 : double ncv = (double) getNumberOfArguments();
981 : double bmin = 1.0e+19;
982 960 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
983 640 : cv[0] = getArgument(i);
984 640 : der[0] = 0.0;
985 640 : bias[i] = getBiasAndDerivatives(i, cv, der);
986 640 : deriv[i] = der[0];
987 640 : if(bias[i] < bmin) bmin = bias[i];
988 : }
989 : double ene = 0.;
990 960 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
991 640 : ene += std::exp((-bias[i]+bmin)/kbt_);
992 : }
993 :
994 : // set Forces - set them to zero if SELECTOR is active
995 320 : if(do_select_) current_value_ = static_cast<unsigned>(plumed.passMap[selector_]);
996 :
997 320 : if(!do_select_ || (do_select_ && select_value_==current_value_)) {
998 960 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
999 640 : const double f = - std::exp((-bias[i]+bmin)/kbt_) / (ene) * deriv[i];
1000 640 : setOutputForce(i, f);
1001 : }
1002 : }
1003 :
1004 320 : if(do_select_ && select_value_!=current_value_) {
1005 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) setOutputForce(i, 0.0);
1006 : }
1007 :
1008 : // set bias
1009 320 : ene = -kbt_ * (std::log(ene) - std::log(ncv)) + bmin;
1010 : setBias(ene);
1011 320 : }
1012 :
1013 320 : void PBMetaD::update()
1014 : {
1015 : bool multivariate;
1016 : // adding hills criteria
1017 : bool nowAddAHill;
1018 320 : if(getStep()%stride_==0 && !isFirstStep_) nowAddAHill=true;
1019 : else {
1020 : nowAddAHill=false;
1021 46 : isFirstStep_=false;
1022 : }
1023 :
1024 : // if you use adaptive, call the FlexibleBin
1025 320 : if(adaptive_!=FlexibleBin::none) {
1026 120 : for(unsigned i=0; i<getNumberOfArguments(); i++) flexbin_[i].update(nowAddAHill,i);
1027 : multivariate=true;
1028 : } else {
1029 : multivariate=false;
1030 : }
1031 :
1032 320 : if(nowAddAHill && (!do_select_ || (do_select_ && select_value_==current_value_))) {
1033 : // get all biases and heights
1034 274 : std::vector<double> cv(getNumberOfArguments());
1035 274 : std::vector<double> bias(getNumberOfArguments());
1036 274 : std::vector<double> thissigma(getNumberOfArguments());
1037 274 : std::vector<double> height(getNumberOfArguments());
1038 274 : std::vector<double> cv_tmp(1);
1039 274 : std::vector<double> sigma_tmp(1);
1040 : double norm = 0.0;
1041 : double bmin = 1.0e+19;
1042 822 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1043 620 : if(adaptive_!=FlexibleBin::none) thissigma[i]=flexbin_[i].getInverseMatrix(i)[0];
1044 476 : else thissigma[i]=sigma0_[i];
1045 548 : cv[i] = getArgument(i);
1046 548 : cv_tmp[0] = getArgument(i);
1047 548 : bias[i] = getBiasAndDerivatives(i, cv_tmp);
1048 548 : if(bias[i] < bmin) bmin = bias[i];
1049 : }
1050 : // calculate heights and norm
1051 822 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1052 548 : double h = std::exp((-bias[i]+bmin)/kbt_);
1053 548 : norm += h;
1054 548 : height[i] = h;
1055 : }
1056 : // normalize and apply welltemp correction
1057 822 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1058 548 : height[i] *= height0_ / norm;
1059 548 : if(welltemp_) height[i] *= std::exp(-bias[i]/(kbt_*(biasf_-1.0)));
1060 : }
1061 :
1062 : // MPI Multiple walkers: share hills and add them all
1063 274 : if(walkers_mpi_) {
1064 : // Allocate arrays to store all walkers hills
1065 258 : std::vector<double> all_cv(mpi_nw_*cv.size(), 0.0);
1066 258 : std::vector<double> all_sigma(mpi_nw_*getNumberOfArguments(), 0.0);
1067 258 : std::vector<double> all_height(mpi_nw_*height.size(), 0.0);
1068 258 : if(comm.Get_rank()==0) {
1069 : // fill in value
1070 390 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1071 260 : unsigned j = mpi_id_ * getNumberOfArguments() + i;
1072 260 : all_cv[j] = cv[i];
1073 260 : all_sigma[j] = thissigma[i];
1074 260 : all_height[j] = height[i];
1075 : }
1076 : // Communicate (only root)
1077 130 : multi_sim_comm.Sum(&all_cv[0], all_cv.size());
1078 130 : multi_sim_comm.Sum(&all_sigma[0], all_sigma.size());
1079 130 : multi_sim_comm.Sum(&all_height[0], all_height.size());
1080 : }
1081 : // Share info with group members
1082 258 : comm.Sum(&all_cv[0], all_cv.size());
1083 258 : comm.Sum(&all_sigma[0], all_sigma.size());
1084 258 : comm.Sum(&all_height[0], all_height.size());
1085 : // now add hills one by one
1086 774 : for(unsigned j=0; j<mpi_nw_; ++j) {
1087 1548 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1088 1032 : cv_tmp[0] = all_cv[j*cv.size()+i];
1089 1032 : double height_tmp = all_height[j*cv.size()+i];
1090 1032 : sigma_tmp[0] = all_sigma[j*cv.size()+i];
1091 1032 : Gaussian newhill = Gaussian(cv_tmp, sigma_tmp, height_tmp, multivariate);
1092 1032 : addGaussian(i, newhill);
1093 1032 : writeGaussian(i, newhill, hillsOfiles_[i].get());
1094 1032 : }
1095 : }
1096 : // just add your own hills
1097 : } else {
1098 48 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1099 32 : cv_tmp[0] = cv[i];
1100 32 : if(adaptive_!=FlexibleBin::none) sigma_tmp[0]=thissigma[i];
1101 32 : else sigma_tmp[0] = sigma0_[i];
1102 32 : Gaussian newhill = Gaussian(cv_tmp, sigma_tmp, height[i], multivariate);
1103 32 : addGaussian(i, newhill);
1104 32 : writeGaussian(i, newhill, hillsOfiles_[i].get());
1105 32 : }
1106 : }
1107 : }
1108 :
1109 : // write grid files
1110 320 : if(wgridstride_>0 && (getStep()%wgridstride_==0 || getCPT())) {
1111 14 : int r = 0;
1112 14 : if(walkers_mpi_) {
1113 4 : if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
1114 4 : comm.Bcast(r,0);
1115 : }
1116 14 : if(r==0) {
1117 36 : for(unsigned i=0; i<gridfiles_.size(); ++i) {
1118 24 : gridfiles_[i]->rewind();
1119 24 : BiasGrids_[i]->writeToFile(*gridfiles_[i]);
1120 24 : gridfiles_[i]->flush();
1121 : }
1122 : }
1123 : }
1124 :
1125 : // if multiple walkers and time to read Gaussians
1126 320 : if(mw_n_>1 && getStep()%mw_rstride_==0) {
1127 0 : for(int j=0; j<mw_n_; ++j) {
1128 0 : for(unsigned i=0; i<hillsfname_.size(); ++i) {
1129 0 : unsigned k=j*hillsfname_.size()+i;
1130 : // don't read your own Gaussians
1131 0 : if(j==mw_id_) continue;
1132 : // if the file is not open yet
1133 0 : if(!(ifiles_[k]->isOpen())) {
1134 : // check if it exists now and open it!
1135 0 : if(ifiles_[k]->FileExist(ifilesnames_[k])) {
1136 0 : ifiles_[k]->open(ifilesnames_[k]);
1137 0 : ifiles_[k]->reset(false);
1138 : }
1139 : // otherwise read the new Gaussians
1140 : } else {
1141 0 : log.printf(" Reading hills from %s:",ifilesnames_[k].c_str());
1142 0 : readGaussians(i,ifiles_[k].get());
1143 0 : ifiles_[k]->reset(false);
1144 : }
1145 : }
1146 : }
1147 : }
1148 :
1149 320 : }
1150 :
1151 : /// takes a pointer to the file and a template string with values v and gives back the next center, sigma and height
1152 10 : bool PBMetaD::scanOneHill(unsigned iarg, IFile *ifile, std::vector<Value> &tmpvalues, std::vector<double> ¢er, std::vector<double> &sigma, double &height, bool &multivariate)
1153 : {
1154 : double dummy;
1155 10 : multivariate=false;
1156 20 : if(ifile->scanField("time",dummy)) {
1157 8 : ifile->scanField( &tmpvalues[0] );
1158 8 : if( tmpvalues[0].isPeriodic() && ! getPntrToArgument(iarg)->isPeriodic() ) {
1159 0 : error("in hills file periodicity for variable " + tmpvalues[0].getName() + " does not match periodicity in input");
1160 8 : } else if( tmpvalues[0].isPeriodic() ) {
1161 0 : std::string imin, imax; tmpvalues[0].getDomain( imin, imax );
1162 0 : std::string rmin, rmax; getPntrToArgument(iarg)->getDomain( rmin, rmax );
1163 0 : if( imin!=rmin || imax!=rmax ) {
1164 0 : error("in hills file periodicity for variable " + tmpvalues[0].getName() + " does not match periodicity in input");
1165 : }
1166 : }
1167 8 : center[0]=tmpvalues[0].get();
1168 8 : std::string ktype="stretched-gaussian";
1169 24 : if( ifile->FieldExist("kerneltype") ) ifile->scanField("kerneltype",ktype);
1170 :
1171 8 : if( ktype=="gaussian" ) {
1172 0 : noStretchWarning();
1173 8 : } else if( ktype!="stretched-gaussian") {
1174 0 : error("non Gaussian kernels are not supported in MetaD");
1175 : }
1176 :
1177 : std::string sss;
1178 16 : ifile->scanField("multivariate",sss);
1179 8 : if(sss=="true") multivariate=true;
1180 8 : else if(sss=="false") multivariate=false;
1181 0 : else plumed_merror("cannot parse multivariate = "+ sss);
1182 8 : if(multivariate) {
1183 0 : ifile->scanField("sigma_"+getPntrToArgument(iarg)->getName()+"_"+
1184 : getPntrToArgument(iarg)->getName(),sigma[0]);
1185 0 : sigma[0] = 1./(sigma[0]*sigma[0]);
1186 : } else {
1187 16 : ifile->scanField("sigma_"+getPntrToArgument(iarg)->getName(),sigma[0]);
1188 : }
1189 8 : ifile->scanField("height",height);
1190 8 : ifile->scanField("biasf",dummy);
1191 16 : if(ifile->FieldExist("clock")) ifile->scanField("clock",dummy);
1192 16 : if(ifile->FieldExist("lower_int")) ifile->scanField("lower_int",dummy);
1193 16 : if(ifile->FieldExist("upper_int")) ifile->scanField("upper_int",dummy);
1194 8 : ifile->scanField();
1195 : return true;
1196 : } else {
1197 : return false;
1198 : }
1199 :
1200 : }
1201 :
1202 320 : bool PBMetaD::checkNeedsGradients()const
1203 : {
1204 320 : if(adaptive_==FlexibleBin::geometry) {
1205 0 : if(getStep()%stride_==0 && !isFirstStep_) return true;
1206 0 : else return false;
1207 : } else return false;
1208 : }
1209 :
1210 : }
1211 : }
1212 :
|
