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