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