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