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