Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2020-2021 of Michele Invernizzi.
3 :
4 : This file is part of the OPES plumed module.
5 :
6 : The OPES plumed module is free software: you can redistribute it and/or modify
7 : it under the terms of the GNU Lesser General Public License as published by
8 : the Free Software Foundation, either version 3 of the License, or
9 : (at your option) any later version.
10 :
11 : The OPES plumed module is distributed in the hope that it will be useful,
12 : but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : GNU Lesser General Public License for more details.
15 :
16 : You should have received a copy of the GNU Lesser General Public License
17 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
18 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
19 : #include "bias/Bias.h"
20 : #include "core/PlumedMain.h"
21 : #include "core/ActionRegister.h"
22 : #include "tools/Communicator.h"
23 : #include "tools/File.h"
24 : #include "tools/OpenMP.h"
25 :
26 : namespace PLMD {
27 : namespace opes {
28 :
29 : //+PLUMEDOC OPES_BIAS OPES_METAD
30 : /*
31 : On-the-fly probability enhanced sampling with metadynamics-like target distribution.
32 :
33 : This On-the-fly probability enhanced sampling (\ref OPES "OPES") method with metadynamics-like target distribution is described in \cite Invernizzi2020rethinking.
34 :
35 : This \ref OPES_METAD action samples target distributions defined via their marginal \f$p^{\text{tg}}(\mathbf{s})\f$ over some collective variables (CVs), \f$\mathbf{s}=\mathbf{s}(\mathbf{x})\f$.
36 : By default \ref OPES_METAD targets the well-tempered distribution, \f$p^{\text{WT}}(\mathbf{s})\propto [P(\mathbf{s})]^{1/\gamma}\f$, where \f$\gamma\f$ is known as BIASFACTOR.
37 : Similarly to \ref METAD, \ref OPES_METAD optimizes the bias on-the-fly, with a given PACE.
38 : It does so by reweighting via kernel density estimation the unbiased distribution in the CV space, \f$P(\mathbf{s})\f$.
39 : A compression algorithm is used to prevent the number of kernels from growing linearly with the simulation time.
40 : The bias at step \f$n\f$ is
41 : \f[
42 : V_n(\mathbf{s}) = (1-1/\gamma)\frac{1}{\beta}\log\left(\frac{P_n(\mathbf{s})}{Z_n}+\epsilon\right)\, .
43 : \f]
44 : See Ref.\cite Invernizzi2020rethinking for a complete description of the method.
45 :
46 : As an intuitive picture, rather than gradually filling the metastable basins, \ref OPES_METAD quickly tries to get a coarse idea of the full free energy surface (FES), and then slowly refines its details.
47 : It has a fast initial exploration phase, and then becomes extremely conservative and does not significantly change the shape of the deposited bias any more, reaching a regime of quasi-static bias.
48 : For this reason, it is possible to use standard umbrella sampling reweighting (see \ref REWEIGHT_BIAS) to analyse the trajectory.
49 : At <a href="https://github.com/invemichele/opes/tree/master/postprocessing">this link</a> you can find some python scripts that work in a similar way to \ref sum_hills, but the preferred way to obtain a FES with OPES is via reweighting (see \ref opes-metad).
50 : The estimated \f$c(t)\f$ is printed for reference only, since it should converge to a fixed value as the bias converges.
51 : This \f$c(t)\f$ should NOT be used for reweighting.
52 : Similarly, the \f$Z_n\f$ factor is printed only for reference, and it should converge when no new region of the CV-space is explored.
53 :
54 : Notice that \ref OPES_METAD is more sensitive to degenerate CVs than \ref METAD.
55 : If the employed CVs map different metastable basins onto the same CV-space region, then \ref OPES_METAD will remain stuck rather than completely reshaping the bias.
56 : This can be useful to diagnose problems with your collective variable.
57 : If it is not possible to improve the set of CVs and remove this degeneracy, then you might instead consider to use \ref OPES_METAD_EXPLORE or \ref METAD.
58 : In this way you will be able to obtain an estimate of the FES, but be aware that you most likely will not reach convergence and thus this estimate will be subjected to systematic errors (see e.g. Fig.3 in \cite Pietrucci2017review).
59 : On the contrary, if your CVs are not degenerate but only suboptimal, you should converge faster by using \ref OPES_METAD instead of \ref METAD \cite Invernizzi2020rethinking.
60 :
61 : The parameter BARRIER should be set to be at least equal to the highest free energy barrier you wish to overcome.
62 : If it is much lower than that, you will not cross the barrier, if it is much higher, convergence might take a little longer.
63 : If the system has a basin that is clearly more stable than the others, it is better to start the simulation from there.
64 :
65 : By default, the kernels SIGMA is adaptive, estimated from the fluctuations over ADAPTIVE_SIGMA_STRIDE simulation steps (similar to \ref METAD ADAPTIVE=DIFF, but contrary to that, no artifacts are introduced and the bias will converge to the correct one).
66 : However, notice that depending on the system this might not be the optimal choice for SIGMA.
67 :
68 : You can target a uniform flat distribution by explicitly setting BIASFACTOR=inf.
69 : However, this should be useful only in very specific cases.
70 :
71 : It is possible to take into account also of other bias potentials besides the one of \ref OPES_METAD during the internal reweighting for \f$P(\mathbf{s})\f$ estimation.
72 : To do so, one has to add those biases with the EXTRA_BIAS keyword, as in the example below.
73 : This allows one to define a custom target distribution by adding another bias potential equal to the desired target free energy and setting BIASFACTOR=inf (see example below).
74 : Another possible usage of EXTRA_BIAS is to make sure that \ref OPES_METAD does not push against another fixed bias added to restrain the CVs range (e.g. \ref UPPER_WALLS).
75 :
76 : Through the EXCLUDED_REGION keywork, it is possible to specify a region of CV space where no kernels will be deposited.
77 : This can be useful for example for making sure the bias does not modify the transition region, thus allowing for rate calculation.
78 : See below for an example of how to use this keyword.
79 :
80 : Restart can be done from a KERNELS file, but it might be not perfect (due to limited precision when printing kernels to file, or if adaptive SIGMA is used).
81 : For an exact restart you must use STATE_RFILE to read a checkpoint with all the needed info.
82 : To save such checkpoints, define a STATE_WFILE and choose how often to print them with STATE_WSTRIDE.
83 : By default this file is overwritten, but you can instead append to it using the flag STORE_STATES.
84 :
85 : Multiple walkers are supported only with MPI communication, via the keyword WALKERS_MPI.
86 :
87 : \par Examples
88 :
89 : Several examples can be found on the <a href="https://www.plumed-nest.org/browse.html">PLUMED-NEST website</a>, by searching for the OPES keyword.
90 : The \ref opes-metad can also be useful to get started with the method.
91 :
92 : The following is a minimal working example:
93 :
94 : \plumedfile
95 : cv: DISTANCE ATOMS=1,2
96 : opes: OPES_METAD ARG=cv PACE=200 BARRIER=40
97 : PRINT STRIDE=200 FILE=COLVAR ARG=*
98 : \endplumedfile
99 :
100 : Another more articulated one:
101 :
102 : \plumedfile
103 : phi: TORSION ATOMS=5,7,9,15
104 : psi: TORSION ATOMS=7,9,15,17
105 : opes: OPES_METAD ...
106 : FILE=Kernels.data
107 : TEMP=300
108 : ARG=phi,psi
109 : PACE=500
110 : BARRIER=50
111 : SIGMA=0.15,0.15
112 : SIGMA_MIN=0.01,0.01
113 : STATE_RFILE=Restart.data
114 : STATE_WFILE=State.data
115 : STATE_WSTRIDE=500*100
116 : STORE_STATES
117 : WALKERS_MPI
118 : NLIST
119 : ...
120 : PRINT FMT=%g STRIDE=500 FILE=Colvar.data ARG=phi,psi,opes.*
121 : \endplumedfile
122 :
123 : Next is an example of how to define a custom target distribution different from the well-tempered one.
124 : Here we chose to focus more on the transition state, that is around \f$\phi=0\f$.
125 : Our target distribution is a Gaussian centered there, thus the target free energy we want to sample is a parabola, \f$F^{\text{tg}}(\mathbf{s})=-\frac{1}{\beta} \log [p^{\text{tg}}(\mathbf{s})]\f$.
126 :
127 : \plumedfile
128 : phi: TORSION ATOMS=5,7,9,15
129 : FtgValue: CUSTOM ARG=phi PERIODIC=NO FUNC=(x/0.4)^2
130 : Ftg: BIASVALUE ARG=FtgValue
131 : opes: OPES_METAD ...
132 : ARG=phi
133 : PACE=500
134 : BARRIER=50
135 : SIGMA=0.2
136 : BIASFACTOR=inf
137 : EXTRA_BIAS=Ftg.bias
138 : ...
139 : PRINT FMT=%g STRIDE=500 FILE=COLVAR ARG=phi,Ftg.bias,opes.bias
140 : \endplumedfile
141 :
142 : Notice that in order to reweight for the unbiased \f$P(\mathbf{s})\f$ during postprocessing, the total bias `Ftg.bias+opes.bias` must be used.
143 :
144 : Finally, an example of how to use the EXCLUDED_REGION keyword.
145 : It expects a characteristic function that is different from zero in the region to be excluded.
146 : You can use \ref CUSTOM and a combination of the step function to define it.
147 : With the following input no kernel is deposited in the transition state region of alanine dipeptide, defined by the interval \f$\phi \in [-0.6, 0.7]\f$:
148 :
149 : \plumedfile
150 : phi: TORSION ATOMS=5,7,9,15
151 : psi: TORSION ATOMS=7,9,15,17
152 : xx: CUSTOM PERIODIC=NO ARG=phi FUNC=step(x+0.6)-step(x-0.7)
153 : opes: OPES_METAD ...
154 : ARG=phi,psi
155 : PACE=500
156 : BARRIER=30
157 : EXCLUDED_REGION=xx
158 : NLIST
159 : ...
160 : PRINT FMT=%g STRIDE=500 FILE=COLVAR ARG=phi,psi,xx,opes.*
161 : \endplumedfile
162 :
163 : */
164 : //+ENDPLUMEDOC
165 :
166 : template <class mode>
167 : class OPESmetad : public bias::Bias {
168 :
169 : private:
170 : bool isFirstStep_;
171 : unsigned NumOMP_;
172 : unsigned NumParallel_;
173 : unsigned rank_;
174 : unsigned NumWalkers_;
175 : unsigned walker_rank_;
176 : unsigned long long counter_;
177 : std::size_t ncv_;
178 :
179 : double kbt_;
180 : double biasfactor_;
181 : double bias_prefactor_;
182 : unsigned stride_;
183 : std::vector<double> sigma0_;
184 : std::vector<double> sigma_min_;
185 : unsigned adaptive_sigma_stride_;
186 : unsigned long long adaptive_counter_;
187 : std::vector<double> av_cv_;
188 : std::vector<double> av_M2_;
189 : bool fixed_sigma_;
190 : bool adaptive_sigma_;
191 : double epsilon_;
192 : double sum_weights_;
193 : double sum_weights2_;
194 :
195 : bool no_Zed_;
196 : double Zed_;
197 : double KDEnorm_;
198 :
199 : double threshold2_;
200 : bool recursive_merge_;
201 : //kernels are truncated diagonal Gaussians
202 1050 : struct kernel
203 : {
204 : double height;
205 : std::vector<double> center;
206 : std::vector<double> sigma;
207 1050 : kernel(double h, const std::vector<double>& c,const std::vector<double>& s):
208 1050 : height(h),center(c),sigma(s) {}
209 : };
210 : double cutoff2_;
211 : double val_at_cutoff_;
212 : void mergeKernels(kernel&,const kernel&); //merge the second one into the first one
213 : double evaluateKernel(const kernel&,const std::vector<double>&) const;
214 : double evaluateKernel(const kernel&,const std::vector<double>&,std::vector<double>&,std::vector<double>&);
215 : std::vector<kernel> kernels_; //all compressed kernels
216 : OFile kernelsOfile_;
217 : //neighbour list stuff
218 : bool nlist_;
219 : double nlist_param_[2];
220 : std::vector<unsigned> nlist_index_;
221 : std::vector<double> nlist_center_;
222 : std::vector<double> nlist_dev2_;
223 : unsigned nlist_steps_;
224 : bool nlist_update_;
225 : bool nlist_pace_reset_;
226 :
227 : bool calc_work_;
228 : double work_;
229 : double old_KDEnorm_;
230 : std::vector<kernel> delta_kernels_;
231 :
232 : Value* excluded_region_;
233 : std::vector<Value*> extra_biases_;
234 :
235 : OFile stateOfile_;
236 : int wStateStride_;
237 : bool storeOldStates_;
238 :
239 : double getProbAndDerivatives(const std::vector<double>&,std::vector<double>&);
240 : void addKernel(const double,const std::vector<double>&,const std::vector<double>&);
241 : void addKernel(const double,const std::vector<double>&,const std::vector<double>&,const double); //also print to file
242 : unsigned getMergeableKernel(const std::vector<double>&,const unsigned);
243 : void updateNlist(const std::vector<double>&);
244 : void dumpStateToFile();
245 :
246 : public:
247 : explicit OPESmetad(const ActionOptions&);
248 : void calculate() override;
249 : void update() override;
250 : static void registerKeywords(Keywords& keys);
251 : };
252 :
253 : struct convergence { static const bool explore=false; };
254 : typedef OPESmetad<convergence> OPESmetad_c;
255 : PLUMED_REGISTER_ACTION(OPESmetad_c,"OPES_METAD")
256 :
257 : //OPES_METAD_EXPLORE is very similar from the point of view of the code,
258 : //but conceptually it is better to make it a separate BIAS action
259 :
260 : //+PLUMEDOC OPES_BIAS OPES_METAD_EXPLORE
261 : /*
262 : On-the-fly probability enhanced sampling with well-tempered target distribution in exploreation mode.
263 :
264 : On-the-fly probability enhanced sampling with well-tempered target distribution (\ref OPES "OPES") with well-tempered target distribution, exploration mode \cite Invernizzi2022explore .
265 :
266 : This \ref OPES_METAD_EXPLORE action samples the well-tempered target distribution, that is defined via its marginal \f$p^{\text{WT}}(\mathbf{s})\propto [P(\mathbf{s})]^{1/\gamma}\f$ over some collective variables (CVs), \f$\mathbf{s}=\mathbf{s}(\mathbf{x})\f$.
267 : While \ref OPES_METAD does so by estimating the unbiased distribution \f$P(\mathbf{s})\f$, \ref OPES_METAD_EXPLORE instead estimates on-the-fly the target \f$p^{\text{WT}}(\mathbf{s})\f$ and uses it to define the bias.
268 : The bias at step \f$n\f$ is
269 : \f[
270 : V_n(\mathbf{s}) = (\gamma-1)\frac{1}{\beta}\log\left(\frac{p^{\text{WT}}_n(\mathbf{s})}{Z_n}+\epsilon\right)\, .
271 : \f]
272 : See Ref.\cite Invernizzi2022explore for a complete description of the method.
273 :
274 : Intuitively, while \ref OPES_METAD aims at quickly converging the reweighted free energy, \ref OPES_METAD_EXPLORE aims at quickly sampling the target well-tempered distribution.
275 : Given enough simulation time, both will converge to the same bias potential but they do so in a qualitatively different way.
276 : Compared to \ref OPES_METAD, \ref OPES_METAD_EXPLORE is more similar to \ref METAD, because it allows the bias to vary significantly, thus enhancing exploration.
277 : This goes at the expenses of a typically slower convergence of the reweight estimate.
278 : \ref OPES_METAD_EXPLORE can be useful e.g.~for simulating a new system with an unknown BARRIER, or for quickly testing the effectiveness of a new CV that might be degenerate.
279 :
280 : Similarly to \ref OPES_METAD, also \ref OPES_METAD_EXPLORE uses a kernel density estimation with the same on-the-fly compression algorithm.
281 : The only difference is that the kernels are not weighted, since it estimates the sampled distribution and not the reweighted unbiased one.
282 :
283 : All the options of \ref OPES_METAD are also available in \ref OPES_METAD_EXPLORE, except for those that modify the target distribution, since only a well-tempered target is allowed in this case.
284 :
285 : \par Examples
286 :
287 : The following is a minimal working example:
288 :
289 : \plumedfile
290 : cv: DISTANCE ATOMS=1,2
291 : opes: OPES_METAD_EXPLORE ARG=cv PACE=500 BARRIER=40
292 : PRINT STRIDE=100 FILE=COLVAR ARG=cv,opes.*
293 : \endplumedfile
294 : */
295 : //+ENDPLUMEDOC
296 :
297 : struct exploration { static const bool explore=true; };
298 : typedef OPESmetad<exploration> OPESmetad_e;
299 : PLUMED_REGISTER_ACTION(OPESmetad_e,"OPES_METAD_EXPLORE")
300 :
301 : template <class mode>
302 25 : void OPESmetad<mode>::registerKeywords(Keywords& keys)
303 : {
304 25 : Bias::registerKeywords(keys);
305 50 : keys.add("compulsory","TEMP","-1","temperature. If not set, it is taken from MD engine, but not all MD codes provide it");
306 50 : keys.add("compulsory","PACE","the frequency for kernel deposition");
307 25 : std::string info_sigma("the initial widths of the kernels");
308 : if(mode::explore)
309 : info_sigma+=", divided by the square root of gamma";
310 : info_sigma+=". If not set, an adaptive sigma will be used with the given ADAPTIVE_SIGMA_STRIDE";
311 50 : keys.add("compulsory","SIGMA","ADAPTIVE",info_sigma);
312 50 : keys.add("compulsory","BARRIER","the free energy barrier to be overcome. It is used to set BIASFACTOR, EPSILON, and KERNEL_CUTOFF to reasonable values");
313 50 : keys.add("compulsory","COMPRESSION_THRESHOLD","1","merge kernels if closer than this threshold, in units of sigma");
314 : //extra options
315 50 : keys.add("optional","ADAPTIVE_SIGMA_STRIDE","number of steps for measuring adaptive sigma. Default is 10xPACE");
316 50 : keys.add("optional","SIGMA_MIN","never reduce SIGMA below this value");
317 25 : std::string info_biasfactor("the gamma bias factor used for the well-tempered target distribution. ");
318 : if(mode::explore)
319 : info_biasfactor+="Cannot be 'inf'";
320 : else
321 : info_biasfactor+="Set to 'inf' for uniform flat target";
322 50 : keys.add("optional","BIASFACTOR",info_biasfactor);
323 50 : keys.add("optional","EPSILON","the value of the regularization constant for the probability");
324 50 : keys.add("optional","KERNEL_CUTOFF","truncate kernels at this distance, in units of sigma");
325 50 : keys.add("optional","NLIST_PARAMETERS","( default=3.0,0.5 ) the two cutoff parameters for the kernels neighbor list");
326 50 : keys.addFlag("NLIST",false,"use neighbor list for kernels summation, faster but experimental");
327 50 : keys.addFlag("NLIST_PACE_RESET",false,"force the reset of the neighbor list at each PACE. Can be useful with WALKERS_MPI");
328 50 : keys.addFlag("FIXED_SIGMA",false,"do not decrease sigma as the simulation proceeds. Can be added in a RESTART, to keep in check the number of compressed kernels");
329 50 : keys.addFlag("RECURSIVE_MERGE_OFF",false,"do not recursively attempt kernel merging when a new one is added");
330 50 : keys.addFlag("NO_ZED",false,"do not normalize over the explored CV space, Z_n=1");
331 : //kernels and state files
332 50 : keys.add("compulsory","FILE","KERNELS","a file in which the list of all deposited kernels is stored");
333 50 : keys.add("optional","FMT","specify format for KERNELS file");
334 50 : keys.add("optional","STATE_RFILE","read from this file the compressed kernels and all the info needed to RESTART the simulation");
335 50 : keys.add("optional","STATE_WFILE","write to this file the compressed kernels and all the info needed to RESTART the simulation");
336 50 : keys.add("optional","STATE_WSTRIDE","number of MD steps between writing the STATE_WFILE. Default is only on CPT events (but not all MD codes set them)");
337 50 : keys.addFlag("STORE_STATES",false,"append to STATE_WFILE instead of ovewriting it each time");
338 : //miscellaneous
339 50 : keys.addInputKeyword("optional","EXCLUDED_REGION","scalar","kernels are not deposited when the action provided here has a nonzero value, see example above");
340 : if(!mode::explore)
341 32 : keys.addInputKeyword("optional","EXTRA_BIAS","scalar","consider also these other bias potentials for the internal reweighting. This can be used e.g. for sampling a custom target distribution (see example above)");
342 50 : keys.addFlag("CALC_WORK",false,"calculate the total accumulated work done by the bias since last restart");
343 50 : keys.addFlag("WALKERS_MPI",false,"switch on MPI version of multiple walkers");
344 50 : keys.addFlag("SERIAL",false,"perform calculations in serial");
345 25 : keys.use("RESTART");
346 25 : keys.use("UPDATE_FROM");
347 25 : keys.use("UPDATE_UNTIL");
348 :
349 : //output components
350 50 : keys.addOutputComponent("rct","default","scalar","estimate of c(t). log(exp(beta V)/beta, should become flat as the simulation converges. Do NOT use for reweighting");
351 50 : keys.addOutputComponent("zed","default","scalar","estimate of Z_n. should become flat once no new CV-space region is explored");
352 50 : keys.addOutputComponent("neff","default","scalar","effective sample size");
353 50 : keys.addOutputComponent("nker","default","scalar","total number of compressed kernels used to represent the bias");
354 50 : keys.addOutputComponent("work","CALC_WORK","scalar","total accumulated work done by the bias");
355 50 : keys.addOutputComponent("nlker","NLIST","scalar","number of kernels in the neighbor list");
356 50 : keys.addOutputComponent("nlsteps","NLIST","scalar","number of steps from last neighbor list update");
357 25 : }
358 :
359 : template <class mode>
360 21 : OPESmetad<mode>::OPESmetad(const ActionOptions& ao)
361 : : PLUMED_BIAS_INIT(ao)
362 21 : , isFirstStep_(true)
363 21 : , counter_(1)
364 21 : , ncv_(getNumberOfArguments())
365 21 : , Zed_(1)
366 21 : , work_(0)
367 21 : , excluded_region_(NULL)
368 : {
369 42 : std::string error_in_input1("Error in input in action "+getName()+" with label "+getLabel()+": the keyword ");
370 21 : std::string error_in_input2(" could not be read correctly");
371 :
372 : //set kbt_
373 21 : const double kB=getKBoltzmann();
374 21 : kbt_=getkBT();
375 :
376 : //other compulsory input
377 21 : parse("PACE",stride_);
378 :
379 21 : double barrier=0;
380 21 : parse("BARRIER",barrier);
381 21 : plumed_massert(barrier>=0,"the BARRIER should be greater than zero");
382 :
383 21 : biasfactor_=barrier/kbt_;
384 : std::string biasfactor_str;
385 42 : parse("BIASFACTOR",biasfactor_str);
386 38 : if(biasfactor_str=="inf" || biasfactor_str=="INF")
387 : {
388 4 : biasfactor_=std::numeric_limits<double>::infinity();
389 4 : bias_prefactor_=1;
390 : }
391 : else
392 : {
393 17 : if(biasfactor_str.length()>0)
394 3 : plumed_massert(Tools::convertNoexcept(biasfactor_str,biasfactor_),error_in_input1+"BIASFACTOR"+error_in_input2);
395 17 : plumed_massert(biasfactor_>1,"BIASFACTOR must be greater than one (use 'inf' for uniform target)");
396 17 : bias_prefactor_=1-1./biasfactor_;
397 : }
398 : if(mode::explore)
399 : {
400 7 : plumed_massert(!std::isinf(biasfactor_),"BIASFACTOR=inf is not compatible with EXPLORE mode");
401 7 : bias_prefactor_=biasfactor_-1;
402 : }
403 :
404 21 : adaptive_sigma_=false;
405 21 : adaptive_sigma_stride_=0;
406 42 : parse("ADAPTIVE_SIGMA_STRIDE",adaptive_sigma_stride_);
407 : std::vector<std::string> sigma_str;
408 21 : parseVector("SIGMA",sigma_str);
409 21 : sigma0_.resize(ncv_);
410 : double dummy;
411 21 : if(sigma_str.size()==1 && !Tools::convertNoexcept(sigma_str[0],dummy))
412 : {
413 11 : plumed_massert(sigma_str[0]=="ADAPTIVE" || sigma_str[0]=="adaptive",error_in_input1+"SIGMA"+error_in_input2);
414 11 : plumed_massert(!std::isinf(biasfactor_),"cannot use BIASFACTOR=inf with adaptive SIGMA");
415 11 : adaptive_counter_=0;
416 11 : if(adaptive_sigma_stride_==0)
417 2 : adaptive_sigma_stride_=10*stride_; //NB: this is arbitrary, chosen from few tests
418 11 : av_cv_.resize(ncv_,0);
419 11 : av_M2_.resize(ncv_,0);
420 11 : plumed_massert(adaptive_sigma_stride_>=stride_,"better to chose ADAPTIVE_SIGMA_STRIDE > PACE");
421 11 : adaptive_sigma_=true;
422 : }
423 : else
424 : {
425 10 : plumed_massert(sigma_str.size()==ncv_,"number of SIGMA parameters does not match number of arguments");
426 10 : plumed_massert(adaptive_sigma_stride_==0,"if SIGMA is not ADAPTIVE you cannot set an ADAPTIVE_SIGMA_STRIDE");
427 29 : for(unsigned i=0; i<ncv_; i++)
428 : {
429 19 : plumed_massert(Tools::convertNoexcept(sigma_str[i],sigma0_[i]),error_in_input1+"SIGMA"+error_in_input2);
430 : if(mode::explore)
431 6 : sigma0_[i]*=std::sqrt(biasfactor_); //the sigma of the target is broader Ftg(s)=1/gamma*F(s)
432 : }
433 : }
434 42 : parseVector("SIGMA_MIN",sigma_min_);
435 21 : plumed_massert(sigma_min_.size()==0 || sigma_min_.size()==ncv_,"number of SIGMA_MIN does not match number of arguments");
436 21 : if(sigma_min_.size()>0 && !adaptive_sigma_)
437 : {
438 3 : for(unsigned i=0; i<ncv_; i++)
439 2 : plumed_massert(sigma_min_[i]<=sigma0_[i],"SIGMA_MIN should be smaller than SIGMA");
440 : }
441 :
442 21 : epsilon_=std::exp(-barrier/bias_prefactor_/kbt_);
443 21 : parse("EPSILON",epsilon_);
444 21 : plumed_massert(epsilon_>0,"you must choose a value for EPSILON greater than zero. Is your BARRIER too high?");
445 21 : sum_weights_=std::pow(epsilon_,bias_prefactor_); //to avoid NANs we start with counter_=1 and w0=exp(beta*V0)
446 21 : sum_weights2_=sum_weights_*sum_weights_;
447 :
448 21 : double cutoff=sqrt(2.*barrier/bias_prefactor_/kbt_);
449 : if(mode::explore)
450 7 : cutoff=sqrt(2.*barrier/kbt_); //otherwise it is too small
451 21 : parse("KERNEL_CUTOFF",cutoff);
452 21 : plumed_massert(cutoff>0,"you must choose a value for KERNEL_CUTOFF greater than zero");
453 21 : cutoff2_=cutoff*cutoff;
454 21 : val_at_cutoff_=std::exp(-0.5*cutoff2_);
455 :
456 21 : threshold2_=1;
457 21 : parse("COMPRESSION_THRESHOLD",threshold2_);
458 21 : threshold2_*=threshold2_;
459 21 : if(threshold2_!=0)
460 21 : plumed_massert(threshold2_>0 && threshold2_<cutoff2_,"COMPRESSION_THRESHOLD cannot be bigger than the KERNEL_CUTOFF");
461 :
462 : //setup neighbor list
463 21 : nlist_=false;
464 21 : parseFlag("NLIST",nlist_);
465 21 : nlist_pace_reset_=false;
466 21 : parseFlag("NLIST_PACE_RESET",nlist_pace_reset_);
467 21 : if(nlist_pace_reset_)
468 2 : nlist_=true;
469 : std::vector<double> nlist_param;
470 42 : parseVector("NLIST_PARAMETERS",nlist_param);
471 21 : if(nlist_param.size()==0)
472 : {
473 17 : nlist_param_[0]=3.0;//*cutoff2_ -> max distance of neighbors
474 17 : nlist_param_[1]=0.5;//*nlist_dev2_[i] -> condition for rebuilding
475 : }
476 : else
477 : {
478 4 : nlist_=true;
479 4 : plumed_massert(nlist_param.size()==2,"two cutoff parameters are needed for the neighbor list");
480 4 : plumed_massert(nlist_param[0]>1.0,"the first of NLIST_PARAMETERS must be greater than 1. The smaller the first, the smaller should be the second as well");
481 4 : const double min_PARAM_1=(1.-1./std::sqrt(nlist_param[0]))+0.16;
482 4 : plumed_massert(nlist_param[1]>0,"the second of NLIST_PARAMETERS must be greater than 0");
483 4 : plumed_massert(nlist_param[1]<=min_PARAM_1,"the second of NLIST_PARAMETERS must be smaller to avoid systematic errors. Largest suggested value is: 1.16-1/sqrt(PARAM_0) = "+std::to_string(min_PARAM_1));
484 4 : nlist_param_[0]=nlist_param[0];
485 4 : nlist_param_[1]=nlist_param[1];
486 : }
487 21 : nlist_center_.resize(ncv_);
488 21 : nlist_dev2_.resize(ncv_,0.);
489 21 : nlist_steps_=0;
490 21 : nlist_update_=true;
491 :
492 : //optional stuff
493 21 : no_Zed_=false;
494 21 : parseFlag("NO_ZED",no_Zed_);
495 21 : if(no_Zed_)
496 : { //this makes it more gentle in the initial phase
497 6 : sum_weights_=1;
498 6 : sum_weights2_=1;
499 : }
500 21 : fixed_sigma_=false;
501 21 : parseFlag("FIXED_SIGMA",fixed_sigma_);
502 21 : bool recursive_merge_off=false;
503 21 : parseFlag("RECURSIVE_MERGE_OFF",recursive_merge_off);
504 21 : recursive_merge_=!recursive_merge_off;
505 42 : parseFlag("CALC_WORK",calc_work_);
506 :
507 : //options involving extra arguments
508 : std::vector<Value*> args;
509 42 : parseArgumentList("EXCLUDED_REGION",args);
510 21 : if(args.size()>0)
511 : {
512 2 : plumed_massert(args.size()==1,"only one characteristic function should define the region to be excluded");
513 2 : requestExtraDependencies(args);
514 2 : excluded_region_=args[0];
515 : }
516 : if(!mode::explore)
517 : {
518 28 : parseArgumentList("EXTRA_BIAS",extra_biases_);
519 14 : if(extra_biases_.size()>0)
520 2 : requestExtraDependencies(extra_biases_);
521 : }
522 :
523 : //kernels file
524 : std::string kernelsFileName;
525 42 : parse("FILE",kernelsFileName);
526 : std::string fmt;
527 42 : parse("FMT",fmt);
528 :
529 : //output checkpoint of current state
530 : std::string restartFileName;
531 42 : parse("STATE_RFILE",restartFileName);
532 : std::string stateFileName;
533 21 : parse("STATE_WFILE",stateFileName);
534 21 : wStateStride_=0;
535 21 : parse("STATE_WSTRIDE",wStateStride_);
536 21 : storeOldStates_=false;
537 21 : parseFlag("STORE_STATES",storeOldStates_);
538 21 : if(wStateStride_!=0 || storeOldStates_)
539 10 : plumed_massert(stateFileName.length()>0,"filename for storing simulation status not specified, use STATE_WFILE");
540 21 : if(wStateStride_>0)
541 10 : plumed_massert(wStateStride_>=(int)stride_,"STATE_WSTRIDE is in units of MD steps, thus it is suggested to use a multiple of PACE");
542 21 : if(stateFileName.length()>0 && wStateStride_==0)
543 1 : wStateStride_=-1;//will print only on CPT events (checkpoints set by some MD engines, like gromacs)
544 :
545 : //multiple walkers //TODO implement also external mw for cp2k
546 21 : bool walkers_mpi=false;
547 21 : parseFlag("WALKERS_MPI",walkers_mpi);
548 21 : if(walkers_mpi)
549 : {
550 : //If this Action is not compiled with MPI the user is informed and we exit gracefully
551 10 : plumed_massert(Communicator::plumedHasMPI(),"Invalid walkers configuration: WALKERS_MPI flag requires MPI compilation");
552 10 : plumed_massert(Communicator::initialized(),"Invalid walkers configuration: WALKERS_MPI needs the communicator correctly initialized.");
553 :
554 10 : if(comm.Get_rank()==0)//multi_sim_comm works on first rank only
555 : {
556 10 : NumWalkers_=multi_sim_comm.Get_size();
557 10 : walker_rank_=multi_sim_comm.Get_rank();
558 : }
559 10 : comm.Bcast(NumWalkers_,0); //if each walker has more than one processor update them all
560 10 : comm.Bcast(walker_rank_,0);
561 : }
562 : else
563 : {
564 11 : NumWalkers_=1;
565 11 : walker_rank_=0;
566 : }
567 :
568 : //parallelization stuff
569 21 : NumOMP_=OpenMP::getNumThreads();
570 21 : NumParallel_=comm.Get_size();
571 21 : rank_=comm.Get_rank();
572 21 : bool serial=false;
573 21 : parseFlag("SERIAL",serial);
574 21 : if(serial)
575 : {
576 4 : NumOMP_=1;
577 4 : NumParallel_=1;
578 4 : rank_=0;
579 : }
580 :
581 21 : checkRead();
582 :
583 : //restart if needed
584 : bool convertKernelsToState=false;
585 21 : if(getRestart())
586 : {
587 : bool stateRestart=true;
588 11 : if(restartFileName.length()==0)
589 : {
590 : stateRestart=false;
591 : restartFileName=kernelsFileName;
592 : }
593 11 : IFile ifile;
594 11 : ifile.link(*this);
595 11 : if(ifile.FileExist(restartFileName))
596 : {
597 11 : bool tmp_nlist=nlist_;
598 11 : nlist_=false; // NLIST is not needed while restarting
599 11 : ifile.open(restartFileName);
600 11 : log.printf(" RESTART - make sure all used options are compatible\n");
601 11 : log.printf(" restarting from: %s\n",restartFileName.c_str());
602 11 : std::string action_name=getName();
603 11 : if(stateRestart)
604 : {
605 6 : log.printf(" it should be a STATE file (not a KERNELS file)\n");
606 : action_name+="_state";
607 : }
608 : else
609 : {
610 5 : log.printf(" +++ WARNING +++ RESTART from KERNELS might be approximate, use STATE_WFILE and STATE_RFILE to restart from the exact state\n");
611 : action_name+="_kernels";
612 : }
613 : std::string old_action_name;
614 11 : ifile.scanField("action",old_action_name);
615 11 : plumed_massert(action_name==old_action_name,"RESTART - mismatch between old and new action name. Expected '"+action_name+"', but found '"+old_action_name+"'");
616 : std::string old_biasfactor_str;
617 22 : ifile.scanField("biasfactor",old_biasfactor_str);
618 21 : if(old_biasfactor_str=="inf" || old_biasfactor_str=="INF")
619 : {
620 1 : if(!std::isinf(biasfactor_))
621 0 : log.printf(" +++ WARNING +++ previous bias factor was inf while now it is %g\n",biasfactor_);
622 : }
623 : else
624 : {
625 : double old_biasfactor;
626 10 : ifile.scanField("biasfactor",old_biasfactor);
627 10 : if(std::abs(biasfactor_-old_biasfactor)>1e-6*biasfactor_)
628 0 : log.printf(" +++ WARNING +++ previous bias factor was %g while now it is %g. diff = %g\n",old_biasfactor,biasfactor_,biasfactor_-old_biasfactor);
629 : }
630 : double old_epsilon;
631 11 : ifile.scanField("epsilon",old_epsilon);
632 11 : if(std::abs(epsilon_-old_epsilon)>1e-6*epsilon_)
633 8 : log.printf(" +++ WARNING +++ previous epsilon was %g while now it is %g. diff = %g\n",old_epsilon,epsilon_,epsilon_-old_epsilon);
634 : double old_cutoff;
635 11 : ifile.scanField("kernel_cutoff",old_cutoff);
636 11 : if(std::abs(cutoff-old_cutoff)>1e-6*cutoff)
637 0 : log.printf(" +++ WARNING +++ previous kernel_cutoff was %g while now it is %g. diff = %g\n",old_cutoff,cutoff,cutoff-old_cutoff);
638 : double old_threshold;
639 11 : const double threshold=sqrt(threshold2_);
640 11 : ifile.scanField("compression_threshold",old_threshold);
641 11 : if(std::abs(threshold-old_threshold)>1e-6*threshold)
642 0 : log.printf(" +++ WARNING +++ previous compression_threshold was %g while now it is %g. diff = %g\n",old_threshold,threshold,threshold-old_threshold);
643 11 : if(stateRestart)
644 : {
645 6 : ifile.scanField("zed",Zed_);
646 6 : ifile.scanField("sum_weights",sum_weights_);
647 6 : ifile.scanField("sum_weights2",sum_weights2_);
648 6 : ifile.scanField("counter",counter_);
649 6 : if(adaptive_sigma_)
650 : {
651 6 : ifile.scanField("adaptive_counter",adaptive_counter_);
652 6 : if(NumWalkers_==1)
653 : {
654 6 : for(unsigned i=0; i<ncv_; i++)
655 : {
656 8 : ifile.scanField("sigma0_"+getPntrToArgument(i)->getName(),sigma0_[i]);
657 8 : ifile.scanField("av_cv_"+getPntrToArgument(i)->getName(),av_cv_[i]);
658 8 : ifile.scanField("av_M2_"+getPntrToArgument(i)->getName(),av_M2_[i]);
659 : }
660 : }
661 : else
662 : {
663 12 : for(unsigned w=0; w<NumWalkers_; w++)
664 24 : for(unsigned i=0; i<ncv_; i++)
665 : {
666 : double tmp0,tmp1,tmp2;
667 32 : const std::string arg_iw=getPntrToArgument(i)->getName()+"_"+std::to_string(w);
668 16 : ifile.scanField("sigma0_"+arg_iw,tmp0);
669 16 : ifile.scanField("av_cv_"+arg_iw,tmp1);
670 16 : ifile.scanField("av_M2_"+arg_iw,tmp2);
671 16 : if(w==walker_rank_)
672 : {
673 8 : sigma0_[i]=tmp0;
674 8 : av_cv_[i]=tmp1;
675 8 : av_M2_[i]=tmp2;
676 : }
677 : }
678 : }
679 : }
680 : }
681 33 : for(unsigned i=0; i<ncv_; i++)
682 : {
683 22 : if(getPntrToArgument(i)->isPeriodic())
684 : {
685 : std::string arg_min,arg_max;
686 22 : getPntrToArgument(i)->getDomain(arg_min,arg_max);
687 : std::string file_min,file_max;
688 44 : ifile.scanField("min_"+getPntrToArgument(i)->getName(),file_min);
689 22 : ifile.scanField("max_"+getPntrToArgument(i)->getName(),file_max);
690 22 : plumed_massert(file_min==arg_min,"RESTART - mismatch between old and new ARG periodicity");
691 22 : plumed_massert(file_max==arg_max,"RESTART - mismatch between old and new ARG periodicity");
692 : }
693 : }
694 11 : if(stateRestart)
695 : {
696 : double time;
697 60 : while(ifile.scanField("time",time))
698 : {
699 24 : std::vector<double> center(ncv_);
700 24 : std::vector<double> sigma(ncv_);
701 : double height;
702 72 : for(unsigned i=0; i<ncv_; i++)
703 48 : ifile.scanField(getPntrToArgument(i)->getName(),center[i]);
704 72 : for(unsigned i=0; i<ncv_; i++)
705 96 : ifile.scanField("sigma_"+getPntrToArgument(i)->getName(),sigma[i]);
706 24 : ifile.scanField("height",height);
707 24 : ifile.scanField();
708 24 : kernels_.emplace_back(height,center,sigma);
709 : }
710 6 : log.printf(" a total of %lu kernels where read\n",kernels_.size());
711 : }
712 : else
713 : {
714 5 : ifile.allowIgnoredFields(); //this allows for multiple restart, but without checking for consistency between them!
715 : double time;
716 270 : while(ifile.scanField("time",time))
717 : {
718 130 : std::vector<double> center(ncv_);
719 130 : std::vector<double> sigma(ncv_);
720 : double height;
721 : double logweight;
722 390 : for(unsigned i=0; i<ncv_; i++)
723 260 : ifile.scanField(getPntrToArgument(i)->getName(),center[i]);
724 390 : for(unsigned i=0; i<ncv_; i++)
725 520 : ifile.scanField("sigma_"+getPntrToArgument(i)->getName(),sigma[i]);
726 130 : if(counter_==(1+walker_rank_) && adaptive_sigma_)
727 0 : sigma0_=sigma;
728 130 : ifile.scanField("height",height);
729 130 : ifile.scanField("logweight",logweight);
730 130 : ifile.scanField();
731 130 : addKernel(height,center,sigma);
732 130 : const double weight=std::exp(logweight);
733 130 : sum_weights_+=weight; //this sum is slightly inaccurate, because when printing some precision is lost
734 130 : sum_weights2_+=weight*weight;
735 130 : counter_++;
736 : }
737 5 : KDEnorm_=mode::explore?counter_:sum_weights_;
738 5 : if(!no_Zed_)
739 : {
740 2 : double sum_uprob=0;
741 48 : for(unsigned k=rank_; k<kernels_.size(); k+=NumParallel_)
742 1104 : for(unsigned kk=0; kk<kernels_.size(); kk++)
743 1058 : sum_uprob+=evaluateKernel(kernels_[kk],kernels_[k].center);
744 2 : if(NumParallel_>1)
745 0 : comm.Sum(sum_uprob);
746 2 : Zed_=sum_uprob/KDEnorm_/kernels_.size();
747 : }
748 5 : log.printf(" a total of %llu kernels where read, and compressed to %lu\n",counter_-1,kernels_.size());
749 : convertKernelsToState=true;
750 : }
751 11 : ifile.reset(false);
752 11 : ifile.close();
753 11 : nlist_=tmp_nlist;
754 : }
755 : else //same behaviour as METAD
756 0 : plumed_merror("RESTART requested, but file '"+restartFileName+"' was not found!\n Set RESTART=NO or provide a restart file");
757 11 : if(NumWalkers_>1) //make sure that all walkers are doing the same thing
758 : {
759 6 : const unsigned kernels_size=kernels_.size();
760 6 : std::vector<unsigned> all_kernels_size(NumWalkers_);
761 6 : if(comm.Get_rank()==0)
762 6 : multi_sim_comm.Allgather(kernels_size,all_kernels_size);
763 6 : comm.Bcast(all_kernels_size,0);
764 : bool same_number_of_kernels=true;
765 12 : for(unsigned w=1; w<NumWalkers_; w++)
766 6 : if(all_kernels_size[0]!=all_kernels_size[w])
767 : same_number_of_kernels=false;
768 6 : plumed_massert(same_number_of_kernels,"RESTART - not all walkers are reading the same file!");
769 : }
770 11 : }
771 10 : else if(restartFileName.length()>0)
772 4 : log.printf(" +++ WARNING +++ the provided STATE_RFILE will be ignored, since RESTART was not requested\n");
773 :
774 : //sync all walkers to avoid opening files before reading is over (see also METAD)
775 21 : comm.Barrier();
776 21 : if(comm.Get_rank()==0 && walkers_mpi)
777 10 : multi_sim_comm.Barrier();
778 :
779 : //setup output kernels file
780 21 : kernelsOfile_.link(*this);
781 21 : if(NumWalkers_>1)
782 : {
783 10 : if(walker_rank_>0)
784 : kernelsFileName="/dev/null"; //only first walker writes on file
785 20 : kernelsOfile_.enforceSuffix("");
786 : }
787 21 : kernelsOfile_.open(kernelsFileName);
788 21 : if(fmt.length()>0)
789 42 : kernelsOfile_.fmtField(" "+fmt);
790 : kernelsOfile_.setHeavyFlush(); //do I need it?
791 : //define and set const fields
792 21 : kernelsOfile_.addConstantField("action");
793 21 : kernelsOfile_.addConstantField("biasfactor");
794 21 : kernelsOfile_.addConstantField("epsilon");
795 21 : kernelsOfile_.addConstantField("kernel_cutoff");
796 21 : kernelsOfile_.addConstantField("compression_threshold");
797 61 : for(unsigned i=0; i<ncv_; i++)
798 40 : kernelsOfile_.setupPrintValue(getPntrToArgument(i));
799 42 : kernelsOfile_.printField("action",getName()+"_kernels");
800 21 : kernelsOfile_.printField("biasfactor",biasfactor_);
801 21 : kernelsOfile_.printField("epsilon",epsilon_);
802 21 : kernelsOfile_.printField("kernel_cutoff",sqrt(cutoff2_));
803 21 : kernelsOfile_.printField("compression_threshold",sqrt(threshold2_));
804 :
805 : //open file for storing state
806 21 : if(wStateStride_!=0)
807 : {
808 11 : stateOfile_.link(*this);
809 11 : if(NumWalkers_>1)
810 : {
811 8 : if(walker_rank_>0)
812 : stateFileName="/dev/null"; //only first walker writes on file
813 16 : stateOfile_.enforceSuffix("");
814 : }
815 11 : stateOfile_.open(stateFileName);
816 11 : if(fmt.length()>0)
817 22 : stateOfile_.fmtField(" "+fmt);
818 11 : if(convertKernelsToState)
819 0 : dumpStateToFile();
820 : }
821 :
822 : //set initial old values
823 21 : KDEnorm_=mode::explore?counter_:sum_weights_;
824 21 : old_KDEnorm_=KDEnorm_;
825 :
826 : //add and set output components
827 42 : addComponent("rct");
828 21 : componentIsNotPeriodic("rct");
829 42 : getPntrToComponent("rct")->set(kbt_*std::log(sum_weights_/counter_));
830 42 : addComponent("zed");
831 21 : componentIsNotPeriodic("zed");
832 42 : getPntrToComponent("zed")->set(Zed_);
833 42 : addComponent("neff");
834 21 : componentIsNotPeriodic("neff");
835 42 : getPntrToComponent("neff")->set(std::pow(1+sum_weights_,2)/(1+sum_weights2_));
836 42 : addComponent("nker");
837 21 : componentIsNotPeriodic("nker");
838 21 : getPntrToComponent("nker")->set(kernels_.size());
839 21 : if(calc_work_)
840 : {
841 14 : addComponent("work");
842 14 : componentIsNotPeriodic("work");
843 : }
844 21 : if(nlist_)
845 : {
846 10 : addComponent("nlker");
847 10 : componentIsNotPeriodic("nlker");
848 10 : addComponent("nlsteps");
849 10 : componentIsNotPeriodic("nlsteps");
850 : }
851 :
852 : //printing some info
853 21 : log.printf(" temperature = %g\n",kbt_/kB);
854 21 : log.printf(" beta = %g\n",1./kbt_);
855 21 : log.printf(" depositing new kernels with PACE = %u\n",stride_);
856 21 : log.printf(" expected BARRIER is %g\n",barrier);
857 21 : log.printf(" using target distribution with BIASFACTOR gamma = %g\n",biasfactor_);
858 21 : if(std::isinf(biasfactor_))
859 4 : log.printf(" (thus a uniform flat target distribution, no well-tempering)\n");
860 21 : if(excluded_region_!=NULL)
861 2 : log.printf(" -- EXCLUDED_REGION: kernels will be deposited only when '%s' is equal to zero\n",excluded_region_->getName().c_str());
862 21 : if(extra_biases_.size()>0)
863 : {
864 2 : log.printf(" -- EXTRA_BIAS: ");
865 5 : for(unsigned e=0; e<extra_biases_.size(); e++)
866 3 : log.printf("%s ",extra_biases_[e]->getName().c_str());
867 2 : log.printf("will be reweighted\n");
868 : }
869 21 : if(adaptive_sigma_)
870 : {
871 11 : log.printf(" adaptive SIGMA will be used, with ADAPTIVE_SIGMA_STRIDE = %u\n",adaptive_sigma_stride_);
872 11 : unsigned x=std::ceil(adaptive_sigma_stride_/stride_);
873 11 : log.printf(" thus the first x kernel depositions will be skipped, x = ADAPTIVE_SIGMA_STRIDE/PACE = %u\n",x);
874 : }
875 : else
876 : {
877 10 : log.printf(" kernels have initial SIGMA = ");
878 29 : for(unsigned i=0; i<ncv_; i++)
879 19 : log.printf(" %g",sigma0_[i]);
880 10 : log.printf("\n");
881 : }
882 21 : if(sigma_min_.size()>0)
883 : {
884 3 : log.printf(" kernels have a SIGMA_MIN = ");
885 9 : for(unsigned i=0; i<ncv_; i++)
886 6 : log.printf(" %g",sigma_min_[i]);
887 3 : log.printf("\n");
888 : }
889 21 : if(fixed_sigma_)
890 6 : log.printf(" -- FIXED_SIGMA: sigma will not decrease as the simulation proceeds\n");
891 21 : log.printf(" kernels are truncated with KERNELS_CUTOFF = %g\n",cutoff);
892 21 : if(cutoff<3.5)
893 0 : log.printf(" +++ WARNING +++ probably kernels are truncated too much\n");
894 21 : log.printf(" the value at cutoff is = %g\n",val_at_cutoff_);
895 21 : log.printf(" regularization EPSILON = %g\n",epsilon_);
896 21 : if(val_at_cutoff_>epsilon_*(1+1e-6))
897 0 : log.printf(" +++ WARNING +++ the KERNEL_CUTOFF might be too small for the given EPSILON\n");
898 21 : log.printf(" kernels will be compressed when closer than COMPRESSION_THRESHOLD = %g\n",sqrt(threshold2_));
899 21 : if(threshold2_==0)
900 0 : log.printf(" +++ WARNING +++ kernels will never merge, expect slowdowns\n");
901 21 : if(!recursive_merge_)
902 6 : log.printf(" -- RECURSIVE_MERGE_OFF: only one merge for each new kernel will be attempted. This is faster only if total number of kernels does not grow too much\n");
903 21 : if(nlist_)
904 5 : log.printf(" -- NLIST: using neighbor list for kernels, with parameters: %g,%g\n",nlist_param_[0],nlist_param_[1]);
905 21 : if(nlist_pace_reset_)
906 2 : log.printf(" -- NLIST_PACE_RESET: forcing the neighbor list to update every PACE\n");
907 21 : if(no_Zed_)
908 6 : log.printf(" -- NO_ZED: using fixed normalization factor = %g\n",Zed_);
909 21 : if(wStateStride_>0)
910 10 : log.printf(" state checkpoints are written on file %s every %d MD steps\n",stateFileName.c_str(),wStateStride_);
911 21 : if(wStateStride_==-1)
912 1 : log.printf(" state checkpoints are written on file %s only on CPT events (or never if MD code does define them!)\n",stateFileName.c_str());
913 21 : if(walkers_mpi)
914 10 : log.printf(" -- WALKERS_MPI: if multiple replicas are present, they will share the same bias via MPI\n");
915 21 : if(NumWalkers_>1)
916 : {
917 10 : log.printf(" using multiple walkers\n");
918 10 : log.printf(" number of walkers: %u\n",NumWalkers_);
919 10 : log.printf(" walker rank: %u\n",walker_rank_);
920 : }
921 21 : int mw_warning=0;
922 21 : if(!walkers_mpi && comm.Get_rank()==0 && multi_sim_comm.Get_size()>(int)NumWalkers_)
923 0 : mw_warning=1;
924 21 : comm.Bcast(mw_warning,0);
925 21 : if(mw_warning) //log.printf messes up with comm, so never use it without Bcast!
926 0 : log.printf(" +++ WARNING +++ multiple replicas will NOT communicate unless the flag WALKERS_MPI is used\n");
927 21 : if(NumParallel_>1)
928 4 : log.printf(" using multiple MPI processes per simulation: %u\n",NumParallel_);
929 21 : if(NumOMP_>1)
930 17 : log.printf(" using multiple OpenMP threads per simulation: %u\n",NumOMP_);
931 21 : if(serial)
932 4 : log.printf(" -- SERIAL: no loop parallelization, despite %d MPI processes and %u OpenMP threads available\n",comm.Get_size(),OpenMP::getNumThreads());
933 21 : log.printf(" Bibliography: ");
934 42 : log<<plumed.cite("M. Invernizzi and M. Parrinello, J. Phys. Chem. Lett. 11, 2731-2736 (2020)");
935 14 : if(mode::explore || adaptive_sigma_)
936 28 : log<<plumed.cite("M. Invernizzi and M. Parrinello, J. Chem. Theory Comput. 18, 3988-3996 (2022)");
937 21 : log.printf("\n");
938 42 : }
939 :
940 : template <class mode>
941 1071 : void OPESmetad<mode>::calculate()
942 : {
943 : //get cv
944 1071 : std::vector<double> cv(ncv_);
945 3111 : for(unsigned i=0; i<ncv_; i++)
946 2040 : cv[i]=getArgument(i);
947 :
948 : //check neighbor list
949 1071 : if(nlist_)
950 : {
951 255 : nlist_steps_++;
952 255 : if(getExchangeStep())
953 0 : nlist_update_=true;
954 : else
955 : {
956 275 : for(unsigned i=0; i<ncv_; i++)
957 : {
958 270 : const double diff_i=difference(i,cv[i],nlist_center_[i]);
959 270 : if(diff_i*diff_i>nlist_param_[1]*nlist_dev2_[i])
960 : {
961 250 : nlist_update_=true;
962 250 : break;
963 : }
964 : }
965 : }
966 255 : if(nlist_update_)
967 250 : updateNlist(cv);
968 : }
969 :
970 : //set bias and forces
971 1071 : std::vector<double> der_prob(ncv_,0);
972 1071 : const double prob=getProbAndDerivatives(cv,der_prob);
973 1071 : const double bias=kbt_*bias_prefactor_*std::log(prob/Zed_+epsilon_);
974 1071 : setBias(bias);
975 3111 : for(unsigned i=0; i<ncv_; i++)
976 2040 : setOutputForce(i,-kbt_*bias_prefactor_/(prob/Zed_+epsilon_)*der_prob[i]/Zed_);
977 1071 : }
978 :
979 : template <class mode>
980 1071 : void OPESmetad<mode>::update()
981 : {
982 1071 : if(isFirstStep_)//same in MetaD, useful for restarts?
983 : {
984 21 : isFirstStep_=false;
985 21 : return;
986 : }
987 :
988 : //update variance if adaptive sigma
989 1050 : if(adaptive_sigma_)
990 : {
991 550 : adaptive_counter_++;
992 550 : unsigned tau=adaptive_sigma_stride_;
993 550 : if(adaptive_counter_<adaptive_sigma_stride_)
994 45 : tau=adaptive_counter_;
995 1600 : for(unsigned i=0; i<ncv_; i++)
996 : { //Welford's online algorithm for standard deviation
997 1050 : const double cv_i=getArgument(i);
998 1050 : const double diff_i=difference(i,av_cv_[i],cv_i);
999 1050 : av_cv_[i]+=diff_i/tau; //exponentially decaying average
1000 1050 : av_M2_[i]+=diff_i*difference(i,av_cv_[i],cv_i);
1001 : }
1002 550 : if(adaptive_counter_<adaptive_sigma_stride_ && counter_==1) //counter_>1 if restarting
1003 : return; //do not apply bias before having measured sigma
1004 : }
1005 :
1006 : //do update
1007 1005 : if(getStep()%stride_==0 && (excluded_region_==NULL || excluded_region_->get()==0))
1008 : {
1009 257 : old_KDEnorm_=KDEnorm_;
1010 257 : delta_kernels_.clear();
1011 257 : unsigned old_nker=kernels_.size();
1012 :
1013 : //get new kernel height
1014 257 : double log_weight=getOutputQuantity(0)/kbt_; //first value is always the current bias
1015 277 : for(unsigned e=0; e<extra_biases_.size(); e++)
1016 20 : log_weight+=extra_biases_[e]->get()/kbt_; //extra biases contribute to the weight
1017 257 : double height=std::exp(log_weight);
1018 :
1019 : //update sum_weights_ and neff
1020 257 : double sum_heights=height;
1021 257 : double sum_heights2=height*height;
1022 257 : if(NumWalkers_>1)
1023 : {
1024 126 : if(comm.Get_rank()==0)
1025 : {
1026 126 : multi_sim_comm.Sum(sum_heights);
1027 126 : multi_sim_comm.Sum(sum_heights2);
1028 : }
1029 126 : comm.Bcast(sum_heights,0);
1030 126 : comm.Bcast(sum_heights2,0);
1031 : }
1032 257 : counter_+=NumWalkers_;
1033 257 : sum_weights_+=sum_heights;
1034 257 : sum_weights2_+=sum_heights2;
1035 257 : const double neff=std::pow(1+sum_weights_,2)/(1+sum_weights2_); //adding 1 makes it more robust at the start
1036 257 : getPntrToComponent("rct")->set(kbt_*std::log(sum_weights_/counter_));
1037 257 : getPntrToComponent("neff")->set(neff);
1038 : if(mode::explore)
1039 : {
1040 68 : KDEnorm_=counter_;
1041 68 : height=1; //plain KDE, bias reweight does not enter here
1042 : }
1043 : else
1044 189 : KDEnorm_=sum_weights_;
1045 :
1046 : //if needed, rescale sigma and height
1047 257 : std::vector<double> sigma=sigma0_;
1048 257 : if(adaptive_sigma_)
1049 : {
1050 93 : const double factor=mode::explore?1:biasfactor_;
1051 131 : if(counter_==1+NumWalkers_) //first time only
1052 : {
1053 14 : for(unsigned i=0; i<ncv_; i++)
1054 9 : av_M2_[i]*=biasfactor_; //from unbiased, thus must be adjusted
1055 14 : for(unsigned i=0; i<ncv_; i++)
1056 9 : sigma0_[i]=std::sqrt(av_M2_[i]/adaptive_counter_/factor);
1057 5 : if(sigma_min_.size()==0)
1058 : {
1059 14 : for(unsigned i=0; i<ncv_; i++)
1060 9 : plumed_massert(sigma0_[i]>1e-6,"ADAPTIVE SIGMA is suspiciously small for CV i="+std::to_string(i)+"\nManually provide SIGMA or set a safe SIGMA_MIN to avoid possible issues");
1061 : }
1062 : else
1063 : {
1064 0 : for(unsigned i=0; i<ncv_; i++)
1065 0 : sigma0_[i]=std::max(sigma0_[i],sigma_min_[i]);
1066 : }
1067 : }
1068 388 : for(unsigned i=0; i<ncv_; i++)
1069 257 : sigma[i]=std::sqrt(av_M2_[i]/adaptive_counter_/factor);
1070 131 : if(sigma_min_.size()==0)
1071 : {
1072 238 : for(unsigned i=0; i<ncv_; i++)
1073 : {
1074 157 : if(sigma[i]<1e-6)
1075 : {
1076 0 : log.printf("+++ WARNING +++ the ADAPTIVE SIGMA is suspiciously small, you should set a safe SIGMA_MIN. 1e-6 will be used here\n");
1077 0 : sigma[i]=1e-6;
1078 0 : sigma_min_.resize(ncv_,1e-6);
1079 : }
1080 : }
1081 : }
1082 : else
1083 : {
1084 150 : for(unsigned i=0; i<ncv_; i++)
1085 100 : sigma[i]=std::max(sigma[i],sigma_min_[i]);
1086 : }
1087 : }
1088 257 : if(!fixed_sigma_)
1089 : {
1090 38 : const double size=mode::explore?counter_:neff; //for EXPLORE neff is not relevant
1091 197 : const double s_rescaling=std::pow(size*(ncv_+2.)/4.,-1./(4+ncv_));
1092 576 : for(unsigned i=0; i<ncv_; i++)
1093 379 : sigma[i]*=s_rescaling;
1094 197 : if(sigma_min_.size()>0)
1095 : {
1096 150 : for(unsigned i=0; i<ncv_; i++)
1097 100 : sigma[i]=std::max(sigma[i],sigma_min_[i]);
1098 : }
1099 : }
1100 : //the height should be divided by sqrt(2*pi)*sigma0_,
1101 : //but this overall factor would be canceled when dividing by Zed
1102 : //thus we skip it altogether, but keep any other sigma rescaling
1103 756 : for(unsigned i=0; i<ncv_; i++)
1104 499 : height*=(sigma0_[i]/sigma[i]);
1105 :
1106 : //get new kernel center
1107 257 : std::vector<double> center(ncv_);
1108 756 : for(unsigned i=0; i<ncv_; i++)
1109 499 : center[i]=getArgument(i);
1110 :
1111 : //add new kernel(s)
1112 257 : if(NumWalkers_==1)
1113 131 : addKernel(height,center,sigma,log_weight);
1114 : else
1115 : {
1116 126 : std::vector<double> all_height(NumWalkers_,0.0);
1117 126 : std::vector<double> all_center(NumWalkers_*ncv_,0.0);
1118 126 : std::vector<double> all_sigma(NumWalkers_*ncv_,0.0);
1119 126 : std::vector<double> all_logweight(NumWalkers_,0.0);
1120 126 : if(comm.Get_rank()==0)
1121 : {
1122 126 : multi_sim_comm.Allgather(height,all_height);
1123 126 : multi_sim_comm.Allgather(center,all_center);
1124 126 : multi_sim_comm.Allgather(sigma,all_sigma);
1125 126 : multi_sim_comm.Allgather(log_weight,all_logweight);
1126 : }
1127 126 : comm.Bcast(all_height,0);
1128 126 : comm.Bcast(all_center,0);
1129 126 : comm.Bcast(all_sigma,0);
1130 126 : comm.Bcast(all_logweight,0);
1131 126 : if(nlist_)
1132 : { //gather all the nlist_index_, so merging can be done using it
1133 50 : std::vector<int> all_nlist_size(NumWalkers_);
1134 50 : if(comm.Get_rank()==0)
1135 : {
1136 50 : all_nlist_size[walker_rank_]=nlist_index_.size();
1137 50 : multi_sim_comm.Sum(all_nlist_size);
1138 : }
1139 50 : comm.Bcast(all_nlist_size,0);
1140 : unsigned tot_size=0;
1141 150 : for(unsigned w=0; w<NumWalkers_; w++)
1142 100 : tot_size+=all_nlist_size[w];
1143 50 : if(tot_size>0)
1144 : {
1145 50 : std::vector<int> disp(NumWalkers_);
1146 100 : for(unsigned w=0; w<NumWalkers_-1; w++)
1147 50 : disp[w+1]=disp[w]+all_nlist_size[w];
1148 50 : std::vector<unsigned> all_nlist_index(tot_size);
1149 50 : if(comm.Get_rank()==0)
1150 50 : multi_sim_comm.Allgatherv(nlist_index_,all_nlist_index,&all_nlist_size[0],&disp[0]);
1151 50 : comm.Bcast(all_nlist_index,0);
1152 50 : std::set<unsigned> nlist_index_set(all_nlist_index.begin(),all_nlist_index.end()); //remove duplicates and sort
1153 50 : nlist_index_.assign(nlist_index_set.begin(),nlist_index_set.end());
1154 : }
1155 : }
1156 378 : for(unsigned w=0; w<NumWalkers_; w++)
1157 : {
1158 252 : std::vector<double> center_w(all_center.begin()+ncv_*w,all_center.begin()+ncv_*(w+1));
1159 252 : std::vector<double> sigma_w(all_sigma.begin()+ncv_*w,all_sigma.begin()+ncv_*(w+1));
1160 252 : addKernel(all_height[w],center_w,sigma_w,all_logweight[w]);
1161 : }
1162 : }
1163 257 : getPntrToComponent("nker")->set(kernels_.size());
1164 257 : if(nlist_)
1165 : {
1166 106 : getPntrToComponent("nlker")->set(nlist_index_.size());
1167 106 : if(nlist_pace_reset_)
1168 50 : nlist_update_=true;
1169 : }
1170 :
1171 : //update Zed_
1172 257 : if(!no_Zed_)
1173 : {
1174 197 : double sum_uprob=0;
1175 197 : const unsigned ks=kernels_.size();
1176 197 : const unsigned ds=delta_kernels_.size();
1177 197 : const bool few_kernels=(ks*ks<(3*ks*ds+2*ds*ds*NumParallel_+100)); //this seems reasonable, but is not rigorous...
1178 197 : if(few_kernels) //really needed? Probably is almost always false
1179 : {
1180 147 : #pragma omp parallel num_threads(NumOMP_)
1181 : {
1182 : #pragma omp for reduction(+:sum_uprob) nowait
1183 : for(unsigned k=rank_; k<kernels_.size(); k+=NumParallel_)
1184 : for(unsigned kk=0; kk<kernels_.size(); kk++)
1185 : sum_uprob+=evaluateKernel(kernels_[kk],kernels_[k].center);
1186 : }
1187 147 : if(NumParallel_>1)
1188 50 : comm.Sum(sum_uprob);
1189 : }
1190 : else
1191 : {
1192 : // Here instead of redoing the full summation, we add only the changes, knowing that
1193 : // uprob = old_uprob + delta_uprob
1194 : // and we also need to consider that in the new sum there are some novel centers and some disappeared ones
1195 50 : double delta_sum_uprob=0;
1196 50 : if(!nlist_)
1197 : {
1198 0 : #pragma omp parallel num_threads(NumOMP_)
1199 : {
1200 : #pragma omp for reduction(+:delta_sum_uprob) nowait
1201 : for(unsigned k=rank_; k<kernels_.size(); k+=NumParallel_)
1202 : {
1203 : for(unsigned d=0; d<delta_kernels_.size(); d++)
1204 : {
1205 : const double sign=delta_kernels_[d].height<0?-1:1; //take away contribution from kernels that are gone, and add the one from new ones
1206 : delta_sum_uprob+=evaluateKernel(delta_kernels_[d],kernels_[k].center)+sign*evaluateKernel(kernels_[k],delta_kernels_[d].center);
1207 : }
1208 : }
1209 : }
1210 : }
1211 : else
1212 : {
1213 50 : #pragma omp parallel num_threads(NumOMP_)
1214 : {
1215 : #pragma omp for reduction(+:delta_sum_uprob) nowait
1216 : for(unsigned nk=rank_; nk<nlist_index_.size(); nk+=NumParallel_)
1217 : {
1218 : const unsigned k=nlist_index_[nk];
1219 : for(unsigned d=0; d<delta_kernels_.size(); d++)
1220 : {
1221 : const double sign=delta_kernels_[d].height<0?-1:1; //take away contribution from kernels that are gone, and add the one from new ones
1222 : delta_sum_uprob+=evaluateKernel(delta_kernels_[d],kernels_[k].center)+sign*evaluateKernel(kernels_[k],delta_kernels_[d].center);
1223 : }
1224 : }
1225 : }
1226 : }
1227 50 : if(NumParallel_>1)
1228 0 : comm.Sum(delta_sum_uprob);
1229 50 : #pragma omp parallel num_threads(NumOMP_)
1230 : {
1231 : #pragma omp for reduction(+:delta_sum_uprob)
1232 : for(unsigned d=0; d<delta_kernels_.size(); d++)
1233 : {
1234 : for(unsigned dd=0; dd<delta_kernels_.size(); dd++)
1235 : { //now subtract the delta_uprob added before, but not needed
1236 : const double sign=delta_kernels_[d].height<0?-1:1;
1237 : delta_sum_uprob-=sign*evaluateKernel(delta_kernels_[dd],delta_kernels_[d].center);
1238 : }
1239 : }
1240 : }
1241 50 : sum_uprob=Zed_*old_KDEnorm_*old_nker+delta_sum_uprob;
1242 : }
1243 197 : Zed_=sum_uprob/KDEnorm_/kernels_.size();
1244 394 : getPntrToComponent("zed")->set(Zed_);
1245 : }
1246 :
1247 : //calculate work if requested
1248 257 : if(calc_work_)
1249 : {
1250 70 : std::vector<double> dummy(ncv_); //derivatives are not actually needed
1251 70 : const double prob=getProbAndDerivatives(center,dummy);
1252 70 : const double new_bias=kbt_*bias_prefactor_*std::log(prob/Zed_+epsilon_);
1253 70 : work_+=new_bias-getOutputQuantity(0);
1254 140 : getPntrToComponent("work")->set(work_);
1255 : }
1256 : }
1257 :
1258 : //dump state if requested
1259 1005 : if( (wStateStride_>0 && getStep()%wStateStride_==0) || (wStateStride_==-1 && getCPT()) )
1260 18 : dumpStateToFile();
1261 : }
1262 :
1263 : template <class mode>
1264 1141 : double OPESmetad<mode>::getProbAndDerivatives(const std::vector<double>& cv,std::vector<double>& der_prob)
1265 : {
1266 1141 : double prob=0.0;
1267 1141 : if(!nlist_)
1268 : {
1269 886 : if(NumOMP_==1 || (unsigned)kernels_.size()<2*NumOMP_*NumParallel_)
1270 : {
1271 : // for performances and thread safety
1272 707 : std::vector<double> dist(ncv_);
1273 2647 : for(unsigned k=rank_; k<kernels_.size(); k+=NumParallel_)
1274 1940 : prob+=evaluateKernel(kernels_[k],cv,der_prob,dist);
1275 : }
1276 : else
1277 : {
1278 179 : #pragma omp parallel num_threads(NumOMP_)
1279 : {
1280 : std::vector<double> omp_deriv(der_prob.size(),0.);
1281 : // for performances and thread safety
1282 : std::vector<double> dist(ncv_);
1283 : #pragma omp for reduction(+:prob) nowait
1284 : for(unsigned k=rank_; k<kernels_.size(); k+=NumParallel_)
1285 : prob+=evaluateKernel(kernels_[k],cv,omp_deriv,dist);
1286 : #pragma omp critical
1287 : for(unsigned i=0; i<ncv_; i++)
1288 : der_prob[i]+=omp_deriv[i];
1289 : }
1290 : }
1291 : }
1292 : else
1293 : {
1294 255 : if(NumOMP_==1 || (unsigned)nlist_index_.size()<2*NumOMP_*NumParallel_)
1295 : {
1296 : // for performances and thread safety
1297 230 : std::vector<double> dist(ncv_);
1298 660 : for(unsigned nk=rank_; nk<nlist_index_.size(); nk+=NumParallel_)
1299 430 : prob+=evaluateKernel(kernels_[nlist_index_[nk]],cv,der_prob,dist);
1300 : }
1301 : else
1302 : {
1303 25 : #pragma omp parallel num_threads(NumOMP_)
1304 : {
1305 : std::vector<double> omp_deriv(der_prob.size(),0.);
1306 : // for performances and thread safety
1307 : std::vector<double> dist(ncv_);
1308 : #pragma omp for reduction(+:prob) nowait
1309 : for(unsigned nk=rank_; nk<nlist_index_.size(); nk+=NumParallel_)
1310 : prob+=evaluateKernel(kernels_[nlist_index_[nk]],cv,omp_deriv,dist);
1311 : #pragma omp critical
1312 : for(unsigned i=0; i<ncv_; i++)
1313 : der_prob[i]+=omp_deriv[i];
1314 : }
1315 : }
1316 : }
1317 1141 : if(NumParallel_>1)
1318 : {
1319 224 : comm.Sum(prob);
1320 224 : comm.Sum(der_prob);
1321 : }
1322 : //normalize the estimate
1323 1141 : prob/=KDEnorm_;
1324 3311 : for(unsigned i=0; i<ncv_; i++)
1325 2170 : der_prob[i]/=KDEnorm_;
1326 :
1327 1141 : return prob;
1328 : }
1329 :
1330 : template <class mode>
1331 513 : void OPESmetad<mode>::addKernel(const double height,const std::vector<double>& center,const std::vector<double>& sigma)
1332 : {
1333 : bool no_match=true;
1334 513 : if(threshold2_!=0)
1335 : {
1336 513 : unsigned taker_k=getMergeableKernel(center,kernels_.size());
1337 513 : if(taker_k<kernels_.size())
1338 : {
1339 : no_match=false;
1340 388 : delta_kernels_.emplace_back(-1*kernels_[taker_k].height,kernels_[taker_k].center,kernels_[taker_k].sigma);
1341 776 : mergeKernels(kernels_[taker_k],kernel(height,center,sigma));
1342 388 : delta_kernels_.push_back(kernels_[taker_k]);
1343 388 : if(recursive_merge_) //the overhead is worth it if it keeps low the total number of kernels
1344 : {
1345 : unsigned giver_k=taker_k;
1346 337 : taker_k=getMergeableKernel(kernels_[giver_k].center,giver_k);
1347 337 : while(taker_k<kernels_.size())
1348 : {
1349 0 : delta_kernels_.pop_back();
1350 0 : delta_kernels_.emplace_back(-1*kernels_[taker_k].height,kernels_[taker_k].center,kernels_[taker_k].sigma);
1351 0 : if(taker_k>giver_k) //saves time when erasing
1352 : std::swap(taker_k,giver_k);
1353 0 : mergeKernels(kernels_[taker_k],kernels_[giver_k]);
1354 0 : delta_kernels_.push_back(kernels_[taker_k]);
1355 0 : kernels_.erase(kernels_.begin()+giver_k);
1356 0 : if(nlist_)
1357 : {
1358 : unsigned giver_nk=0;
1359 : bool found_giver=false;
1360 0 : for(unsigned nk=0; nk<nlist_index_.size(); nk++)
1361 : {
1362 0 : if(found_giver)
1363 0 : nlist_index_[nk]--; //all the indexes shift due to erase
1364 0 : if(nlist_index_[nk]==giver_k)
1365 : {
1366 : giver_nk=nk;
1367 : found_giver=true;
1368 : }
1369 : }
1370 : plumed_dbg_massert(found_giver,"problem with merging and NLIST");
1371 0 : nlist_index_.erase(nlist_index_.begin()+giver_nk);
1372 : }
1373 : giver_k=taker_k;
1374 0 : taker_k=getMergeableKernel(kernels_[giver_k].center,giver_k);
1375 : }
1376 : }
1377 : }
1378 : }
1379 : if(no_match)
1380 : {
1381 125 : kernels_.emplace_back(height,center,sigma);
1382 125 : delta_kernels_.emplace_back(height,center,sigma);
1383 125 : if(nlist_)
1384 12 : nlist_index_.push_back(kernels_.size()-1);
1385 : }
1386 513 : }
1387 :
1388 : template <class mode>
1389 383 : void OPESmetad<mode>::addKernel(const double height,const std::vector<double>& center,const std::vector<double>& sigma,const double logweight)
1390 : {
1391 383 : addKernel(height,center,sigma);
1392 : //write to file
1393 383 : kernelsOfile_.printField("time",getTime());
1394 1134 : for(unsigned i=0; i<ncv_; i++)
1395 751 : kernelsOfile_.printField(getPntrToArgument(i),center[i]);
1396 1134 : for(unsigned i=0; i<ncv_; i++)
1397 1502 : kernelsOfile_.printField("sigma_"+getPntrToArgument(i)->getName(),sigma[i]);
1398 383 : kernelsOfile_.printField("height",height);
1399 383 : kernelsOfile_.printField("logweight",logweight);
1400 383 : kernelsOfile_.printField();
1401 383 : }
1402 :
1403 : template <class mode>
1404 850 : unsigned OPESmetad<mode>::getMergeableKernel(const std::vector<double>& giver_center,const unsigned giver_k)
1405 : { //returns kernels_.size() if no match is found
1406 850 : unsigned min_k=kernels_.size();
1407 850 : double min_norm2=threshold2_;
1408 850 : if(!nlist_)
1409 : {
1410 550 : #pragma omp parallel num_threads(NumOMP_)
1411 : {
1412 : unsigned min_k_omp = min_k;
1413 : double min_norm2_omp = threshold2_;
1414 : #pragma omp for nowait
1415 : for(unsigned k=rank_; k<kernels_.size(); k+=NumParallel_)
1416 : {
1417 : if(k==giver_k) //a kernel should not be merged with itself
1418 : continue;
1419 : double norm2=0;
1420 : for(unsigned i=0; i<ncv_; i++)
1421 : {
1422 : const double dist_i=difference(i,giver_center[i],kernels_[k].center[i])/kernels_[k].sigma[i];
1423 : norm2+=dist_i*dist_i;
1424 : if(norm2>=min_norm2_omp)
1425 : break;
1426 : }
1427 : if(norm2<min_norm2_omp)
1428 : {
1429 : min_norm2_omp=norm2;
1430 : min_k_omp=k;
1431 : }
1432 : }
1433 : #pragma omp critical
1434 : {
1435 : if(min_norm2_omp < min_norm2)
1436 : {
1437 : min_norm2 = min_norm2_omp;
1438 : min_k = min_k_omp;
1439 : }
1440 : }
1441 : }
1442 : }
1443 : else
1444 : {
1445 300 : #pragma omp parallel num_threads(NumOMP_)
1446 : {
1447 : unsigned min_k_omp = min_k;
1448 : double min_norm2_omp = threshold2_;
1449 : #pragma omp for nowait
1450 : for(unsigned nk=rank_; nk<nlist_index_.size(); nk+=NumParallel_)
1451 : {
1452 : const unsigned k=nlist_index_[nk];
1453 : if(k==giver_k) //a kernel should not be merged with itself
1454 : continue;
1455 : double norm2=0;
1456 : for(unsigned i=0; i<ncv_; i++)
1457 : {
1458 : const double dist_i=difference(i,giver_center[i],kernels_[k].center[i])/kernels_[k].sigma[i];
1459 : norm2+=dist_i*dist_i;
1460 : if(norm2>=min_norm2)
1461 : break;
1462 : }
1463 : if(norm2<min_norm2_omp)
1464 : {
1465 : min_norm2_omp=norm2;
1466 : min_k_omp=k;
1467 : }
1468 : }
1469 : #pragma omp critical
1470 : {
1471 : if(min_norm2_omp < min_norm2)
1472 : {
1473 : min_norm2 = min_norm2_omp;
1474 : min_k = min_k_omp;
1475 : }
1476 : }
1477 : }
1478 : }
1479 850 : if(NumParallel_>1)
1480 : {
1481 134 : std::vector<double> all_min_norm2(NumParallel_);
1482 134 : std::vector<unsigned> all_min_k(NumParallel_);
1483 134 : comm.Allgather(min_norm2,all_min_norm2);
1484 134 : comm.Allgather(min_k,all_min_k);
1485 : const unsigned best=std::distance(std::begin(all_min_norm2),std::min_element(std::begin(all_min_norm2),std::end(all_min_norm2)));
1486 134 : min_k=all_min_k[best];
1487 : }
1488 850 : return min_k;
1489 : }
1490 :
1491 : template <class mode>
1492 250 : void OPESmetad<mode>::updateNlist(const std::vector<double>& new_center)
1493 : {
1494 250 : if(kernels_.size()==0) //no need to check for neighbors
1495 6 : return;
1496 :
1497 244 : nlist_center_=new_center;
1498 244 : nlist_index_.clear();
1499 : //first we gather all the nlist_index
1500 244 : if(NumOMP_==1 || (unsigned)kernels_.size()<2*NumOMP_*NumParallel_)
1501 : {
1502 198 : for(unsigned k=rank_; k<kernels_.size(); k+=NumParallel_)
1503 : {
1504 : double norm2_k=0;
1505 444 : for(unsigned i=0; i<ncv_; i++)
1506 : {
1507 296 : const double dist_ik=difference(i,nlist_center_[i],kernels_[k].center[i])/kernels_[k].sigma[i];
1508 296 : norm2_k+=dist_ik*dist_ik;
1509 : }
1510 148 : if(norm2_k<=nlist_param_[0]*cutoff2_)
1511 104 : nlist_index_.push_back(k);
1512 : }
1513 : }
1514 : else
1515 : {
1516 194 : #pragma omp parallel num_threads(NumOMP_)
1517 : {
1518 : std::vector<unsigned> private_nlist_index;
1519 : #pragma omp for nowait
1520 : for(unsigned k=rank_; k<kernels_.size(); k+=NumParallel_)
1521 : {
1522 : double norm2_k=0;
1523 : for(unsigned i=0; i<ncv_; i++)
1524 : {
1525 : const double dist_ik=difference(i,nlist_center_[i],kernels_[k].center[i])/kernels_[k].sigma[i];
1526 : norm2_k+=dist_ik*dist_ik;
1527 : }
1528 : if(norm2_k<=nlist_param_[0]*cutoff2_)
1529 : private_nlist_index.push_back(k);
1530 : }
1531 : #pragma omp critical
1532 : nlist_index_.insert(nlist_index_.end(),private_nlist_index.begin(),private_nlist_index.end());
1533 : }
1534 194 : if(recursive_merge_)
1535 194 : std::sort(nlist_index_.begin(),nlist_index_.end());
1536 : }
1537 244 : if(NumParallel_>1)
1538 : {
1539 100 : std::vector<int> all_nlist_size(NumParallel_);
1540 100 : all_nlist_size[rank_]=nlist_index_.size();
1541 100 : comm.Sum(all_nlist_size);
1542 : unsigned tot_size=0;
1543 300 : for(unsigned r=0; r<NumParallel_; r++)
1544 200 : tot_size+=all_nlist_size[r];
1545 100 : if(tot_size>0)
1546 : {
1547 100 : std::vector<int> disp(NumParallel_);
1548 200 : for(unsigned r=0; r<NumParallel_-1; r++)
1549 100 : disp[r+1]=disp[r]+all_nlist_size[r];
1550 100 : std::vector<unsigned> local_nlist_index=nlist_index_;
1551 100 : nlist_index_.resize(tot_size);
1552 100 : comm.Allgatherv(local_nlist_index,nlist_index_,&all_nlist_size[0],&disp[0]);
1553 100 : if(recursive_merge_)
1554 100 : std::sort(nlist_index_.begin(),nlist_index_.end());
1555 : }
1556 : }
1557 : //calculate the square deviation
1558 244 : std::vector<double> dev2(ncv_,0.);
1559 773 : for(unsigned k=rank_; k<nlist_index_.size(); k+=NumParallel_)
1560 : {
1561 1587 : for(unsigned i=0; i<ncv_; i++)
1562 : {
1563 1058 : const double diff_ik=difference(i,nlist_center_[i],kernels_[nlist_index_[k]].center[i]);
1564 1058 : dev2[i]+=diff_ik*diff_ik;
1565 : }
1566 : }
1567 244 : if(NumParallel_>1)
1568 100 : comm.Sum(dev2);
1569 732 : for(unsigned i=0; i<ncv_; i++)
1570 : {
1571 488 : if(dev2[i]==0) //e.g. if nlist_index_.size()==0
1572 56 : nlist_dev2_[i]=std::pow(kernels_.back().sigma[i],2);
1573 : else
1574 432 : nlist_dev2_[i]=dev2[i]/nlist_index_.size();
1575 : }
1576 244 : getPntrToComponent("nlker")->set(nlist_index_.size());
1577 244 : getPntrToComponent("nlsteps")->set(nlist_steps_);
1578 244 : nlist_steps_=0;
1579 244 : nlist_update_=false;
1580 : }
1581 :
1582 : template <class mode>
1583 18 : void OPESmetad<mode>::dumpStateToFile()
1584 : {
1585 : //gather adaptive sigma info if needed
1586 : //doing this while writing to file can lead to misterious slowdowns
1587 : std::vector<double> all_sigma0;
1588 : std::vector<double> all_av_cv;
1589 : std::vector<double> all_av_M2;
1590 18 : if(adaptive_sigma_ && NumWalkers_>1)
1591 : {
1592 16 : all_sigma0.resize(NumWalkers_*ncv_);
1593 16 : all_av_cv.resize(NumWalkers_*ncv_);
1594 16 : all_av_M2.resize(NumWalkers_*ncv_);
1595 16 : if(comm.Get_rank()==0)
1596 : {
1597 16 : multi_sim_comm.Allgather(sigma0_,all_sigma0);
1598 16 : multi_sim_comm.Allgather(av_cv_,all_av_cv);
1599 16 : multi_sim_comm.Allgather(av_M2_,all_av_M2);
1600 : }
1601 16 : comm.Bcast(all_sigma0,0);
1602 16 : comm.Bcast(all_av_cv,0);
1603 16 : comm.Bcast(all_av_M2,0);
1604 : }
1605 :
1606 : //rewrite header or rewind file
1607 18 : if(storeOldStates_)
1608 16 : stateOfile_.clearFields();
1609 2 : else if(walker_rank_==0)
1610 2 : stateOfile_.rewind();
1611 : //define fields
1612 18 : stateOfile_.addConstantField("action");
1613 18 : stateOfile_.addConstantField("biasfactor");
1614 18 : stateOfile_.addConstantField("epsilon");
1615 18 : stateOfile_.addConstantField("kernel_cutoff");
1616 18 : stateOfile_.addConstantField("compression_threshold");
1617 18 : stateOfile_.addConstantField("zed");
1618 18 : stateOfile_.addConstantField("sum_weights");
1619 18 : stateOfile_.addConstantField("sum_weights2");
1620 18 : stateOfile_.addConstantField("counter");
1621 18 : if(adaptive_sigma_)
1622 : {
1623 18 : stateOfile_.addConstantField("adaptive_counter");
1624 18 : if(NumWalkers_==1)
1625 : {
1626 6 : for(unsigned i=0; i<ncv_; i++)
1627 : {
1628 8 : stateOfile_.addConstantField("sigma0_"+getPntrToArgument(i)->getName());
1629 8 : stateOfile_.addConstantField("av_cv_"+getPntrToArgument(i)->getName());
1630 8 : stateOfile_.addConstantField("av_M2_"+getPntrToArgument(i)->getName());
1631 : }
1632 : }
1633 : else
1634 : {
1635 48 : for(unsigned w=0; w<NumWalkers_; w++)
1636 96 : for(unsigned i=0; i<ncv_; i++)
1637 : {
1638 128 : const std::string arg_iw=getPntrToArgument(i)->getName()+"_"+std::to_string(w);
1639 64 : stateOfile_.addConstantField("sigma0_"+arg_iw);
1640 64 : stateOfile_.addConstantField("av_cv_"+arg_iw);
1641 128 : stateOfile_.addConstantField("av_M2_"+arg_iw);
1642 : }
1643 : }
1644 : }
1645 : //print fields
1646 54 : for(unsigned i=0; i<ncv_; i++) //periodicity of CVs
1647 36 : stateOfile_.setupPrintValue(getPntrToArgument(i));
1648 36 : stateOfile_.printField("action",getName()+"_state");
1649 18 : stateOfile_.printField("biasfactor",biasfactor_);
1650 18 : stateOfile_.printField("epsilon",epsilon_);
1651 18 : stateOfile_.printField("kernel_cutoff",sqrt(cutoff2_));
1652 18 : stateOfile_.printField("compression_threshold",sqrt(threshold2_));
1653 18 : stateOfile_.printField("zed",Zed_);
1654 18 : stateOfile_.printField("sum_weights",sum_weights_);
1655 18 : stateOfile_.printField("sum_weights2",sum_weights2_);
1656 18 : stateOfile_.printField("counter",counter_);
1657 18 : if(adaptive_sigma_)
1658 : {
1659 18 : stateOfile_.printField("adaptive_counter",adaptive_counter_);
1660 18 : if(NumWalkers_==1)
1661 : {
1662 6 : for(unsigned i=0; i<ncv_; i++)
1663 : {
1664 8 : stateOfile_.printField("sigma0_"+getPntrToArgument(i)->getName(),sigma0_[i]);
1665 8 : stateOfile_.printField("av_cv_"+getPntrToArgument(i)->getName(),av_cv_[i]);
1666 8 : stateOfile_.printField("av_M2_"+getPntrToArgument(i)->getName(),av_M2_[i]);
1667 : }
1668 : }
1669 : else
1670 : {
1671 48 : for(unsigned w=0; w<NumWalkers_; w++)
1672 96 : for(unsigned i=0; i<ncv_; i++)
1673 : {
1674 128 : const std::string arg_iw=getPntrToArgument(i)->getName()+"_"+std::to_string(w);
1675 64 : stateOfile_.printField("sigma0_"+arg_iw,all_sigma0[w*ncv_+i]);
1676 64 : stateOfile_.printField("av_cv_"+arg_iw,all_av_cv[w*ncv_+i]);
1677 128 : stateOfile_.printField("av_M2_"+arg_iw,all_av_M2[w*ncv_+i]);
1678 : }
1679 : }
1680 : }
1681 : //print kernels
1682 82 : for(unsigned k=0; k<kernels_.size(); k++)
1683 : {
1684 64 : stateOfile_.printField("time",getTime()); //this is not very usefull
1685 192 : for(unsigned i=0; i<ncv_; i++)
1686 128 : stateOfile_.printField(getPntrToArgument(i),kernels_[k].center[i]);
1687 192 : for(unsigned i=0; i<ncv_; i++)
1688 256 : stateOfile_.printField("sigma_"+getPntrToArgument(i)->getName(),kernels_[k].sigma[i]);
1689 64 : stateOfile_.printField("height",kernels_[k].height);
1690 64 : stateOfile_.printField();
1691 : }
1692 : //make sure file is written even if small
1693 18 : if(!storeOldStates_)
1694 2 : stateOfile_.flush();
1695 18 : }
1696 :
1697 : template <class mode>
1698 5969 : inline double OPESmetad<mode>::evaluateKernel(const kernel& G,const std::vector<double>& x) const
1699 : { //NB: cannot be a method of kernel class, because uses external variables (for cutoff)
1700 : double norm2=0;
1701 13923 : for(unsigned i=0; i<ncv_; i++)
1702 : {
1703 10288 : const double dist_i=difference(i,G.center[i],x[i])/G.sigma[i];
1704 10288 : norm2+=dist_i*dist_i;
1705 10288 : if(norm2>=cutoff2_)
1706 : return 0;
1707 : }
1708 3635 : return G.height*(std::exp(-0.5*norm2)-val_at_cutoff_);
1709 : }
1710 :
1711 : template <class mode>
1712 3364 : inline double OPESmetad<mode>::evaluateKernel(const kernel& G,const std::vector<double>& x, std::vector<double>& acc_der, std::vector<double>& dist)
1713 : { //NB: cannot be a method of kernel class, because uses external variables (for cutoff)
1714 : double norm2=0;
1715 7483 : for(unsigned i=0; i<ncv_; i++)
1716 : {
1717 5578 : dist[i]=difference(i,G.center[i],x[i])/G.sigma[i];
1718 5578 : norm2+=dist[i]*dist[i];
1719 5578 : if(norm2>=cutoff2_)
1720 : return 0;
1721 : }
1722 1905 : const double val=G.height*(std::exp(-0.5*norm2)-val_at_cutoff_);
1723 5576 : for(unsigned i=0; i<ncv_; i++)
1724 3671 : acc_der[i]-=dist[i]/G.sigma[i]*val; //NB: we accumulate the derivative into der
1725 : return val;
1726 : }
1727 :
1728 : template <class mode>
1729 388 : inline void OPESmetad<mode>::mergeKernels(kernel& k1,const kernel& k2)
1730 : {
1731 388 : const double h=k1.height+k2.height;
1732 1159 : for(unsigned i=0; i<k1.center.size(); i++)
1733 : {
1734 771 : const bool isPeriodic_i=getPntrToArgument(i)->isPeriodic();
1735 771 : if(isPeriodic_i)
1736 771 : k1.center[i]=k2.center[i]+difference(i,k2.center[i],k1.center[i]); //fix PBC
1737 771 : const double c_i=(k1.height*k1.center[i]+k2.height*k2.center[i])/h;
1738 771 : const double ss_k1_part=k1.height*(k1.sigma[i]*k1.sigma[i]+k1.center[i]*k1.center[i]);
1739 771 : const double ss_k2_part=k2.height*(k2.sigma[i]*k2.sigma[i]+k2.center[i]*k2.center[i]);
1740 771 : const double ss_i=(ss_k1_part+ss_k2_part)/h-c_i*c_i;
1741 771 : if(isPeriodic_i)
1742 771 : k1.center[i]=bringBackInPbc(i,c_i);
1743 : else
1744 0 : k1.center[i]=c_i;
1745 771 : k1.sigma[i]=sqrt(ss_i);
1746 : }
1747 388 : k1.height=h;
1748 388 : }
1749 :
1750 : }
1751 : }
|
