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