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 : To continue a simulation with a static bias, for example once the bias has converged, restart from the STATE file and add UPDATE_UNTIL=0.
86 : No new kernels are then deposited, and the bias remains the one read from STATE_RFILE (no STATE_WFILE is needed).
87 :
88 : Multiple walkers are supported only with MPI communication, via the keyword WALKERS_MPI.
89 :
90 : \par Examples
91 :
92 : 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.
93 : The \ref opes-metad can also be useful to get started with the method.
94 :
95 : The following is a minimal working example:
96 :
97 : \plumedfile
98 : cv: DISTANCE ATOMS=1,2
99 : opes: OPES_METAD ARG=cv PACE=200 BARRIER=40
100 : PRINT STRIDE=200 FILE=COLVAR ARG=*
101 : \endplumedfile
102 :
103 : Another more articulated one:
104 :
105 : \plumedfile
106 : phi: TORSION ATOMS=5,7,9,15
107 : psi: TORSION ATOMS=7,9,15,17
108 : opes: OPES_METAD ...
109 : FILE=Kernels.data
110 : TEMP=300
111 : ARG=phi,psi
112 : PACE=500
113 : BARRIER=50
114 : SIGMA=0.15,0.15
115 : SIGMA_MIN=0.01,0.01
116 : STATE_RFILE=Restart.data
117 : STATE_WFILE=State.data
118 : STATE_WSTRIDE=500*100
119 : STORE_STATES
120 : WALKERS_MPI
121 : NLIST
122 : ...
123 : PRINT FMT=%g STRIDE=500 FILE=Colvar.data ARG=phi,psi,opes.*
124 : \endplumedfile
125 :
126 : Next is an example of how to define a custom target distribution different from the well-tempered one.
127 : Here we chose to focus more on the transition state, that is around \f$\phi=0\f$.
128 : 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$.
129 :
130 : \plumedfile
131 : phi: TORSION ATOMS=5,7,9,15
132 : FtgValue: CUSTOM ARG=phi PERIODIC=NO FUNC=(x/0.4)^2
133 : Ftg: BIASVALUE ARG=FtgValue
134 : opes: OPES_METAD ...
135 : ARG=phi
136 : PACE=500
137 : BARRIER=50
138 : SIGMA=0.2
139 : BIASFACTOR=inf
140 : EXTRA_BIAS=Ftg.bias
141 : ...
142 : PRINT FMT=%g STRIDE=500 FILE=COLVAR ARG=phi,Ftg.bias,opes.bias
143 : \endplumedfile
144 :
145 : 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.
146 :
147 : Finally, an example of how to use the EXCLUDED_REGION keyword.
148 : It expects a characteristic function that is different from zero in the region to be excluded.
149 : You can use \ref CUSTOM and a combination of the step function to define it.
150 : 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$:
151 :
152 : \plumedfile
153 : phi: TORSION ATOMS=5,7,9,15
154 : psi: TORSION ATOMS=7,9,15,17
155 : xx: CUSTOM PERIODIC=NO ARG=phi FUNC=step(x+0.6)-step(x-0.7)
156 : opes: OPES_METAD ...
157 : ARG=phi,psi
158 : PACE=500
159 : BARRIER=30
160 : EXCLUDED_REGION=xx
161 : NLIST
162 : ...
163 : PRINT FMT=%g STRIDE=500 FILE=COLVAR ARG=phi,psi,xx,opes.*
164 : \endplumedfile
165 :
166 : */
167 : //+ENDPLUMEDOC
168 :
169 : template <class mode>
170 : class OPESmetad : public bias::Bias {
171 :
172 : private:
173 : bool isFirstStep_;
174 : unsigned NumOMP_;
175 : unsigned NumParallel_;
176 : unsigned rank_;
177 : unsigned NumWalkers_;
178 : unsigned walker_rank_;
179 : unsigned long long counter_;
180 : std::size_t ncv_;
181 :
182 : double kbt_;
183 : double biasfactor_;
184 : double bias_prefactor_;
185 : unsigned stride_;
186 : std::vector<double> sigma0_;
187 : std::vector<double> sigma_min_;
188 : unsigned adaptive_sigma_stride_;
189 : unsigned long long adaptive_counter_;
190 : std::vector<double> av_cv_;
191 : std::vector<double> av_M2_;
192 : bool fixed_sigma_;
193 : bool adaptive_sigma_;
194 : double epsilon_;
195 : double sum_weights_;
196 : double sum_weights2_;
197 :
198 : bool no_Zed_;
199 : double Zed_;
200 : double KDEnorm_;
201 :
202 : double threshold2_;
203 : bool recursive_merge_;
204 : //kernels are truncated diagonal Gaussians
205 26734 : struct kernel {
206 : double height;
207 : std::vector<double> center;
208 : std::vector<double> sigma;
209 27111 : kernel(double h, const std::vector<double>& c,const std::vector<double>& s):
210 27111 : height(h),center(c),sigma(s) {}
211 : };
212 : double cutoff2_;
213 : double val_at_cutoff_;
214 : void mergeKernels(kernel&,const kernel&); //merge the second one into the first one
215 : double evaluateKernel(const kernel&,const std::vector<double>&) const;
216 : double shiftedGaussian(const double) const; //for a squared distance within the cutoff
217 : std::vector<kernel> kernels_; //all compressed kernels
218 : //kernels summed at every step (all of them, or those in the neighbor list), with each CV in a
219 : //contiguous array, so that the distances to all of them are computed in a loop without branches
220 : bool sum_outdated_; //kernels_ or nlist_index_ changed since the last updateSum()
221 : std::vector<double> sum_center_; //CV i of the k-th summed kernel is sum_center_[i*size+k]
222 : std::vector<double> sum_inv_sigma_; //same layout, inverse of sigma
223 : std::vector<double> sum_height_;
224 : std::vector<double> period_; //period of each CV, zero if not periodic
225 : std::vector<double> inv_period_;
226 : std::vector<double> cv_; //buffers for calculate(), to avoid allocations at every step
227 : std::vector<double> der_prob_;
228 : void updateSum();
229 : void sumKernels(const std::vector<double>&,const unsigned,const unsigned,double&,std::vector<double>&) const;
230 : OFile kernelsOfile_;
231 : //neighbour list stuff
232 : bool nlist_;
233 : double nlist_param_; //kernels within sqrt(nlist_param_)*cutoff are in the list
234 : std::vector<unsigned> nlist_index_;
235 : std::vector<double> nlist_center_;
236 : unsigned nlist_steps_;
237 : bool nlist_update_;
238 : bool nlist_pace_reset_;
239 : double nlist_skin2_; //square of the largest displacement allowed, in units of nlist_sigma_min_
240 : std::vector<double> nlist_sigma_min_; //smallest sigma of all kernels, for each CV, at the last rebuild
241 : //all kernels, with the same layout as the summed ones, to rebuild the list and update Zed quickly
242 : bool nlist_all_outdated_; //kernels_ changed since the last updateNlistAll()
243 : std::vector<double> nlist_all_center_;
244 : std::vector<double> nlist_all_inv_sigma_;
245 : std::vector<double> nlist_all_height_;
246 : std::vector<double> nlist_all_sigma_min_; //smallest sigma of all kernels, for each CV
247 : void updateNlistAll();
248 : void listKernels(const unsigned,const unsigned,std::vector<unsigned>&) const;
249 :
250 : bool calc_work_;
251 : double work_;
252 : double old_KDEnorm_;
253 : std::vector<kernel> delta_kernels_;
254 : double deltaSumUprob(const unsigned,const unsigned,const bool) const;
255 :
256 : Value* excluded_region_;
257 : std::vector<Value*> extra_biases_;
258 :
259 : OFile stateOfile_;
260 : int wStateStride_;
261 : bool storeOldStates_;
262 :
263 : double getProbAndDerivatives(const std::vector<double>&,std::vector<double>&);
264 : void addKernel(const double,const std::vector<double>&,const std::vector<double>&);
265 : void addKernel(const double,const std::vector<double>&,const std::vector<double>&,const double); //also print to file
266 : unsigned getMergeableKernel(const std::vector<double>&,const unsigned,const bool all_kernels=false);
267 : void updateNlist(const std::vector<double>&);
268 : void dumpStateToFile();
269 :
270 : public:
271 : explicit OPESmetad(const ActionOptions&);
272 : void calculate() override;
273 : void update() override;
274 : static void registerKeywords(Keywords& keys);
275 : };
276 :
277 : struct convergence {
278 : static const bool explore=false;
279 : };
280 : typedef OPESmetad<convergence> OPESmetad_c;
281 : PLUMED_REGISTER_ACTION(OPESmetad_c,"OPES_METAD")
282 :
283 : //OPES_METAD_EXPLORE is very similar from the point of view of the code,
284 : //but conceptually it is better to make it a separate BIAS action
285 :
286 : //+PLUMEDOC OPES_BIAS OPES_METAD_EXPLORE
287 : /*
288 : On-the-fly probability enhanced sampling with well-tempered target distribution in exploreation mode.
289 :
290 : On-the-fly probability enhanced sampling with well-tempered target distribution (\ref OPES "OPES") with well-tempered target distribution, exploration mode \cite Invernizzi2022explore .
291 :
292 : 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$.
293 : 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.
294 : The bias at step \f$n\f$ is
295 : \f[
296 : V_n(\mathbf{s}) = (\gamma-1)\frac{1}{\beta}\log\left(\frac{p^{\text{WT}}_n(\mathbf{s})}{Z_n}+\epsilon\right)\, .
297 : \f]
298 : See Ref.\cite Invernizzi2022explore for a complete description of the method.
299 :
300 : 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.
301 : Given enough simulation time, both will converge to the same bias potential but they do so in a qualitatively different way.
302 : 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.
303 : This goes at the expenses of a typically slower convergence of the reweight estimate.
304 : \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.
305 :
306 : Similarly to \ref OPES_METAD, also \ref OPES_METAD_EXPLORE uses a kernel density estimation with the same on-the-fly compression algorithm.
307 : The only difference is that the kernels are not weighted, since it estimates the sampled distribution and not the reweighted unbiased one.
308 :
309 : 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.
310 :
311 : \par Examples
312 :
313 : The following is a minimal working example:
314 :
315 : \plumedfile
316 : cv: DISTANCE ATOMS=1,2
317 : opes: OPES_METAD_EXPLORE ARG=cv PACE=500 BARRIER=40
318 : PRINT STRIDE=100 FILE=COLVAR ARG=cv,opes.*
319 : \endplumedfile
320 : */
321 : //+ENDPLUMEDOC
322 :
323 : struct exploration {
324 : static const bool explore=true;
325 : };
326 : typedef OPESmetad<exploration> OPESmetad_e;
327 : PLUMED_REGISTER_ACTION(OPESmetad_e,"OPES_METAD_EXPLORE")
328 :
329 : template <class mode>
330 61 : void OPESmetad<mode>::registerKeywords(Keywords& keys) {
331 61 : Bias::registerKeywords(keys);
332 61 : keys.use("ARG");
333 122 : keys.add("compulsory","TEMP","-1","temperature. If not set, it is taken from MD engine, but not all MD codes provide it");
334 122 : keys.add("compulsory","PACE","the frequency for kernel deposition");
335 61 : std::string info_sigma("the initial widths of the kernels");
336 : if(mode::explore) {
337 : info_sigma+=", divided by the square root of gamma";
338 : }
339 : info_sigma+=". If not set, an adaptive sigma will be used with the given ADAPTIVE_SIGMA_STRIDE";
340 122 : keys.add("compulsory","SIGMA","ADAPTIVE",info_sigma);
341 122 : keys.add("compulsory","BARRIER","the free energy barrier to be overcome. It is used to set BIASFACTOR, EPSILON, and KERNEL_CUTOFF to reasonable values");
342 122 : keys.add("compulsory","COMPRESSION_THRESHOLD","1","merge kernels if closer than this threshold, in units of sigma");
343 : //extra options
344 122 : keys.add("optional","ADAPTIVE_SIGMA_STRIDE","number of steps for measuring adaptive sigma. Default is 10xPACE");
345 122 : keys.add("optional","SIGMA_MIN","never reduce SIGMA below this value");
346 61 : std::string info_biasfactor("the gamma bias factor used for the well-tempered target distribution. ");
347 : if(mode::explore) {
348 : info_biasfactor+="Cannot be 'inf'";
349 : } else {
350 : info_biasfactor+="Set to 'inf' for uniform flat target";
351 : }
352 122 : keys.add("optional","BIASFACTOR",info_biasfactor);
353 122 : keys.add("optional","EPSILON","the value of the regularization constant for the probability");
354 122 : keys.add("optional","KERNEL_CUTOFF","truncate kernels at this distance, in units of sigma");
355 122 : keys.add("optional","NLIST_PARAMETERS","( default=3.0 ) the kernels within sqrt(NLIST_PARAMETERS)*KERNEL_CUTOFF are in the neighbor list. "
356 : "Larger values give a larger list, rebuilt less often. A second value is accepted for backward compatibility, but is not used");
357 122 : keys.addFlag("NLIST",false,"use a neighbor list for the kernels summation. The result is the same as without it, and it is faster when there are many kernels");
358 122 : keys.addFlag("NLIST_PACE_RESET",false,"force the reset of the neighbor list at each PACE. Can be useful with WALKERS_MPI");
359 122 : 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");
360 122 : keys.addFlag("RECURSIVE_MERGE_OFF",false,"do not recursively attempt kernel merging when a new one is added");
361 122 : keys.addFlag("NO_ZED",false,"do not normalize over the explored CV space, Z_n=1");
362 : //kernels and state files
363 122 : keys.add("compulsory","FILE","KERNELS","a file in which the list of all deposited kernels is stored");
364 122 : keys.add("optional","FMT","specify format for KERNELS file");
365 122 : keys.add("optional","STATE_RFILE","read from this file the compressed kernels and all the info needed to RESTART the simulation");
366 122 : keys.add("optional","STATE_WFILE","write to this file the compressed kernels and all the info needed to RESTART the simulation");
367 122 : 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)");
368 122 : keys.addFlag("STORE_STATES",false,"append to STATE_WFILE instead of ovewriting it each time");
369 : //miscellaneous
370 122 : keys.add("optional","EXCLUDED_REGION","kernels are not deposited when the action provided here has a nonzero value, see example above");
371 : if(!mode::explore) {
372 68 : 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)");
373 : }
374 122 : keys.addFlag("CALC_WORK",false,"calculate the total accumulated work done by the bias since last restart");
375 122 : keys.addFlag("WALKERS_MPI",false,"switch on MPI version of multiple walkers");
376 122 : keys.addFlag("SERIAL",false,"perform calculations in serial");
377 61 : keys.use("RESTART");
378 61 : keys.use("UPDATE_FROM");
379 61 : keys.use("UPDATE_UNTIL");
380 :
381 : //output components
382 122 : keys.addOutputComponent("rct","default","estimate of c(t). log(exp(beta V)/beta, should become flat as the simulation converges. Do NOT use for reweighting");
383 122 : keys.addOutputComponent("zed","default","estimate of Z_n. should become flat once no new CV-space region is explored");
384 122 : keys.addOutputComponent("neff","default","effective sample size");
385 122 : keys.addOutputComponent("nker","default","total number of compressed kernels used to represent the bias");
386 122 : keys.addOutputComponent("work","CALC_WORK","total accumulated work done by the bias");
387 122 : keys.addOutputComponent("nlker","NLIST","number of kernels in the neighbor list");
388 122 : keys.addOutputComponent("nlsteps","NLIST","number of steps from last neighbor list update");
389 61 : }
390 :
391 : template <class mode>
392 53 : OPESmetad<mode>::OPESmetad(const ActionOptions& ao)
393 : : PLUMED_BIAS_INIT(ao)
394 53 : , isFirstStep_(true)
395 53 : , counter_(1)
396 53 : , ncv_(getNumberOfArguments())
397 53 : , Zed_(1)
398 53 : , work_(0)
399 53 : , excluded_region_(NULL) {
400 106 : std::string error_in_input1("Error in input in action "+getName()+" with label "+getLabel()+": the keyword ");
401 65 : std::string error_in_input2(" could not be read correctly");
402 :
403 : //set kbt_
404 53 : const double kB=getKBoltzmann();
405 53 : kbt_=getkBT();
406 :
407 : //other compulsory input
408 53 : long long pace=0;
409 53 : parse("PACE",pace);
410 61 : plumed_massert(pace>0 && static_cast<unsigned long long>(pace)<=std::numeric_limits<unsigned>::max(),
411 : "PACE must be greater than zero and no larger than "+std::to_string(std::numeric_limits<unsigned>::max()));
412 49 : stride_=static_cast<unsigned>(pace);
413 :
414 49 : double barrier=0;
415 61 : parse("BARRIER",barrier);
416 49 : plumed_massert(barrier>=0,"the BARRIER should be greater than zero");
417 :
418 49 : biasfactor_=barrier/kbt_;
419 : std::string biasfactor_str;
420 98 : parse("BIASFACTOR",biasfactor_str);
421 49 : if(biasfactor_str=="inf" || biasfactor_str=="INF") {
422 4 : biasfactor_=std::numeric_limits<double>::infinity();
423 4 : bias_prefactor_=1;
424 : } else {
425 45 : if(biasfactor_str.length()>0) {
426 3 : plumed_massert(Tools::convertNoexcept(biasfactor_str,biasfactor_),error_in_input1+"BIASFACTOR"+error_in_input2);
427 : }
428 45 : plumed_massert(biasfactor_>1,"BIASFACTOR must be greater than one (use 'inf' for uniform target)");
429 45 : bias_prefactor_=1-1./biasfactor_;
430 : }
431 : if(mode::explore) {
432 21 : plumed_massert(!std::isinf(biasfactor_),"BIASFACTOR=inf is not compatible with EXPLORE mode");
433 21 : bias_prefactor_=biasfactor_-1;
434 : }
435 :
436 49 : adaptive_sigma_=false;
437 49 : adaptive_sigma_stride_=0;
438 106 : parse("ADAPTIVE_SIGMA_STRIDE",adaptive_sigma_stride_);
439 : std::vector<std::string> sigma_str;
440 49 : parseVector("SIGMA",sigma_str);
441 49 : sigma0_.resize(ncv_);
442 : double dummy;
443 49 : if(sigma_str.size()==1 && !Tools::convertNoexcept(sigma_str[0],dummy)) {
444 21 : plumed_massert(sigma_str[0]=="ADAPTIVE" || sigma_str[0]=="adaptive",error_in_input1+"SIGMA"+error_in_input2);
445 21 : plumed_massert(!std::isinf(biasfactor_),"cannot use BIASFACTOR=inf with adaptive SIGMA");
446 21 : adaptive_counter_=0;
447 21 : if(adaptive_sigma_stride_==0) {
448 10 : adaptive_sigma_stride_=10*stride_; //NB: this is arbitrary, chosen from few tests
449 : }
450 21 : av_cv_.resize(ncv_,0);
451 21 : av_M2_.resize(ncv_,0);
452 21 : plumed_massert(adaptive_sigma_stride_>=stride_,"better to chose ADAPTIVE_SIGMA_STRIDE > PACE");
453 21 : adaptive_sigma_=true;
454 : } else {
455 28 : plumed_massert(sigma_str.size()==ncv_,"number of SIGMA parameters does not match number of arguments");
456 28 : plumed_massert(adaptive_sigma_stride_==0,"if SIGMA is not ADAPTIVE you cannot set an ADAPTIVE_SIGMA_STRIDE");
457 61 : for(unsigned i=0; i<ncv_; i++) {
458 37 : plumed_massert(Tools::convertNoexcept(sigma_str[i],sigma0_[i]),error_in_input1+"SIGMA"+error_in_input2);
459 41 : plumed_massert(std::isfinite(sigma0_[i]) && sigma0_[i]>0,"SIGMA values must be finite and greater than zero");
460 : if(mode::explore) {
461 13 : sigma0_[i]*=std::sqrt(biasfactor_); //the sigma of the target is broader Ftg(s)=1/gamma*F(s)
462 : }
463 : }
464 : }
465 90 : parseVector("SIGMA_MIN",sigma_min_);
466 45 : plumed_massert(sigma_min_.size()==0 || sigma_min_.size()==ncv_,"number of SIGMA_MIN does not match number of arguments");
467 51 : for(double sigma : sigma_min_) {
468 14 : plumed_massert(std::isfinite(sigma) && sigma>0,"SIGMA_MIN values must be finite and greater than zero");
469 : }
470 41 : if(sigma_min_.size()>0 && !adaptive_sigma_) {
471 3 : for(unsigned i=0; i<ncv_; i++) {
472 2 : plumed_massert(sigma_min_[i]<=sigma0_[i],"SIGMA_MIN should be smaller than SIGMA");
473 : }
474 : }
475 :
476 41 : epsilon_=std::exp(-barrier/bias_prefactor_/kbt_);
477 41 : parse("EPSILON",epsilon_);
478 41 : plumed_massert(epsilon_>0,"you must choose a value for EPSILON greater than zero. Is your BARRIER too high?");
479 41 : sum_weights_=std::pow(epsilon_,bias_prefactor_); //to avoid NANs we start with counter_=1 and w0=exp(beta*V0)
480 41 : sum_weights2_=sum_weights_*sum_weights_;
481 :
482 41 : double cutoff=sqrt(2.*barrier/bias_prefactor_/kbt_);
483 : if(mode::explore) {
484 17 : cutoff=sqrt(2.*barrier/kbt_); //otherwise it is too small
485 : }
486 41 : parse("KERNEL_CUTOFF",cutoff);
487 41 : plumed_massert(cutoff>0,"you must choose a value for KERNEL_CUTOFF greater than zero");
488 41 : cutoff2_=cutoff*cutoff;
489 41 : val_at_cutoff_=std::exp(-0.5*cutoff2_);
490 :
491 41 : threshold2_=1;
492 41 : parse("COMPRESSION_THRESHOLD",threshold2_);
493 41 : threshold2_*=threshold2_;
494 41 : if(threshold2_!=0) {
495 41 : plumed_massert(threshold2_>0 && threshold2_<cutoff2_,"COMPRESSION_THRESHOLD cannot be bigger than the KERNEL_CUTOFF");
496 : }
497 :
498 : //setup neighbor list
499 41 : nlist_=false;
500 41 : parseFlag("NLIST",nlist_);
501 41 : nlist_pace_reset_=false;
502 41 : parseFlag("NLIST_PACE_RESET",nlist_pace_reset_);
503 41 : if(nlist_pace_reset_) {
504 2 : nlist_=true;
505 : }
506 : std::vector<double> nlist_param;
507 41 : parseVector("NLIST_PARAMETERS",nlist_param);
508 41 : nlist_param_=3.0;//*cutoff2_ -> max distance of neighbors
509 41 : if(nlist_param.size()>0) {
510 6 : nlist_=true;
511 6 : plumed_massert(nlist_param.size()<=2,"NLIST_PARAMETERS takes one value");
512 6 : plumed_massert(nlist_param[0]>1.0,"NLIST_PARAMETERS must be greater than 1");
513 6 : nlist_param_=nlist_param[0];
514 : }
515 41 : nlist_center_.resize(ncv_);
516 41 : nlist_steps_=0;
517 41 : nlist_update_=true;
518 : //the list is rebuilt before any kernel outside it can get within the cutoff,
519 : //i.e. when the CV moves out of a skin of (sqrt(nlist_param_)-1)*cutoff, in units of the smallest sigma
520 41 : nlist_skin2_=std::pow(std::sqrt(nlist_param_)-1,2)*cutoff2_;
521 41 : nlist_sigma_min_.resize(ncv_,std::numeric_limits<double>::max());
522 41 : nlist_all_outdated_=true;
523 :
524 : //setup the sum of the kernels
525 41 : sum_outdated_=true;
526 41 : period_.resize(ncv_,0.);
527 41 : inv_period_.resize(ncv_,0.);
528 41 : cv_.resize(ncv_);
529 41 : der_prob_.resize(ncv_);
530 111 : for(unsigned i=0; i<ncv_; i++) {
531 70 : if(getPntrToArgument(i)->isPeriodic()) {
532 44 : period_[i]=getPntrToArgument(i)->getMaxMinusMin();
533 44 : inv_period_[i]=1./period_[i];
534 : }
535 : }
536 :
537 : //optional stuff
538 41 : no_Zed_=false;
539 41 : parseFlag("NO_ZED",no_Zed_);
540 41 : if(no_Zed_) {
541 : //this makes it more gentle in the initial phase
542 6 : sum_weights_=1;
543 6 : sum_weights2_=1;
544 : }
545 41 : fixed_sigma_=false;
546 41 : parseFlag("FIXED_SIGMA",fixed_sigma_);
547 41 : bool recursive_merge_off=false;
548 41 : parseFlag("RECURSIVE_MERGE_OFF",recursive_merge_off);
549 41 : recursive_merge_=!recursive_merge_off;
550 82 : parseFlag("CALC_WORK",calc_work_);
551 :
552 : //options involving extra arguments
553 : std::vector<Value*> args;
554 82 : parseArgumentList("EXCLUDED_REGION",args);
555 41 : if(args.size()>0) {
556 2 : plumed_massert(args.size()==1,"only one characteristic function should define the region to be excluded");
557 2 : requestExtraDependencies(args);
558 2 : excluded_region_=args[0];
559 : }
560 : if(!mode::explore) {
561 48 : parseArgumentList("EXTRA_BIAS",extra_biases_);
562 24 : if(extra_biases_.size()>0) {
563 2 : requestExtraDependencies(extra_biases_);
564 : }
565 : }
566 :
567 : //kernels file
568 : std::string kernelsFileName;
569 82 : parse("FILE",kernelsFileName);
570 : std::string fmt;
571 82 : parse("FMT",fmt);
572 :
573 : //output checkpoint of current state
574 : std::string restartFileName;
575 82 : parse("STATE_RFILE",restartFileName);
576 : std::string stateFileName;
577 41 : parse("STATE_WFILE",stateFileName);
578 41 : wStateStride_=0;
579 41 : parse("STATE_WSTRIDE",wStateStride_);
580 41 : storeOldStates_=false;
581 41 : parseFlag("STORE_STATES",storeOldStates_);
582 41 : if(wStateStride_!=0 || storeOldStates_) {
583 12 : plumed_massert(stateFileName.length()>0,"filename for storing simulation status not specified, use STATE_WFILE");
584 : }
585 41 : if(wStateStride_>0) {
586 12 : plumed_massert(wStateStride_>=(int)stride_,"STATE_WSTRIDE is in units of MD steps, thus it is suggested to use a multiple of PACE");
587 : }
588 41 : if(stateFileName.length()>0 && wStateStride_==0) {
589 1 : wStateStride_=-1; //will print only on CPT events (checkpoints set by some MD engines, like gromacs)
590 : }
591 :
592 : //multiple walkers //TODO implement also external mw for cp2k
593 41 : bool walkers_mpi=false;
594 41 : parseFlag("WALKERS_MPI",walkers_mpi);
595 41 : if(walkers_mpi) {
596 : //If this Action is not compiled with MPI the user is informed and we exit gracefully
597 10 : plumed_massert(Communicator::plumedHasMPI(),"Invalid walkers configuration: WALKERS_MPI flag requires MPI compilation");
598 10 : plumed_massert(Communicator::initialized(),"Invalid walkers configuration: WALKERS_MPI needs the communicator correctly initialized.");
599 :
600 10 : if(comm.Get_rank()==0) { //multi_sim_comm works on first rank only
601 10 : NumWalkers_=multi_sim_comm.Get_size();
602 10 : walker_rank_=multi_sim_comm.Get_rank();
603 : }
604 10 : comm.Bcast(NumWalkers_,0); //if each walker has more than one processor update them all
605 10 : comm.Bcast(walker_rank_,0);
606 : } else {
607 31 : NumWalkers_=1;
608 31 : walker_rank_=0;
609 : }
610 :
611 : //parallelization stuff
612 41 : NumOMP_=OpenMP::getNumThreads();
613 41 : NumParallel_=comm.Get_size();
614 41 : rank_=comm.Get_rank();
615 41 : bool serial=false;
616 41 : parseFlag("SERIAL",serial);
617 41 : if(serial) {
618 4 : NumOMP_=1;
619 4 : NumParallel_=1;
620 4 : rank_=0;
621 : }
622 :
623 41 : checkRead();
624 :
625 : //restart if needed
626 : bool convertKernelsToState=false;
627 41 : if(getRestart()) {
628 : bool stateRestart=true;
629 13 : if(restartFileName.length()==0) {
630 : stateRestart=false;
631 : restartFileName=kernelsFileName;
632 : }
633 13 : IFile ifile;
634 13 : ifile.link(*this);
635 13 : if(ifile.FileExist(restartFileName)) {
636 13 : bool tmp_nlist=nlist_;
637 13 : nlist_=false; // NLIST is not needed while restarting
638 13 : ifile.open(restartFileName);
639 13 : log.printf(" RESTART - make sure all used options are compatible\n");
640 13 : log.printf(" restarting from: %s\n",restartFileName.c_str());
641 13 : std::string action_name=getName();
642 13 : if(stateRestart) {
643 6 : log.printf(" it should be a STATE file (not a KERNELS file)\n");
644 : action_name+="_state";
645 : } else {
646 7 : log.printf(" +++ WARNING +++ RESTART from KERNELS might be approximate, use STATE_WFILE and STATE_RFILE to restart from the exact state\n");
647 : action_name+="_kernels";
648 : }
649 : std::string old_action_name;
650 13 : ifile.scanField("action",old_action_name);
651 13 : plumed_massert(action_name==old_action_name,"RESTART - mismatch between old and new action name. Expected '"+action_name+"', but found '"+old_action_name+"'");
652 : std::string old_biasfactor_str;
653 26 : ifile.scanField("biasfactor",old_biasfactor_str);
654 13 : if(old_biasfactor_str=="inf" || old_biasfactor_str=="INF") {
655 1 : if(!std::isinf(biasfactor_)) {
656 0 : log.printf(" +++ WARNING +++ previous bias factor was inf while now it is %g\n",biasfactor_);
657 : }
658 : } else {
659 : double old_biasfactor;
660 12 : ifile.scanField("biasfactor",old_biasfactor);
661 12 : if(std::abs(biasfactor_-old_biasfactor)>1e-6*biasfactor_) {
662 0 : log.printf(" +++ WARNING +++ previous bias factor was %g while now it is %g. diff = %g\n",old_biasfactor,biasfactor_,biasfactor_-old_biasfactor);
663 : }
664 : }
665 : double old_epsilon;
666 13 : ifile.scanField("epsilon",old_epsilon);
667 13 : if(std::abs(epsilon_-old_epsilon)>1e-6*epsilon_) {
668 8 : log.printf(" +++ WARNING +++ previous epsilon was %g while now it is %g. diff = %g\n",old_epsilon,epsilon_,epsilon_-old_epsilon);
669 : }
670 : double old_cutoff;
671 13 : ifile.scanField("kernel_cutoff",old_cutoff);
672 13 : if(std::abs(cutoff-old_cutoff)>1e-6*cutoff) {
673 0 : log.printf(" +++ WARNING +++ previous kernel_cutoff was %g while now it is %g. diff = %g\n",old_cutoff,cutoff,cutoff-old_cutoff);
674 : }
675 : double old_threshold;
676 13 : const double threshold=sqrt(threshold2_);
677 13 : ifile.scanField("compression_threshold",old_threshold);
678 13 : if(std::abs(threshold-old_threshold)>1e-6*threshold) {
679 0 : log.printf(" +++ WARNING +++ previous compression_threshold was %g while now it is %g. diff = %g\n",old_threshold,threshold,threshold-old_threshold);
680 : }
681 13 : if(stateRestart) {
682 6 : ifile.scanField("zed",Zed_);
683 6 : ifile.scanField("sum_weights",sum_weights_);
684 6 : ifile.scanField("sum_weights2",sum_weights2_);
685 6 : ifile.scanField("counter",counter_);
686 6 : if(adaptive_sigma_) {
687 6 : ifile.scanField("adaptive_counter",adaptive_counter_);
688 6 : if(NumWalkers_==1) {
689 6 : for(unsigned i=0; i<ncv_; i++) {
690 8 : ifile.scanField("sigma0_"+getPntrToArgument(i)->getName(),sigma0_[i]);
691 8 : ifile.scanField("av_cv_"+getPntrToArgument(i)->getName(),av_cv_[i]);
692 8 : ifile.scanField("av_M2_"+getPntrToArgument(i)->getName(),av_M2_[i]);
693 : }
694 : } else {
695 12 : for(unsigned w=0; w<NumWalkers_; w++)
696 24 : for(unsigned i=0; i<ncv_; i++) {
697 : double tmp0,tmp1,tmp2;
698 32 : const std::string arg_iw=getPntrToArgument(i)->getName()+"_"+std::to_string(w);
699 16 : ifile.scanField("sigma0_"+arg_iw,tmp0);
700 16 : ifile.scanField("av_cv_"+arg_iw,tmp1);
701 16 : ifile.scanField("av_M2_"+arg_iw,tmp2);
702 16 : if(w==walker_rank_) {
703 8 : sigma0_[i]=tmp0;
704 8 : av_cv_[i]=tmp1;
705 8 : av_M2_[i]=tmp2;
706 : }
707 : }
708 : }
709 : }
710 : }
711 39 : for(unsigned i=0; i<ncv_; i++) {
712 26 : if(getPntrToArgument(i)->isPeriodic()) {
713 : std::string arg_min,arg_max;
714 26 : getPntrToArgument(i)->getDomain(arg_min,arg_max);
715 : std::string file_min,file_max;
716 52 : ifile.scanField("min_"+getPntrToArgument(i)->getName(),file_min);
717 26 : ifile.scanField("max_"+getPntrToArgument(i)->getName(),file_max);
718 26 : plumed_massert(file_min==arg_min,"RESTART - mismatch between old and new ARG periodicity");
719 26 : plumed_massert(file_max==arg_max,"RESTART - mismatch between old and new ARG periodicity");
720 : }
721 : }
722 13 : if(stateRestart) {
723 : double time;
724 60 : while(ifile.scanField("time",time)) {
725 24 : std::vector<double> center(ncv_);
726 24 : std::vector<double> sigma(ncv_);
727 : double height;
728 72 : for(unsigned i=0; i<ncv_; i++) {
729 48 : ifile.scanField(getPntrToArgument(i)->getName(),center[i]);
730 : }
731 72 : for(unsigned i=0; i<ncv_; i++) {
732 96 : ifile.scanField("sigma_"+getPntrToArgument(i)->getName(),sigma[i]);
733 : }
734 24 : ifile.scanField("height",height);
735 24 : ifile.scanField();
736 24 : kernels_.emplace_back(height,center,sigma);
737 : }
738 6 : log.printf(" a total of %lu kernels where read\n",kernels_.size());
739 : } else {
740 7 : ifile.allowIgnoredFields(); //this allows for multiple restart, but without checking for consistency between them!
741 : double time;
742 358 : while(ifile.scanField("time",time)) {
743 172 : std::vector<double> center(ncv_);
744 172 : std::vector<double> sigma(ncv_);
745 : double height;
746 : double logweight;
747 516 : for(unsigned i=0; i<ncv_; i++) {
748 344 : ifile.scanField(getPntrToArgument(i)->getName(),center[i]);
749 : }
750 516 : for(unsigned i=0; i<ncv_; i++) {
751 688 : ifile.scanField("sigma_"+getPntrToArgument(i)->getName(),sigma[i]);
752 : }
753 172 : ifile.scanField("height",height);
754 172 : ifile.scanField("logweight",logweight);
755 172 : ifile.scanField();
756 172 : if(counter_==(1+walker_rank_) && adaptive_sigma_) {
757 : //the first kernel has sigma=s*sigma0, where s is the rescaling at the first deposition,
758 : //and height=w*prod_i(sigma0_i/sigma_i)=w/s^ncv, with w=1 for EXPLORE, so s can be recovered
759 : //(exactly, unless SIGMA_MIN was used)
760 2 : const double prod_ratio=mode::explore?height:height*std::exp(-logweight);
761 2 : const double s=std::pow(prod_ratio,-1./ncv_);
762 6 : for(unsigned i=0; i<ncv_; i++) {
763 4 : sigma0_[i]=sigma[i]/s;
764 : }
765 : }
766 172 : addKernel(height,center,sigma);
767 172 : const double weight=std::exp(logweight);
768 172 : sum_weights_+=weight; //this sum is slightly inaccurate, because when printing some precision is lost
769 172 : sum_weights2_+=weight*weight;
770 172 : counter_++;
771 : }
772 7 : KDEnorm_=mode::explore?counter_:sum_weights_;
773 7 : if(!no_Zed_) {
774 : //all kernels evaluated at all centers, with the kernel sum used for the bias (nlist_ is off here).
775 : //Each MPI process takes a contiguous share of the centers, split among the OpenMP threads
776 4 : updateSum();
777 4 : const unsigned size=kernels_.size();
778 4 : const unsigned share=(size+NumParallel_-1)/NumParallel_;
779 4 : const unsigned begin=std::min(rank_*share,size);
780 4 : const unsigned end=std::min(begin+share,size);
781 4 : double sum_uprob=0;
782 4 : #pragma omp parallel num_threads(NumOMP_)
783 : {
784 : std::vector<double> unused_der(ncv_); //derivatives are not needed
785 : #pragma omp for reduction(+:sum_uprob) nowait
786 : for(unsigned k=begin; k<end; k++) {
787 : sumKernels(kernels_[k].center,0,size,sum_uprob,unused_der);
788 : }
789 : }
790 4 : if(NumParallel_>1) {
791 0 : comm.Sum(sum_uprob);
792 : }
793 4 : Zed_=sum_uprob/KDEnorm_/kernels_.size();
794 4 : sum_outdated_=true; //the neighbor list, if any, is not set yet
795 : }
796 7 : log.printf(" a total of %llu kernels where read, and compressed to %lu\n",counter_-1,kernels_.size());
797 : convertKernelsToState=true;
798 : }
799 13 : ifile.reset(false);
800 13 : ifile.close();
801 13 : nlist_=tmp_nlist;
802 : } else { //same behaviour as METAD
803 0 : plumed_merror("RESTART requested, but file '"+restartFileName+"' was not found!\n Set RESTART=NO or provide a restart file");
804 : }
805 13 : if(NumWalkers_>1) { //make sure that all walkers are doing the same thing
806 6 : const unsigned kernels_size=kernels_.size();
807 6 : std::vector<unsigned> all_kernels_size(NumWalkers_);
808 6 : if(comm.Get_rank()==0) {
809 6 : multi_sim_comm.Allgather(kernels_size,all_kernels_size);
810 : }
811 6 : comm.Bcast(all_kernels_size,0);
812 : bool same_number_of_kernels=true;
813 12 : for(unsigned w=1; w<NumWalkers_; w++)
814 6 : if(all_kernels_size[0]!=all_kernels_size[w]) {
815 : same_number_of_kernels=false;
816 : }
817 6 : plumed_massert(same_number_of_kernels,"RESTART - not all walkers are reading the same file!");
818 : }
819 41 : } else if(restartFileName.length()>0) {
820 4 : log.printf(" +++ WARNING +++ the provided STATE_RFILE will be ignored, since RESTART was not requested\n");
821 : }
822 :
823 : //sync all walkers to avoid opening files before reading is over (see also METAD)
824 41 : comm.Barrier();
825 41 : if(comm.Get_rank()==0 && walkers_mpi) {
826 10 : multi_sim_comm.Barrier();
827 : }
828 :
829 : //setup output kernels file
830 41 : kernelsOfile_.link(*this);
831 41 : if(NumWalkers_>1) {
832 10 : if(walker_rank_>0) {
833 : kernelsFileName="/dev/null"; //only first walker writes on file
834 : }
835 20 : kernelsOfile_.enforceSuffix("");
836 : }
837 41 : kernelsOfile_.open(kernelsFileName);
838 41 : if(fmt.length()>0) {
839 46 : kernelsOfile_.fmtField(" "+fmt);
840 : }
841 : kernelsOfile_.setHeavyFlush(); //do I need it?
842 : //define and set const fields
843 41 : kernelsOfile_.addConstantField("action");
844 41 : kernelsOfile_.addConstantField("biasfactor");
845 41 : kernelsOfile_.addConstantField("epsilon");
846 41 : kernelsOfile_.addConstantField("kernel_cutoff");
847 41 : kernelsOfile_.addConstantField("compression_threshold");
848 111 : for(unsigned i=0; i<ncv_; i++) {
849 70 : kernelsOfile_.setupPrintValue(getPntrToArgument(i));
850 : }
851 82 : kernelsOfile_.printField("action",getName()+"_kernels");
852 41 : kernelsOfile_.printField("biasfactor",biasfactor_);
853 41 : kernelsOfile_.printField("epsilon",epsilon_);
854 41 : kernelsOfile_.printField("kernel_cutoff",sqrt(cutoff2_));
855 41 : kernelsOfile_.printField("compression_threshold",sqrt(threshold2_));
856 :
857 : //open file for storing state
858 41 : if(wStateStride_!=0) {
859 13 : stateOfile_.link(*this);
860 13 : if(NumWalkers_>1) {
861 8 : if(walker_rank_>0) {
862 : stateFileName="/dev/null"; //only first walker writes on file
863 : }
864 16 : stateOfile_.enforceSuffix("");
865 : }
866 13 : stateOfile_.open(stateFileName);
867 13 : if(fmt.length()>0) {
868 26 : stateOfile_.fmtField(" "+fmt);
869 : }
870 13 : if(convertKernelsToState) {
871 2 : dumpStateToFile();
872 : }
873 : }
874 :
875 : //set initial old values
876 41 : KDEnorm_=mode::explore?counter_:sum_weights_;
877 41 : old_KDEnorm_=KDEnorm_;
878 :
879 : //add and set output components
880 82 : addComponent("rct");
881 41 : componentIsNotPeriodic("rct");
882 82 : getPntrToComponent("rct")->set(kbt_*std::log(sum_weights_/counter_));
883 82 : addComponent("zed");
884 41 : componentIsNotPeriodic("zed");
885 82 : getPntrToComponent("zed")->set(Zed_);
886 82 : addComponent("neff");
887 41 : componentIsNotPeriodic("neff");
888 82 : getPntrToComponent("neff")->set(std::pow(1+sum_weights_,2)/(1+sum_weights2_));
889 82 : addComponent("nker");
890 41 : componentIsNotPeriodic("nker");
891 41 : getPntrToComponent("nker")->set(kernels_.size());
892 41 : if(calc_work_) {
893 14 : addComponent("work");
894 14 : componentIsNotPeriodic("work");
895 : }
896 41 : if(nlist_) {
897 18 : addComponent("nlker");
898 18 : componentIsNotPeriodic("nlker");
899 18 : addComponent("nlsteps");
900 18 : componentIsNotPeriodic("nlsteps");
901 : }
902 :
903 : //printing some info
904 41 : log.printf(" temperature = %g\n",kbt_/kB);
905 41 : log.printf(" beta = %g\n",1./kbt_);
906 41 : log.printf(" depositing new kernels with PACE = %u\n",stride_);
907 41 : log.printf(" expected BARRIER is %g\n",barrier);
908 41 : log.printf(" using target distribution with BIASFACTOR gamma = %g\n",biasfactor_);
909 41 : if(std::isinf(biasfactor_)) {
910 4 : log.printf(" (thus a uniform flat target distribution, no well-tempering)\n");
911 : }
912 41 : if(excluded_region_!=NULL) {
913 2 : log.printf(" -- EXCLUDED_REGION: kernels will be deposited only when '%s' is equal to zero\n",excluded_region_->getName().c_str());
914 : }
915 41 : if(extra_biases_.size()>0) {
916 2 : log.printf(" -- EXTRA_BIAS: ");
917 5 : for(unsigned e=0; e<extra_biases_.size(); e++) {
918 3 : log.printf("%s ",extra_biases_[e]->getName().c_str());
919 : }
920 2 : log.printf("will be reweighted\n");
921 : }
922 41 : if(adaptive_sigma_) {
923 21 : log.printf(" adaptive SIGMA will be used, with ADAPTIVE_SIGMA_STRIDE = %u\n",adaptive_sigma_stride_);
924 21 : unsigned x=std::ceil(adaptive_sigma_stride_/stride_);
925 21 : log.printf(" thus the first x kernel depositions will be skipped, x = ADAPTIVE_SIGMA_STRIDE/PACE = %u\n",x);
926 : } else {
927 20 : log.printf(" kernels have initial SIGMA = ");
928 49 : for(unsigned i=0; i<ncv_; i++) {
929 29 : log.printf(" %g",sigma0_[i]);
930 : }
931 20 : log.printf("\n");
932 : }
933 41 : if(sigma_min_.size()>0) {
934 3 : log.printf(" kernels have a SIGMA_MIN = ");
935 9 : for(unsigned i=0; i<ncv_; i++) {
936 6 : log.printf(" %g",sigma_min_[i]);
937 : }
938 3 : log.printf("\n");
939 : }
940 41 : if(fixed_sigma_) {
941 14 : log.printf(" -- FIXED_SIGMA: sigma will not decrease as the simulation proceeds\n");
942 : }
943 41 : log.printf(" kernels are truncated with KERNELS_CUTOFF = %g\n",cutoff);
944 41 : if(cutoff<3.5) {
945 10 : log.printf(" +++ WARNING +++ probably kernels are truncated too much\n");
946 : }
947 41 : log.printf(" the value at cutoff is = %g\n",val_at_cutoff_);
948 41 : log.printf(" regularization EPSILON = %g\n",epsilon_);
949 41 : if(val_at_cutoff_>epsilon_*(1+1e-6)) {
950 7 : log.printf(" +++ WARNING +++ the KERNEL_CUTOFF might be too small for the given EPSILON\n");
951 : }
952 41 : log.printf(" kernels will be compressed when closer than COMPRESSION_THRESHOLD = %g\n",sqrt(threshold2_));
953 41 : if(threshold2_==0) {
954 0 : log.printf(" +++ WARNING +++ kernels will never merge, expect slowdowns\n");
955 : }
956 41 : if(!recursive_merge_) {
957 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");
958 : }
959 41 : if(nlist_) {
960 9 : log.printf(" -- NLIST: using neighbor list for kernels, with parameter: %g\n",nlist_param_);
961 9 : if(nlist_param.size()==2) {
962 4 : log.printf(" +++ WARNING +++ the second of NLIST_PARAMETERS is not used anymore, the neighbor list is now exact\n");
963 : }
964 : }
965 41 : if(nlist_pace_reset_) {
966 2 : log.printf(" -- NLIST_PACE_RESET: forcing the neighbor list to update every PACE\n");
967 : }
968 41 : if(no_Zed_) {
969 6 : log.printf(" -- NO_ZED: using fixed normalization factor = %g\n",Zed_);
970 : }
971 41 : if(wStateStride_>0) {
972 12 : log.printf(" state checkpoints are written on file %s every %d MD steps\n",stateFileName.c_str(),wStateStride_);
973 : }
974 41 : if(wStateStride_==-1) {
975 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());
976 : }
977 41 : if(walkers_mpi) {
978 10 : log.printf(" -- WALKERS_MPI: if multiple replicas are present, they will share the same bias via MPI\n");
979 : }
980 41 : if(NumWalkers_>1) {
981 10 : log.printf(" using multiple walkers\n");
982 10 : log.printf(" number of walkers: %u\n",NumWalkers_);
983 10 : log.printf(" walker rank: %u\n",walker_rank_);
984 : }
985 41 : int mw_warning=0;
986 41 : if(!walkers_mpi && comm.Get_rank()==0 && multi_sim_comm.Get_size()>(int)NumWalkers_) {
987 0 : mw_warning=1;
988 : }
989 41 : comm.Bcast(mw_warning,0);
990 41 : if(mw_warning) { //log.printf messes up with comm, so never use it without Bcast!
991 0 : log.printf(" +++ WARNING +++ multiple replicas will NOT communicate unless the flag WALKERS_MPI is used\n");
992 : }
993 41 : if(NumParallel_>1) {
994 4 : log.printf(" using multiple MPI processes per simulation: %u\n",NumParallel_);
995 : }
996 41 : if(NumOMP_>1) {
997 37 : log.printf(" using multiple OpenMP threads per simulation: %u\n",NumOMP_);
998 : }
999 41 : if(serial) {
1000 4 : log.printf(" -- SERIAL: no loop parallelization, despite %d MPI processes and %u OpenMP threads available\n",comm.Get_size(),OpenMP::getNumThreads());
1001 : }
1002 41 : log.printf(" Bibliography: ");
1003 82 : log<<plumed.cite("M. Invernizzi and M. Parrinello, J. Phys. Chem. Lett. 11, 2731-2736 (2020)");
1004 24 : if(mode::explore || adaptive_sigma_) {
1005 58 : log<<plumed.cite("M. Invernizzi and M. Parrinello, J. Chem. Theory Comput. 18, 3988-3996 (2022)");
1006 : }
1007 41 : log.printf("\n");
1008 162 : }
1009 :
1010 : template <class mode>
1011 65247 : void OPESmetad<mode>::calculate() {
1012 : //get cv
1013 65247 : std::vector<double>& cv=cv_;
1014 195565 : for(unsigned i=0; i<ncv_; i++) {
1015 130318 : cv[i]=getArgument(i);
1016 : }
1017 :
1018 : //check neighbor list
1019 65247 : if(nlist_) {
1020 32255 : nlist_steps_++;
1021 32255 : if(getExchangeStep()) {
1022 0 : nlist_update_=true;
1023 : } else {
1024 : //kernels outside the list cannot be within the cutoff as long as this displacement is within the skin
1025 : double norm2=0;
1026 96765 : for(unsigned i=0; i<ncv_; i++) {
1027 64510 : double diff=cv[i]-nlist_center_[i];
1028 64510 : diff-=period_[i]*std::nearbyint(diff*inv_period_[i]); //zero for non-periodic CVs
1029 64510 : const double dist_i=diff/nlist_sigma_min_[i];
1030 64510 : norm2+=dist_i*dist_i;
1031 : }
1032 32255 : if(norm2>nlist_skin2_) {
1033 16446 : nlist_update_=true;
1034 : }
1035 : }
1036 32255 : if(nlist_update_) {
1037 16665 : updateNlist(cv);
1038 : }
1039 : }
1040 :
1041 : //set bias and forces
1042 65247 : std::vector<double>& der_prob=der_prob_;
1043 : std::fill(der_prob.begin(),der_prob.end(),0.);
1044 65247 : const double prob=getProbAndDerivatives(cv,der_prob);
1045 65247 : const double bias=kbt_*bias_prefactor_*std::log(prob/Zed_+epsilon_);
1046 65247 : setBias(bias);
1047 195565 : for(unsigned i=0; i<ncv_; i++) {
1048 130318 : setOutputForce(i,-kbt_*bias_prefactor_/(prob/Zed_+epsilon_)*der_prob[i]/Zed_);
1049 : }
1050 65247 : }
1051 :
1052 : template <class mode>
1053 65213 : void OPESmetad<mode>::update() {
1054 65213 : if(isFirstStep_) { //same in MetaD, useful for restarts?
1055 41 : isFirstStep_=false;
1056 41 : return;
1057 : }
1058 :
1059 : //update variance if adaptive sigma
1060 65172 : if(adaptive_sigma_) {
1061 64642 : adaptive_counter_++;
1062 64642 : unsigned tau=adaptive_sigma_stride_;
1063 64642 : if(adaptive_counter_<adaptive_sigma_stride_) {
1064 455 : tau=adaptive_counter_;
1065 : }
1066 193876 : for(unsigned i=0; i<ncv_; i++) {
1067 : //Welford's online algorithm for standard deviation
1068 129234 : const double cv_i=getArgument(i);
1069 129234 : const double diff_i=difference(i,av_cv_[i],cv_i);
1070 129234 : av_cv_[i]+=diff_i/tau; //exponentially decaying average
1071 129234 : av_M2_[i]+=diff_i*difference(i,av_cv_[i],cv_i);
1072 : }
1073 64642 : if(adaptive_counter_<adaptive_sigma_stride_ && counter_==1) { //counter_>1 if restarting
1074 : return; //do not apply bias before having measured sigma
1075 : }
1076 : }
1077 :
1078 : //do update
1079 64735 : if(getStep()%stride_==0 && (excluded_region_==NULL || excluded_region_->get()==0)) {
1080 13057 : old_KDEnorm_=KDEnorm_;
1081 13057 : delta_kernels_.clear();
1082 13057 : unsigned old_nker=kernels_.size();
1083 :
1084 : //get new kernel height
1085 13057 : double log_weight=getOutputQuantity(0)/kbt_; //first value is always the current bias
1086 13077 : for(unsigned e=0; e<extra_biases_.size(); e++) {
1087 20 : log_weight+=extra_biases_[e]->get()/kbt_; //extra biases contribute to the weight
1088 : }
1089 13057 : double height=std::exp(log_weight);
1090 :
1091 : //update sum_weights_ and neff
1092 13057 : double sum_heights=height;
1093 13057 : double sum_heights2=height*height;
1094 13057 : if(NumWalkers_>1) {
1095 126 : if(comm.Get_rank()==0) {
1096 126 : multi_sim_comm.Sum(sum_heights);
1097 126 : multi_sim_comm.Sum(sum_heights2);
1098 : }
1099 126 : comm.Bcast(sum_heights,0);
1100 126 : comm.Bcast(sum_heights2,0);
1101 : }
1102 13057 : counter_+=NumWalkers_;
1103 13057 : sum_weights_+=sum_heights;
1104 13057 : sum_weights2_+=sum_heights2;
1105 13057 : const double neff=std::pow(1+sum_weights_,2)/(1+sum_weights2_); //adding 1 makes it more robust at the start
1106 13057 : getPntrToComponent("rct")->set(kbt_*std::log(sum_weights_/counter_));
1107 13057 : getPntrToComponent("neff")->set(neff);
1108 : if(mode::explore) {
1109 6468 : KDEnorm_=counter_;
1110 6468 : height=1; //plain KDE, bias reweight does not enter here
1111 : } else {
1112 6589 : KDEnorm_=sum_weights_;
1113 : }
1114 :
1115 : //if needed, rescale sigma and height
1116 13057 : std::vector<double> sigma=sigma0_;
1117 13057 : if(adaptive_sigma_) {
1118 6478 : const double factor=mode::explore?1:biasfactor_;
1119 12901 : if(counter_==1+NumWalkers_) { //first time only
1120 38 : for(unsigned i=0; i<ncv_; i++) {
1121 25 : av_M2_[i]*=biasfactor_; //from unbiased, thus must be adjusted
1122 : }
1123 38 : for(unsigned i=0; i<ncv_; i++) {
1124 25 : sigma0_[i]=std::sqrt(av_M2_[i]/adaptive_counter_/factor);
1125 : }
1126 13 : if(sigma_min_.size()==0) {
1127 38 : for(unsigned i=0; i<ncv_; i++) {
1128 25 : 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");
1129 : }
1130 : } else {
1131 0 : for(unsigned i=0; i<ncv_; i++) {
1132 0 : sigma0_[i]=std::max(sigma0_[i],sigma_min_[i]);
1133 : }
1134 : }
1135 : }
1136 38698 : for(unsigned i=0; i<ncv_; i++) {
1137 25797 : sigma[i]=std::sqrt(av_M2_[i]/adaptive_counter_/factor);
1138 : }
1139 12901 : if(sigma_min_.size()==0) {
1140 38548 : for(unsigned i=0; i<ncv_; i++) {
1141 25697 : if(sigma[i]<1e-6) {
1142 0 : log.printf("+++ WARNING +++ the ADAPTIVE SIGMA is suspiciously small, you should set a safe SIGMA_MIN. 1e-6 will be used here\n");
1143 0 : sigma[i]=1e-6;
1144 0 : sigma_min_.resize(ncv_,1e-6);
1145 : }
1146 : }
1147 : } else {
1148 150 : for(unsigned i=0; i<ncv_; i++) {
1149 100 : sigma[i]=std::max(sigma[i],sigma_min_[i]);
1150 : }
1151 : }
1152 : }
1153 13057 : if(!fixed_sigma_) {
1154 6426 : const double size=mode::explore?counter_:neff; //for EXPLORE neff is not relevant
1155 12973 : const double s_rescaling=std::pow(size*(ncv_+2.)/4.,-1./(4+ncv_));
1156 38898 : for(unsigned i=0; i<ncv_; i++) {
1157 25925 : sigma[i]*=s_rescaling;
1158 : }
1159 12973 : if(sigma_min_.size()>0) {
1160 150 : for(unsigned i=0; i<ncv_; i++) {
1161 100 : sigma[i]=std::max(sigma[i],sigma_min_[i]);
1162 : }
1163 : }
1164 : }
1165 : //the height should be divided by sqrt(2*pi)*sigma0_,
1166 : //but this overall factor would be canceled when dividing by Zed
1167 : //thus we skip it altogether, but keep any other sigma rescaling
1168 39126 : for(unsigned i=0; i<ncv_; i++) {
1169 26069 : height*=(sigma0_[i]/sigma[i]);
1170 : }
1171 :
1172 : //get new kernel center
1173 13057 : std::vector<double> center(ncv_);
1174 39126 : for(unsigned i=0; i<ncv_; i++) {
1175 26069 : center[i]=getArgument(i);
1176 : }
1177 :
1178 : //add new kernel(s)
1179 13057 : if(NumWalkers_==1) {
1180 12931 : addKernel(height,center,sigma,log_weight);
1181 : } else {
1182 126 : std::vector<double> all_height(NumWalkers_,0.0);
1183 126 : std::vector<double> all_center(NumWalkers_*ncv_,0.0);
1184 126 : std::vector<double> all_sigma(NumWalkers_*ncv_,0.0);
1185 126 : std::vector<double> all_logweight(NumWalkers_,0.0);
1186 126 : if(comm.Get_rank()==0) {
1187 126 : multi_sim_comm.Allgather(height,all_height);
1188 126 : multi_sim_comm.Allgather(center,all_center);
1189 126 : multi_sim_comm.Allgather(sigma,all_sigma);
1190 126 : multi_sim_comm.Allgather(log_weight,all_logweight);
1191 : }
1192 126 : comm.Bcast(all_height,0);
1193 126 : comm.Bcast(all_center,0);
1194 126 : comm.Bcast(all_sigma,0);
1195 126 : comm.Bcast(all_logweight,0);
1196 126 : if(nlist_) {
1197 : //gather all the nlist_index_, so merging can be done using it
1198 50 : std::vector<int> all_nlist_size(NumWalkers_);
1199 50 : if(comm.Get_rank()==0) {
1200 50 : all_nlist_size[walker_rank_]=nlist_index_.size();
1201 50 : multi_sim_comm.Sum(all_nlist_size);
1202 : }
1203 50 : comm.Bcast(all_nlist_size,0);
1204 : unsigned tot_size=0;
1205 150 : for(unsigned w=0; w<NumWalkers_; w++) {
1206 100 : tot_size+=all_nlist_size[w];
1207 : }
1208 50 : if(tot_size>0) {
1209 50 : std::vector<int> disp(NumWalkers_);
1210 100 : for(unsigned w=0; w<NumWalkers_-1; w++) {
1211 50 : disp[w+1]=disp[w]+all_nlist_size[w];
1212 : }
1213 50 : std::vector<unsigned> all_nlist_index(tot_size);
1214 50 : if(comm.Get_rank()==0) {
1215 50 : multi_sim_comm.Allgatherv(nlist_index_,all_nlist_index,&all_nlist_size[0],&disp[0]);
1216 : }
1217 50 : comm.Bcast(all_nlist_index,0);
1218 50 : std::set<unsigned> nlist_index_set(all_nlist_index.begin(),all_nlist_index.end()); //remove duplicates and sort
1219 50 : nlist_index_.assign(nlist_index_set.begin(),nlist_index_set.end());
1220 50 : sum_outdated_=true;
1221 : }
1222 : }
1223 378 : for(unsigned w=0; w<NumWalkers_; w++) {
1224 252 : std::vector<double> center_w(all_center.begin()+ncv_*w,all_center.begin()+ncv_*(w+1));
1225 252 : std::vector<double> sigma_w(all_sigma.begin()+ncv_*w,all_sigma.begin()+ncv_*(w+1));
1226 252 : addKernel(all_height[w],center_w,sigma_w,all_logweight[w]);
1227 : }
1228 : }
1229 13057 : getPntrToComponent("nker")->set(kernels_.size());
1230 13057 : if(nlist_) {
1231 6466 : getPntrToComponent("nlker")->set(nlist_index_.size());
1232 6466 : if(nlist_pace_reset_) {
1233 50 : nlist_update_=true;
1234 : }
1235 : }
1236 :
1237 : //update Zed_
1238 13057 : if(!no_Zed_) {
1239 12997 : double sum_uprob=0;
1240 12997 : const double ks=kernels_.size(); //in double, since ks*ks overflows an unsigned for more than 65535 kernels
1241 12997 : const double ds=delta_kernels_.size();
1242 12997 : const bool few_kernels=(ks*ks<(3*ks*ds+2*ds*ds*NumParallel_+100)); //this seems reasonable, but is not rigorous...
1243 12997 : if(few_kernels) { //really needed? Probably is almost always false
1244 4145 : #pragma omp parallel num_threads(NumOMP_)
1245 : {
1246 : #pragma omp for reduction(+:sum_uprob) nowait
1247 : for(unsigned k=rank_; k<kernels_.size(); k+=NumParallel_)
1248 : for(unsigned kk=0; kk<kernels_.size(); kk++) {
1249 : sum_uprob+=evaluateKernel(kernels_[kk],kernels_[k].center);
1250 : }
1251 : }
1252 4145 : if(NumParallel_>1) {
1253 50 : comm.Sum(sum_uprob);
1254 : }
1255 : } else {
1256 : // Here instead of redoing the full summation, we add only the changes, knowing that
1257 : // uprob = old_uprob + delta_uprob
1258 : // and we also need to consider that in the new sum there are some novel centers and some disappeared ones
1259 : // All kernels are used: the summed ones without the neighbor list, otherwise the copy used to rebuild the list,
1260 : // since a merged kernel changes the sum also at centers outside the list.
1261 : // Each MPI process takes a contiguous share of them, split in blocks among the OpenMP threads
1262 8852 : const bool all=nlist_;
1263 8852 : if(all && nlist_all_outdated_) {
1264 4447 : updateNlistAll();
1265 : }
1266 8852 : if(!all && sum_outdated_) {
1267 4405 : updateSum();
1268 : }
1269 8852 : const unsigned size=all?nlist_all_height_.size():sum_height_.size();
1270 8852 : const unsigned share=(size+NumParallel_-1)/NumParallel_;
1271 8852 : const unsigned begin=std::min(rank_*share,size);
1272 8852 : const unsigned end=std::min(begin+share,size);
1273 : constexpr unsigned omp_block=64;
1274 8852 : double delta_sum_uprob=0;
1275 8852 : if(NumOMP_==1 || end-begin<2*NumOMP_*omp_block) {
1276 4768 : delta_sum_uprob=deltaSumUprob(begin,end,all);
1277 : } else {
1278 4084 : const unsigned nblocks=(end-begin+omp_block-1)/omp_block;
1279 4084 : #pragma omp parallel num_threads(NumOMP_)
1280 : {
1281 : #pragma omp for reduction(+:delta_sum_uprob) nowait
1282 : for(unsigned n=0; n<nblocks; n++) {
1283 : const unsigned b=begin+n*omp_block;
1284 : delta_sum_uprob+=deltaSumUprob(b,std::min(b+omp_block,end),all);
1285 : }
1286 : }
1287 : }
1288 8852 : if(NumParallel_>1) {
1289 0 : comm.Sum(delta_sum_uprob);
1290 : }
1291 8852 : #pragma omp parallel num_threads(NumOMP_)
1292 : {
1293 : #pragma omp for reduction(+:delta_sum_uprob)
1294 : for(unsigned d=0; d<delta_kernels_.size(); d++) {
1295 : for(unsigned dd=0; dd<delta_kernels_.size(); dd++) {
1296 : //now subtract the delta_uprob added before, but not needed
1297 : const double sign=delta_kernels_[d].height<0?-1:1;
1298 : delta_sum_uprob-=sign*evaluateKernel(delta_kernels_[dd],delta_kernels_[d].center);
1299 : }
1300 : }
1301 : }
1302 8852 : sum_uprob=Zed_*old_KDEnorm_*old_nker+delta_sum_uprob;
1303 : }
1304 12997 : Zed_=sum_uprob/KDEnorm_/kernels_.size();
1305 25994 : getPntrToComponent("zed")->set(Zed_);
1306 : }
1307 :
1308 : //calculate work if requested
1309 13057 : if(calc_work_) {
1310 70 : std::vector<double> dummy(ncv_); //derivatives are not actually needed
1311 70 : const double prob=getProbAndDerivatives(center,dummy);
1312 70 : const double new_bias=kbt_*bias_prefactor_*std::log(prob/Zed_+epsilon_);
1313 70 : work_+=new_bias-getOutputQuantity(0);
1314 140 : getPntrToComponent("work")->set(work_);
1315 : }
1316 : }
1317 :
1318 : //dump state if requested
1319 64735 : if( (wStateStride_>0 && getStep()%wStateStride_==0) || (wStateStride_==-1 && getCPT()) ) {
1320 18 : dumpStateToFile();
1321 : }
1322 : }
1323 :
1324 : template <class mode>
1325 65317 : double OPESmetad<mode>::getProbAndDerivatives(const std::vector<double>& cv,std::vector<double>& der_prob) {
1326 65317 : if(sum_outdated_) {
1327 21768 : updateSum();
1328 : }
1329 : //each MPI process sums a contiguous share of the kernels, split in blocks among the OpenMP threads
1330 65317 : const unsigned size=sum_height_.size();
1331 65317 : const unsigned share=(size+NumParallel_-1)/NumParallel_;
1332 65317 : const unsigned begin=std::min(rank_*share,size);
1333 65317 : const unsigned end=std::min(begin+share,size);
1334 : constexpr unsigned omp_block=64;
1335 65317 : double prob=0.0;
1336 65317 : if(NumOMP_==1 || end-begin<2*NumOMP_*omp_block) {
1337 55109 : sumKernels(cv,begin,end,prob,der_prob);
1338 : } else {
1339 10208 : const unsigned nblocks=(end-begin+omp_block-1)/omp_block;
1340 10208 : #pragma omp parallel num_threads(NumOMP_)
1341 : {
1342 : double omp_prob=0;
1343 : std::vector<double> omp_deriv(ncv_,0.);
1344 : #pragma omp for nowait
1345 : for(unsigned n=0; n<nblocks; n++) {
1346 : const unsigned b=begin+n*omp_block;
1347 : sumKernels(cv,b,std::min(b+omp_block,end),omp_prob,omp_deriv);
1348 : }
1349 : #pragma omp critical
1350 : {
1351 : prob+=omp_prob;
1352 : for(unsigned i=0; i<ncv_; i++) {
1353 : der_prob[i]+=omp_deriv[i];
1354 : }
1355 : }
1356 : }
1357 : }
1358 65317 : if(NumParallel_>1) {
1359 224 : comm.Sum(prob);
1360 224 : comm.Sum(der_prob);
1361 : }
1362 : //normalize the estimate
1363 65317 : prob/=KDEnorm_;
1364 195765 : for(unsigned i=0; i<ncv_; i++) {
1365 130448 : der_prob[i]/=KDEnorm_;
1366 : }
1367 :
1368 65317 : return prob;
1369 : }
1370 :
1371 : template <class mode>
1372 26177 : void OPESmetad<mode>::updateSum() {
1373 26177 : const unsigned size=nlist_?nlist_index_.size():kernels_.size();
1374 26177 : sum_center_.resize(ncv_*size);
1375 26177 : sum_inv_sigma_.resize(ncv_*size);
1376 26177 : sum_height_.resize(size);
1377 1499515 : for(unsigned k=0; k<size; k++) {
1378 1473338 : const kernel& G=kernels_[nlist_?nlist_index_[k]:k];
1379 1473338 : sum_height_[k]=G.height;
1380 4419927 : for(unsigned i=0; i<ncv_; i++) {
1381 2946589 : sum_center_[i*size+k]=G.center[i];
1382 2946589 : sum_inv_sigma_[i*size+k]=1./G.sigma[i];
1383 : }
1384 : }
1385 26177 : sum_outdated_=false;
1386 26177 : }
1387 :
1388 : template <class mode>
1389 150460 : void OPESmetad<mode>::sumKernels(const std::vector<double>& cv,const unsigned begin,const unsigned end,double& prob,std::vector<double>& der_prob) const {
1390 : //adds the summed kernels from begin to end to prob and der_prob (not normalized).
1391 : //For a block of kernels, the squared distances are first computed for all of them, without branches,
1392 : //then only the few within the cutoff are evaluated. The result is the same as summing one kernel at a time
1393 150460 : constexpr unsigned block=128;
1394 : double norm2[block];
1395 150460 : const unsigned size=sum_height_.size();
1396 302249 : for(unsigned b=begin; b<end; b+=block) {
1397 151789 : const unsigned nb=std::min(block,end-b);
1398 7316928 : for(unsigned k=0; k<nb; k++) {
1399 7165139 : norm2[k]=0;
1400 : }
1401 455218 : for(unsigned i=0; i<ncv_; i++) {
1402 303429 : const double* center_i=&sum_center_[i*size+b];
1403 : const double* inv_sigma_i=&sum_inv_sigma_[i*size+b];
1404 303429 : const double cv_i=cv[i];
1405 303429 : if(period_[i]>0) {
1406 : const double period=period_[i];
1407 2277 : const double inv_period=inv_period_[i];
1408 : //unlike Tools::pbc, nearbyint has no branches: it is a single instruction on aarch64 and on x86 with SSE4.1
1409 12465 : for(unsigned k=0; k<nb; k++) {
1410 10188 : double diff=cv_i-center_i[k];
1411 10188 : diff-=period*std::nearbyint(diff*inv_period);
1412 10188 : const double dist=diff*inv_sigma_i[k];
1413 10188 : norm2[k]+=dist*dist;
1414 : }
1415 : } else {
1416 14620868 : for(unsigned k=0; k<nb; k++) {
1417 14319716 : const double dist=(cv_i-center_i[k])*inv_sigma_i[k];
1418 14319716 : norm2[k]+=dist*dist;
1419 : }
1420 : }
1421 : }
1422 7316928 : for(unsigned k=0; k<nb; k++) {
1423 7165139 : if(norm2[k]>=cutoff2_) {
1424 6890455 : continue;
1425 : }
1426 274684 : const double gaussian=sum_height_[b+k]*std::exp(-0.5*norm2[k]);
1427 274684 : prob+=gaussian-sum_height_[b+k]*val_at_cutoff_;
1428 823859 : for(unsigned i=0; i<ncv_; i++) {
1429 549175 : double diff=cv[i]-sum_center_[i*size+b+k];
1430 549175 : if(period_[i]>0) {
1431 4677 : diff-=period_[i]*std::nearbyint(diff*inv_period_[i]);
1432 : }
1433 549175 : const double inv_sigma=sum_inv_sigma_[i*size+b+k];
1434 : //NB: we accumulate the derivative into der_prob. The cutoff shift is constant and does not contribute to it
1435 549175 : der_prob[i]-=diff*inv_sigma*inv_sigma*gaussian;
1436 : }
1437 : }
1438 : }
1439 150460 : }
1440 :
1441 : template <class mode>
1442 42896 : double OPESmetad<mode>::deltaSumUprob(const unsigned begin,const unsigned end,const bool all) const {
1443 : //change of the sum of all kernels evaluated at all centers due to the delta_kernels_, for the kernels from
1444 : //begin to end: the delta kernels evaluated at the centers of these kernels, plus these kernels evaluated at the
1445 : //centers of the delta kernels (minus, for those that are gone). The kernels are the summed ones, or with all=true
1446 : //the nlist_all_ ones. As in sumKernels(), the squared distances of a block of kernels are computed without
1447 : //branches, both in units of the delta kernel sigma and of their own sigma
1448 42896 : const std::vector<double>& center=all?nlist_all_center_:sum_center_;
1449 42896 : const std::vector<double>& inv_sigma=all?nlist_all_inv_sigma_:sum_inv_sigma_;
1450 42896 : const std::vector<double>& height=all?nlist_all_height_:sum_height_;
1451 42896 : constexpr unsigned block=128;
1452 : double norm2_d[block]; //in units of the delta kernel sigma
1453 : double norm2_k[block]; //in units of their own sigma
1454 42896 : const unsigned size=height.size();
1455 : double delta_sum_uprob=0;
1456 86666 : for(unsigned b=begin; b<end; b+=block) {
1457 86660 : const unsigned nb=std::min(block,end-b);
1458 109516 : for(const kernel& D : delta_kernels_) {
1459 4159852 : for(unsigned k=0; k<nb; k++) {
1460 4094106 : norm2_d[k]=0;
1461 4094106 : norm2_k[k]=0;
1462 : }
1463 197238 : for(unsigned i=0; i<ncv_; i++) {
1464 131492 : const double* center_i=¢er[i*size+b];
1465 : const double* inv_sigma_i=&inv_sigma[i*size+b];
1466 131492 : const double c_i=D.center[i];
1467 131492 : const double inv_sigma_d=1./D.sigma[i];
1468 131492 : if(period_[i]>0) {
1469 : const double period=period_[i];
1470 424 : const double inv_period=inv_period_[i];
1471 9964 : for(unsigned k=0; k<nb; k++) {
1472 9540 : double diff=c_i-center_i[k];
1473 9540 : diff-=period*std::nearbyint(diff*inv_period);
1474 9540 : const double dist_d=diff*inv_sigma_d;
1475 9540 : const double dist_k=diff*inv_sigma_i[k];
1476 9540 : norm2_d[k]+=dist_d*dist_d;
1477 9540 : norm2_k[k]+=dist_k*dist_k;
1478 : }
1479 : } else {
1480 8309740 : for(unsigned k=0; k<nb; k++) {
1481 8178672 : const double diff=c_i-center_i[k];
1482 8178672 : const double dist_d=diff*inv_sigma_d;
1483 8178672 : const double dist_k=diff*inv_sigma_i[k];
1484 8178672 : norm2_d[k]+=dist_d*dist_d;
1485 8178672 : norm2_k[k]+=dist_k*dist_k;
1486 : }
1487 : }
1488 : }
1489 65746 : const double sign=D.height<0?-1:1; //take away contribution from kernels that are gone, and add the one from new ones
1490 4159852 : for(unsigned k=0; k<nb; k++) {
1491 4094106 : if(norm2_d[k]<cutoff2_) {
1492 83899 : delta_sum_uprob+=D.height*shiftedGaussian(norm2_d[k]);
1493 : }
1494 4094106 : if(norm2_k[k]<cutoff2_) {
1495 84257 : delta_sum_uprob+=sign*height[b+k]*shiftedGaussian(norm2_k[k]);
1496 : }
1497 : }
1498 : }
1499 : }
1500 42896 : return delta_sum_uprob;
1501 : }
1502 :
1503 : template <class mode>
1504 13355 : void OPESmetad<mode>::addKernel(const double height,const std::vector<double>& center,const std::vector<double>& sigma) {
1505 13355 : sum_outdated_=true;
1506 13355 : nlist_all_outdated_=true;
1507 : bool no_match=true;
1508 13355 : if(threshold2_!=0) {
1509 13355 : unsigned taker_k=getMergeableKernel(center,kernels_.size());
1510 13355 : if(taker_k<kernels_.size()) {
1511 : no_match=false;
1512 9462 : delta_kernels_.emplace_back(-1*kernels_[taker_k].height,kernels_[taker_k].center,kernels_[taker_k].sigma);
1513 9462 : mergeKernels(kernels_[taker_k],kernel(height,center,sigma));
1514 9462 : delta_kernels_.push_back(kernels_[taker_k]);
1515 9462 : if(recursive_merge_) { //the overhead is worth it if it keeps low the total number of kernels
1516 : //the neighbor list holds all kernels within the cutoff of the current CV, thus also those that can
1517 : //merge with the new kernel, but the merged kernel has moved: search all kernels
1518 9411 : unsigned giver_k=taker_k;
1519 9411 : taker_k=getMergeableKernel(kernels_[giver_k].center,giver_k,true);
1520 9788 : while(taker_k<kernels_.size()) {
1521 377 : delta_kernels_.pop_back();
1522 377 : delta_kernels_.emplace_back(-1*kernels_[taker_k].height,kernels_[taker_k].center,kernels_[taker_k].sigma);
1523 377 : if(taker_k>giver_k) { //saves time when erasing
1524 : std::swap(taker_k,giver_k);
1525 : }
1526 377 : mergeKernels(kernels_[taker_k],kernels_[giver_k]);
1527 377 : delta_kernels_.push_back(kernels_[taker_k]);
1528 377 : kernels_.erase(kernels_.begin()+giver_k);
1529 377 : if(nlist_) { //remove the giver from the sorted list, shift the indexes after it, and make sure the taker is listed
1530 : std::vector<unsigned> new_nlist_index;
1531 188 : new_nlist_index.reserve(nlist_index_.size()+1);
1532 2688 : for(unsigned k : nlist_index_) {
1533 2500 : if(k!=giver_k) {
1534 2313 : new_nlist_index.push_back(k<giver_k?k:k-1);
1535 : }
1536 : }
1537 188 : const auto taker_pos=std::lower_bound(new_nlist_index.begin(),new_nlist_index.end(),taker_k); //taker_k<giver_k
1538 188 : if(taker_pos==new_nlist_index.end() || *taker_pos!=taker_k) { //possible, since all kernels are searched
1539 2 : new_nlist_index.insert(taker_pos,taker_k);
1540 : }
1541 : nlist_index_.swap(new_nlist_index);
1542 : }
1543 377 : giver_k=taker_k;
1544 377 : taker_k=getMergeableKernel(kernels_[giver_k].center,giver_k,true);
1545 : }
1546 : }
1547 : }
1548 : }
1549 : if(no_match) {
1550 3893 : kernels_.emplace_back(height,center,sigma);
1551 3893 : delta_kernels_.emplace_back(height,center,sigma);
1552 3893 : if(nlist_) {
1553 1881 : nlist_index_.push_back(kernels_.size()-1);
1554 : }
1555 : }
1556 13355 : }
1557 :
1558 : template <class mode>
1559 13183 : void OPESmetad<mode>::addKernel(const double height,const std::vector<double>& center,const std::vector<double>& sigma,const double logweight) {
1560 13183 : addKernel(height,center,sigma);
1561 : //write to file
1562 13183 : kernelsOfile_.printField("time",getTime());
1563 39504 : for(unsigned i=0; i<ncv_; i++) {
1564 26321 : kernelsOfile_.printField(getPntrToArgument(i),center[i]);
1565 : }
1566 39504 : for(unsigned i=0; i<ncv_; i++) {
1567 52642 : kernelsOfile_.printField("sigma_"+getPntrToArgument(i)->getName(),sigma[i]);
1568 : }
1569 13183 : kernelsOfile_.printField("height",height);
1570 13183 : kernelsOfile_.printField("logweight",logweight);
1571 13183 : kernelsOfile_.printField();
1572 13183 : }
1573 :
1574 : template <class mode>
1575 23143 : unsigned OPESmetad<mode>::getMergeableKernel(const std::vector<double>& giver_center,const unsigned giver_k,const bool all_kernels) {
1576 : //returns kernels_.size() if no match is found
1577 : //with NLIST only the kernels in the list are searched, unless all_kernels is true
1578 23143 : const bool use_nlist=(nlist_ && !all_kernels);
1579 23143 : const unsigned size=use_nlist?nlist_index_.size():kernels_.size();
1580 23143 : unsigned min_k=kernels_.size();
1581 23143 : double min_norm2=threshold2_;
1582 23143 : #pragma omp parallel num_threads(NumOMP_)
1583 : {
1584 : unsigned min_k_omp = min_k;
1585 : double min_norm2_omp = threshold2_;
1586 : #pragma omp for nowait
1587 : for(unsigned nk=rank_; nk<size; nk+=NumParallel_) {
1588 : const unsigned k=use_nlist?nlist_index_[nk]:nk;
1589 : if(k==giver_k) { //a kernel should not be merged with itself
1590 : continue;
1591 : }
1592 : double norm2=0;
1593 : for(unsigned i=0; i<ncv_; i++) {
1594 : const double dist_i=difference(i,giver_center[i],kernels_[k].center[i])/kernels_[k].sigma[i];
1595 : norm2+=dist_i*dist_i;
1596 : if(norm2>=min_norm2_omp) {
1597 : break;
1598 : }
1599 : }
1600 : if(norm2<min_norm2_omp) {
1601 : min_norm2_omp=norm2;
1602 : min_k_omp=k;
1603 : }
1604 : }
1605 : #pragma omp critical
1606 : {
1607 : if(min_norm2_omp < min_norm2) {
1608 : min_norm2 = min_norm2_omp;
1609 : min_k = min_k_omp;
1610 : }
1611 : }
1612 : }
1613 23143 : if(NumParallel_>1) {
1614 134 : std::vector<double> all_min_norm2(NumParallel_);
1615 134 : std::vector<unsigned> all_min_k(NumParallel_);
1616 134 : comm.Allgather(min_norm2,all_min_norm2);
1617 134 : comm.Allgather(min_k,all_min_k);
1618 : const unsigned best=std::distance(std::begin(all_min_norm2),std::min_element(std::begin(all_min_norm2),std::end(all_min_norm2)));
1619 134 : min_k=all_min_k[best];
1620 : }
1621 23143 : return min_k;
1622 : }
1623 :
1624 : template <class mode>
1625 16665 : void OPESmetad<mode>::updateNlist(const std::vector<double>& new_center) {
1626 16665 : if(kernels_.size()==0) { //no need to check for neighbors
1627 210 : return;
1628 : }
1629 :
1630 16455 : nlist_center_=new_center;
1631 16455 : nlist_index_.clear();
1632 16455 : if(nlist_all_outdated_) {
1633 1923 : updateNlistAll();
1634 : }
1635 : //smallest sigma of all kernels: it bounds the distance of the kernels outside the list from below.
1636 : //Using only those outside would never trigger a rebuild when all kernels are in the list
1637 16455 : nlist_sigma_min_=nlist_all_sigma_min_;
1638 : //each MPI process lists a contiguous share of the kernels, split in blocks among the OpenMP threads
1639 16455 : const unsigned size=kernels_.size();
1640 16455 : const unsigned share=(size+NumParallel_-1)/NumParallel_;
1641 16455 : const unsigned begin=std::min(rank_*share,size);
1642 16455 : const unsigned end=std::min(begin+share,size);
1643 : constexpr unsigned omp_block=64;
1644 16455 : if(NumOMP_==1 || end-begin<2*NumOMP_*omp_block) {
1645 14938 : listKernels(begin,end,nlist_index_);
1646 : } else {
1647 1517 : const unsigned nblocks=(end-begin+omp_block-1)/omp_block;
1648 1517 : #pragma omp parallel num_threads(NumOMP_)
1649 : {
1650 : std::vector<unsigned> private_nlist_index;
1651 : #pragma omp for nowait
1652 : for(unsigned n=0; n<nblocks; n++) {
1653 : const unsigned b=begin+n*omp_block;
1654 : listKernels(b,std::min(b+omp_block,end),private_nlist_index);
1655 : }
1656 : #pragma omp critical
1657 : nlist_index_.insert(nlist_index_.end(),private_nlist_index.begin(),private_nlist_index.end());
1658 : }
1659 1517 : if(recursive_merge_) {
1660 1517 : std::sort(nlist_index_.begin(),nlist_index_.end());
1661 : }
1662 : }
1663 16455 : if(NumParallel_>1) {
1664 102 : std::vector<int> all_nlist_size(NumParallel_);
1665 102 : all_nlist_size[rank_]=nlist_index_.size();
1666 102 : comm.Sum(all_nlist_size);
1667 : unsigned tot_size=0;
1668 306 : for(unsigned r=0; r<NumParallel_; r++) {
1669 204 : tot_size+=all_nlist_size[r];
1670 : }
1671 102 : if(tot_size>0) {
1672 102 : std::vector<int> disp(NumParallel_);
1673 204 : for(unsigned r=0; r<NumParallel_-1; r++) {
1674 102 : disp[r+1]=disp[r]+all_nlist_size[r];
1675 : }
1676 102 : std::vector<unsigned> local_nlist_index=nlist_index_;
1677 102 : nlist_index_.resize(tot_size);
1678 102 : comm.Allgatherv(local_nlist_index,nlist_index_,&all_nlist_size[0],&disp[0]);
1679 102 : if(recursive_merge_) {
1680 102 : std::sort(nlist_index_.begin(),nlist_index_.end());
1681 : }
1682 : }
1683 : }
1684 16455 : getPntrToComponent("nlker")->set(nlist_index_.size());
1685 16455 : getPntrToComponent("nlsteps")->set(nlist_steps_);
1686 16455 : nlist_steps_=0;
1687 16455 : nlist_update_=false;
1688 16455 : sum_outdated_=true;
1689 : }
1690 :
1691 : template <class mode>
1692 6370 : void OPESmetad<mode>::updateNlistAll() {
1693 : //all kernels with each CV in a contiguous array, as for the sum, refreshed only when they changed
1694 6370 : const unsigned size=kernels_.size();
1695 6370 : nlist_all_center_.resize(ncv_*size);
1696 6370 : nlist_all_inv_sigma_.resize(ncv_*size);
1697 6370 : nlist_all_height_.resize(size);
1698 6370 : nlist_all_sigma_min_.assign(ncv_,std::numeric_limits<double>::max());
1699 1379736 : for(unsigned k=0; k<size; k++) {
1700 1373366 : nlist_all_height_[k]=kernels_[k].height;
1701 4120098 : for(unsigned i=0; i<ncv_; i++) {
1702 2746732 : nlist_all_center_[i*size+k]=kernels_[k].center[i];
1703 2746732 : nlist_all_inv_sigma_[i*size+k]=1./kernels_[k].sigma[i];
1704 2746732 : nlist_all_sigma_min_[i]=std::min(nlist_all_sigma_min_[i],kernels_[k].sigma[i]);
1705 : }
1706 : }
1707 6370 : nlist_all_outdated_=false;
1708 6370 : }
1709 :
1710 : template <class mode>
1711 30381 : void OPESmetad<mode>::listKernels(const unsigned begin,const unsigned end,std::vector<unsigned>& index) const {
1712 : //appends to index the kernels from begin to end within the neighbor list cutoff.
1713 : //As in sumKernels(), the squared distances of a block of kernels are computed without branches
1714 30381 : constexpr unsigned block=128;
1715 : double norm2[block];
1716 30381 : const unsigned size=kernels_.size();
1717 61022 : for(unsigned b=begin; b<end; b+=block) {
1718 30641 : const unsigned nb=std::min(block,end-b);
1719 1478164 : for(unsigned k=0; k<nb; k++) {
1720 1447523 : norm2[k]=0;
1721 : }
1722 91923 : for(unsigned i=0; i<ncv_; i++) {
1723 61282 : const double* center_i=&nlist_all_center_[i*size+b];
1724 : const double* inv_sigma_i=&nlist_all_inv_sigma_[i*size+b];
1725 61282 : const double c_i=nlist_center_[i];
1726 61282 : if(period_[i]>0) {
1727 : const double period=period_[i];
1728 488 : const double inv_period=inv_period_[i];
1729 6312 : for(unsigned k=0; k<nb; k++) {
1730 5824 : double diff=c_i-center_i[k];
1731 5824 : diff-=period*std::nearbyint(diff*inv_period);
1732 5824 : const double dist=diff*inv_sigma_i[k];
1733 5824 : norm2[k]+=dist*dist;
1734 : }
1735 : } else {
1736 2950016 : for(unsigned k=0; k<nb; k++) {
1737 2889222 : const double dist=(c_i-center_i[k])*inv_sigma_i[k];
1738 2889222 : norm2[k]+=dist*dist;
1739 : }
1740 : }
1741 : }
1742 1478164 : for(unsigned k=0; k<nb; k++) {
1743 1447523 : if(norm2[k]<=nlist_param_*cutoff2_) {
1744 49176 : index.push_back(b+k);
1745 : }
1746 : }
1747 : }
1748 30381 : }
1749 :
1750 : template <class mode>
1751 20 : void OPESmetad<mode>::dumpStateToFile() {
1752 : //gather adaptive sigma info if needed
1753 : //doing this while writing to file can lead to misterious slowdowns
1754 : std::vector<double> all_sigma0;
1755 : std::vector<double> all_av_cv;
1756 : std::vector<double> all_av_M2;
1757 20 : if(adaptive_sigma_ && NumWalkers_>1) {
1758 16 : all_sigma0.resize(NumWalkers_*ncv_);
1759 16 : all_av_cv.resize(NumWalkers_*ncv_);
1760 16 : all_av_M2.resize(NumWalkers_*ncv_);
1761 16 : if(comm.Get_rank()==0) {
1762 16 : multi_sim_comm.Allgather(sigma0_,all_sigma0);
1763 16 : multi_sim_comm.Allgather(av_cv_,all_av_cv);
1764 16 : multi_sim_comm.Allgather(av_M2_,all_av_M2);
1765 : }
1766 16 : comm.Bcast(all_sigma0,0);
1767 16 : comm.Bcast(all_av_cv,0);
1768 16 : comm.Bcast(all_av_M2,0);
1769 : }
1770 :
1771 : //only the first walker writes the state, the others would write on /dev/null,
1772 : //and their fields would pile up since they never rewind
1773 20 : if(walker_rank_>0) {
1774 : return;
1775 : }
1776 : //rewrite header or rewind file
1777 12 : if(storeOldStates_) {
1778 8 : stateOfile_.clearFields();
1779 : } else {
1780 4 : stateOfile_.rewind();
1781 : }
1782 : //define fields
1783 12 : stateOfile_.addConstantField("action");
1784 12 : stateOfile_.addConstantField("biasfactor");
1785 12 : stateOfile_.addConstantField("epsilon");
1786 12 : stateOfile_.addConstantField("kernel_cutoff");
1787 12 : stateOfile_.addConstantField("compression_threshold");
1788 12 : stateOfile_.addConstantField("zed");
1789 12 : stateOfile_.addConstantField("sum_weights");
1790 12 : stateOfile_.addConstantField("sum_weights2");
1791 12 : stateOfile_.addConstantField("counter");
1792 12 : if(adaptive_sigma_) {
1793 12 : stateOfile_.addConstantField("adaptive_counter");
1794 12 : if(NumWalkers_==1) {
1795 12 : for(unsigned i=0; i<ncv_; i++) {
1796 16 : stateOfile_.addConstantField("sigma0_"+getPntrToArgument(i)->getName());
1797 16 : stateOfile_.addConstantField("av_cv_"+getPntrToArgument(i)->getName());
1798 16 : stateOfile_.addConstantField("av_M2_"+getPntrToArgument(i)->getName());
1799 : }
1800 : } else {
1801 24 : for(unsigned w=0; w<NumWalkers_; w++)
1802 48 : for(unsigned i=0; i<ncv_; i++) {
1803 64 : const std::string arg_iw=getPntrToArgument(i)->getName()+"_"+std::to_string(w);
1804 32 : stateOfile_.addConstantField("sigma0_"+arg_iw);
1805 32 : stateOfile_.addConstantField("av_cv_"+arg_iw);
1806 64 : stateOfile_.addConstantField("av_M2_"+arg_iw);
1807 : }
1808 : }
1809 : }
1810 : //print fields
1811 36 : for(unsigned i=0; i<ncv_; i++) { //periodicity of CVs
1812 24 : stateOfile_.setupPrintValue(getPntrToArgument(i));
1813 : }
1814 24 : stateOfile_.printField("action",getName()+"_state");
1815 12 : stateOfile_.printField("biasfactor",biasfactor_);
1816 12 : stateOfile_.printField("epsilon",epsilon_);
1817 12 : stateOfile_.printField("kernel_cutoff",sqrt(cutoff2_));
1818 12 : stateOfile_.printField("compression_threshold",sqrt(threshold2_));
1819 12 : stateOfile_.printField("zed",Zed_);
1820 12 : stateOfile_.printField("sum_weights",sum_weights_);
1821 12 : stateOfile_.printField("sum_weights2",sum_weights2_);
1822 12 : stateOfile_.printField("counter",counter_);
1823 12 : if(adaptive_sigma_) {
1824 12 : stateOfile_.printField("adaptive_counter",adaptive_counter_);
1825 12 : if(NumWalkers_==1) {
1826 12 : for(unsigned i=0; i<ncv_; i++) {
1827 16 : stateOfile_.printField("sigma0_"+getPntrToArgument(i)->getName(),sigma0_[i]);
1828 16 : stateOfile_.printField("av_cv_"+getPntrToArgument(i)->getName(),av_cv_[i]);
1829 16 : stateOfile_.printField("av_M2_"+getPntrToArgument(i)->getName(),av_M2_[i]);
1830 : }
1831 : } else {
1832 24 : for(unsigned w=0; w<NumWalkers_; w++)
1833 48 : for(unsigned i=0; i<ncv_; i++) {
1834 64 : const std::string arg_iw=getPntrToArgument(i)->getName()+"_"+std::to_string(w);
1835 32 : stateOfile_.printField("sigma0_"+arg_iw,all_sigma0[w*ncv_+i]);
1836 32 : stateOfile_.printField("av_cv_"+arg_iw,all_av_cv[w*ncv_+i]);
1837 64 : stateOfile_.printField("av_M2_"+arg_iw,all_av_M2[w*ncv_+i]);
1838 : }
1839 : }
1840 : }
1841 : //print kernels
1842 63 : for(unsigned k=0; k<kernels_.size(); k++) {
1843 51 : stateOfile_.printField("time",getTime()); //this is not very usefull
1844 153 : for(unsigned i=0; i<ncv_; i++) {
1845 102 : stateOfile_.printField(getPntrToArgument(i),kernels_[k].center[i]);
1846 : }
1847 153 : for(unsigned i=0; i<ncv_; i++) {
1848 204 : stateOfile_.printField("sigma_"+getPntrToArgument(i)->getName(),kernels_[k].sigma[i]);
1849 : }
1850 51 : stateOfile_.printField("height",kernels_[k].height);
1851 51 : stateOfile_.printField();
1852 : }
1853 : //make sure file is written even if small
1854 12 : if(!storeOldStates_) {
1855 4 : stateOfile_.flush();
1856 : }
1857 : }
1858 :
1859 : template <class mode>
1860 130797 : inline double OPESmetad<mode>::evaluateKernel(const kernel& G,const std::vector<double>& x) const {
1861 : //NB: cannot be a method of kernel class, because uses external variables (for cutoff)
1862 : double norm2=0;
1863 253124 : for(unsigned i=0; i<ncv_; i++) {
1864 197603 : const double dist_i=difference(i,G.center[i],x[i])/G.sigma[i];
1865 197603 : norm2+=dist_i*dist_i;
1866 197603 : if(norm2>=cutoff2_) {
1867 : return 0;
1868 : }
1869 : }
1870 55521 : return G.height*shiftedGaussian(norm2);
1871 : }
1872 :
1873 : template <class mode>
1874 : inline double OPESmetad<mode>::shiftedGaussian(const double norm2) const {
1875 : //Gaussian at a squared distance norm2 within the cutoff, shifted to be zero at the cutoff
1876 223677 : return std::exp(-0.5*norm2)-val_at_cutoff_;
1877 : }
1878 :
1879 : template <class mode>
1880 9839 : inline void OPESmetad<mode>::mergeKernels(kernel& k1,const kernel& k2) {
1881 9839 : const double h=k1.height+k2.height;
1882 29492 : for(unsigned i=0; i<k1.center.size(); i++) {
1883 19653 : const bool isPeriodic_i=getPntrToArgument(i)->isPeriodic();
1884 19653 : if(isPeriodic_i) {
1885 917 : k1.center[i]=k2.center[i]+difference(i,k2.center[i],k1.center[i]); //fix PBC
1886 : }
1887 19653 : const double w1=k1.height/h;
1888 19653 : const double w2=k2.height/h;
1889 19653 : const double delta=k2.center[i]-k1.center[i];
1890 19653 : const double c_i=k1.center[i]+w2*delta;
1891 : // Merge centered second moments. Subtracting E[x]^2 from E[x^2]
1892 : // loses the width for narrow kernels far from the CV origin and can
1893 : // produce zero or negative variance from entirely finite input.
1894 19653 : const double ss_i=w1*k1.sigma[i]*k1.sigma[i]+
1895 19653 : w2*k2.sigma[i]*k2.sigma[i]+w1*w2*delta*delta;
1896 19653 : if(isPeriodic_i) {
1897 917 : k1.center[i]=bringBackInPbc(i,c_i);
1898 : } else {
1899 18736 : k1.center[i]=c_i;
1900 : }
1901 19653 : k1.sigma[i]=sqrt(ss_i);
1902 : }
1903 9839 : k1.height=h;
1904 9839 : }
1905 :
1906 : }
1907 : }
|
