Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017 of Glen Hocky and Andrew White
3 :
4 : The eds module is free software: you can redistribute it and/or modify
5 : it under the terms of the GNU Lesser General Public License as published by
6 : the Free Software Foundation, either version 3 of the License, or
7 : (at your option) any later version.
8 :
9 : The eds module is distributed in the hope that it will be useful,
10 : but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 : GNU Lesser General Public License for more details.
13 :
14 : You should have received a copy of the GNU Lesser General Public License
15 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
16 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
17 : #include "bias/Bias.h"
18 : #include "bias/ReweightBase.h"
19 : #include "core/ActionAtomistic.h"
20 : #include "core/ActionRegister.h"
21 : #include "core/Atoms.h"
22 : #include "core/PlumedMain.h"
23 : #include "tools/File.h"
24 : #include "tools/Matrix.h"
25 : #include "tools/Random.h"
26 :
27 : #include <iostream>
28 :
29 : using namespace PLMD;
30 : using namespace bias;
31 :
32 : // namespace is lowercase to match
33 : // module names being all lowercase
34 :
35 : namespace PLMD
36 : {
37 : namespace eds
38 : {
39 :
40 : //+PLUMEDOC EDSMOD_BIAS EDS
41 : /*
42 : Add a linear bias on a set of observables.
43 :
44 : This force is the same as the linear part of the bias in \ref
45 : RESTRAINT, but this bias has the ability to compute the prefactors
46 : adaptively using the scheme of White and Voth \cite white2014efficient
47 : in order to match target observable values for a set of CVs.
48 : Further updates to the algorithm are described in \cite hocky2017cgds
49 : and you can read a review on the method and its applications here: \cite Amirkulova2019Recent.
50 :
51 : You can
52 : see a tutorial on EDS specifically for biasing coordination number at
53 : <a
54 : href="http://thewhitelab.org/blog/tutorial/2017/05/10/lammps-coordination-number-tutorial/">
55 : Andrew White's webpage</a>.
56 :
57 : The addition to the potential is of the form
58 : \f[
59 : \sum_i \frac{\alpha_i}{s_i} x_i
60 : \f]
61 :
62 : where for CV \f$x_i\f$, a coupling constant \f${\alpha}_i\f$ is determined
63 : adaptively or set by the user to match a target value for
64 : \f$x_i\f$. \f$s_i\f$ is a scale parameter, which by default is set to
65 : the target value. It may also be set separately.
66 :
67 : \warning
68 : It is not possible to set the target value of the observable
69 : to zero with the default value of \f$s_i\f$ as this will cause a
70 : divide-by-zero error. Instead, set \f$s_i=1\f$ or modify the CV so the
71 : desired target value is no longer zero.
72 :
73 : Notice that a similar method is available as \ref MAXENT, although with different features and using a different optimization algorithm.
74 :
75 : \par Virial
76 :
77 : The bias forces modify the virial and this can change your simulation density if the bias is used in an NPT simulation.
78 : One way to avoid changing the virial contribution from the bias is to add the keyword VIRIAL=1.0, which changes how the bias
79 : is computed to minimize its contribution to the virial. This can also lead to smaller magnitude biases that are more robust if
80 : transferred to other systems. VIRIAL=1.0 can be a reasonable starting point and increasing the value changes the balance between matching
81 : the set-points and minimizing the virial. See \cite Amirkulova2019Recent for details on the equations. Since the coupling constants
82 : are unique with a single CV, VIRIAL is not applicable with a single CV. When used with multiple CVs, the CVs should be correlated
83 : which is almost always the case.
84 :
85 : \par Weighting
86 :
87 : EDS computes means and variances as part of its algorithm. If you are
88 : also using a biasing method like metadynamics, you may wish to remove
89 : the effect of this bias in your EDS computations so that EDS works on
90 : the canonical values (reweighted to be unbiased). For example, you may be using
91 : metadynamics to bias a dihedral angle to enhance sampling and be using
92 : EDS to set the average distance between two particular atoms. Specifically:
93 :
94 : \plumedfile
95 : # set-up metadynamics
96 : t: TORSION ATOMS=1,2,3,4
97 : md: METAD ARG=d SIGMA=0.2 HEIGHT=0.3 PACE=500 TEMP=300
98 : # compute bias weights
99 : bias: REWEIGHT_METAD TEMP=300
100 : # now do EDS on distance while removing effect of metadynamics
101 : d: DISTANCE ATOMS=4,7
102 : eds: EDS ARG=d CENTER=3.0 PERIOD=100 TEMP=300 LOGWEIGHTS=bias
103 : \endplumedfile
104 :
105 : This is an approximation though because EDS uses a finite samples while running to get means/variances.
106 : At the end of a run,
107 : you should ensure this approach worked and indeed your reweighted CV matches the target value.
108 :
109 : \par Examples
110 :
111 : The following input for a harmonic oscillator of two beads will
112 : adaptively find a linear bias to change the mean and variance to the
113 : target values. The PRINT line shows how to access the value of the
114 : coupling constants.
115 :
116 : \plumedfile
117 : dist: DISTANCE ATOMS=1,2
118 : # this is the squared of the distance
119 : dist2: COMBINE ARG=dist POWERS=2 PERIODIC=NO
120 :
121 : # bias mean and variance
122 : eds: EDS ARG=dist,dist2 CENTER=2.0,1.0 PERIOD=100 TEMP=1.0
123 : PRINT ARG=dist,dist2,eds.dist_coupling,eds.dist2_coupling,eds.bias,eds.force2 FILE=colvars.dat STRIDE=100
124 : \endplumedfile
125 :
126 : <hr>
127 :
128 : Rather than trying to find the coupling constants adaptively, one can ramp up to a constant value.
129 : \plumedfile
130 : dist: DISTANCE ATOMS=1,2
131 : dist2: COMBINE ARG=dist POWERS=2 PERIODIC=NO
132 :
133 : # ramp couplings from 0,0 to -1,1 over 50000 steps
134 : eds: EDS ARG=dist,dist2 CENTER=2.0,1.0 FIXED=-1,1 RAMP PERIOD=50000 TEMP=1.0
135 :
136 : # same as above, except starting at -0.5,0.5 rather than default of 0,0
137 : eds2: EDS ARG=dist,dist2 CENTER=2.0,1.0 FIXED=-1,1 INIT=-0.5,0.5 RAMP PERIOD=50000 TEMP=1.0
138 : \endplumedfile
139 :
140 : <hr>
141 : A restart file can be added to dump information needed to restart/continue simulation using these parameters every PERIOD.
142 : \plumedfile
143 : dist: DISTANCE ATOMS=1,2
144 : dist2: COMBINE ARG=dist POWERS=2 PERIODIC=NO
145 :
146 : # add the option to write to a restart file
147 : eds: EDS ARG=dist,dist2 CENTER=2.0,1.0 PERIOD=100 TEMP=1.0 OUT_RESTART=checkpoint.eds
148 : \endplumedfile
149 :
150 : The first few lines of the restart file that is output if we run a calculation with one CV will look something like this:
151 :
152 : \auxfile{restart.eds}
153 : #! FIELDS time d1_center d1_set d1_target d1_coupling d1_maxrange d1_maxgrad d1_accum d1_mean d1_std
154 : #! SET adaptive 1
155 : #! SET update_period 1
156 : #! SET seed 0
157 : #! SET kbt 2.4943
158 : 0.0000 1.0000 0.0000 0.0000 0.0000 7.4830 0.1497 0.0000 0.0000 0.0000
159 : 1.0000 1.0000 0.0000 0.0000 0.0000 7.4830 0.1497 0.0000 0.0000 0.0000
160 : 2.0000 1.0000 -7.4830 0.0000 0.0000 7.4830 0.1497 0.0224 0.0000 0.0000
161 : 3.0000 1.0000 -7.4830 0.0000 -7.4830 7.4830 0.1497 0.0224 0.0000 0.0000
162 : 4.0000 1.0000 -7.4830 0.0000 -7.4830 7.4830 0.1497 0.0224 0.0000 0.0000
163 : \endauxfile
164 :
165 : <hr>
166 :
167 : Read in a previous restart file. Adding RESTART flag makes output append
168 : \plumedfile
169 : d1: DISTANCE ATOMS=1,2
170 :
171 : eds: EDS ARG=d1 CENTER=2.0 PERIOD=100 TEMP=1.0 IN_RESTART=restart.eds RESTART=YES
172 : \endplumedfile
173 :
174 : <hr>
175 :
176 : Read in a previous restart file and freeze the bias at the final level from the previous simulation
177 : \plumedfile
178 : d1: DISTANCE ATOMS=1,2
179 :
180 : eds: EDS ARG=d1 CENTER=2.0 TEMP=1.0 IN_RESTART=restart.eds FREEZE
181 : \endplumedfile
182 :
183 : <hr>
184 :
185 : Read in a previous restart file and freeze the bias at the mean from the previous simulation
186 : \plumedfile
187 : d1: DISTANCE ATOMS=1,2
188 :
189 : eds: EDS ARG=d1 CENTER=2.0 TEMP=1.0 IN_RESTART=restart.eds FREEZE MEAN
190 : \endplumedfile
191 :
192 :
193 : */
194 : //+ENDPLUMEDOC
195 :
196 : class EDS : public Bias
197 : {
198 :
199 : private:
200 : /*We will get this and store it once, since on-the-fly changing number of CVs will be fatal*/
201 : const unsigned int ncvs_;
202 : std::vector<double> center_;
203 : std::vector<Value *> center_values_;
204 : ReweightBase *logweights_; // weights to use if reweighting averages
205 : std::vector<double> scale_;
206 : std::vector<double> current_coupling_; // actually current coupling
207 : std::vector<double> set_coupling_; // what our coupling is ramping up to. Equal to current_coupling when gathering stats
208 : std::vector<double> target_coupling_; // used when loaded to reach a value
209 : std::vector<double> max_coupling_range_; // used for adaptive range
210 : std::vector<double> max_coupling_grad_; // maximum allowed gradient
211 : std::vector<double> coupling_rate_;
212 : std::vector<double> coupling_accum_;
213 : std::vector<double> means_;
214 : std::vector<double> differences_;
215 : std::vector<double> alpha_vector_;
216 : std::vector<double> alpha_vector_2_;
217 : std::vector<double> ssds_;
218 : std::vector<double> step_size_;
219 : std::vector<double> pseudo_virial_;
220 : std::vector<Value *> out_coupling_;
221 : Matrix<double> covar_;
222 : Matrix<double> covar2_;
223 : Matrix<double> lm_inv_;
224 : std::string in_restart_name_;
225 : std::string out_restart_name_;
226 : std::string fmt_;
227 : OFile out_restart_;
228 : IFile in_restart_;
229 : bool b_c_values_;
230 : bool b_adaptive_;
231 : bool b_freeze_;
232 : bool b_equil_;
233 : bool b_ramp_;
234 : bool b_covar_;
235 : bool b_restart_;
236 : bool b_write_restart_;
237 : bool b_lm_;
238 : bool b_virial_;
239 : bool b_update_virial_;
240 : bool b_weights_;
241 : int seed_;
242 : int update_period_;
243 : int avg_coupling_count_;
244 : int update_calls_;
245 : double kbt_;
246 : double multi_prop_;
247 : double lm_mixing_par_;
248 : double virial_scaling_;
249 : double pseudo_virial_sum_; // net virial for all cvs in current period
250 : double max_logweight_; // maximum observed max logweight for period
251 : double wsum_; // sum of weights thus far
252 : Random rand_;
253 : Value *value_force2_;
254 : Value *value_pressure_;
255 :
256 : /*read input restart. b_mean sets if we use mean or final value for freeze*/
257 : void readInRestart(const bool b_mean);
258 : /*setup output restart*/
259 : void setupOutRestart();
260 : /*write output restart*/
261 : void writeOutRestart();
262 : void update_statistics();
263 : void update_pseudo_virial();
264 : void calc_lm_step_size();
265 : void calc_covar_step_size();
266 : void calc_ssd_step_size();
267 : void reset_statistics();
268 : void update_bias();
269 : void apply_bias();
270 :
271 : public:
272 : explicit EDS(const ActionOptions &);
273 : void calculate();
274 : void update();
275 : void turnOnDerivatives();
276 : static void registerKeywords(Keywords &keys);
277 : ~EDS();
278 : };
279 :
280 10435 : PLUMED_REGISTER_ACTION(EDS, "EDS")
281 :
282 9 : void EDS::registerKeywords(Keywords &keys)
283 : {
284 9 : Bias::registerKeywords(keys);
285 9 : keys.use("ARG");
286 18 : keys.add("optional", "CENTER", "The desired centers (equilibrium values) which will be sought during the adaptive linear biasing. This is for fixed centers");
287 18 : keys.add("optional", "CENTER_ARG", "The desired centers (equilibrium values) which will be sought during the adaptive linear biasing. "
288 : "CENTER_ARG is for calculated centers, e.g. from a CV or analysis. ");
289 :
290 18 : keys.add("optional", "PERIOD", "Steps over which to adjust bias for adaptive or ramping");
291 18 : keys.add("compulsory", "RANGE", "25.0", "The (starting) maximum increase in coupling constant per PERIOD (in \\f$k_B T\\f$/[BIAS_SCALE unit]) for each CV based");
292 18 : keys.add("compulsory", "SEED", "0", "Seed for random order of changing bias");
293 18 : keys.add("compulsory", "INIT", "0", "Starting value for coupling constant");
294 18 : keys.add("compulsory", "FIXED", "0", "Fixed target values for coupling constant. Non-adaptive.");
295 18 : keys.add("optional", "BIAS_SCALE", "A divisor to set the units of the bias. "
296 : "If not set, this will be the CENTER value by default (as is done in White and Voth 2014).");
297 18 : keys.add("optional", "TEMP", "The system temperature. If not provided will be taken from MD code (if available)");
298 18 : keys.add("optional", "MULTI_PROP", "What proportion of dimensions to update at each step. "
299 : "Must be in interval [1,0), where 1 indicates all and any other indicates a stochastic update. "
300 : "If not set, default is 1 / N, where N is the number of CVs. ");
301 18 : keys.add("optional", "VIRIAL", "Add an update penalty for having non-zero virial contributions. Only makes sense with multiple correlated CVs.");
302 18 : keys.add("optional", "LOGWEIGHTS", "Add weights to use for computing statistics. For example, if biasing with metadynamics.");
303 18 : keys.addFlag("LM", false, "Use Levenberg-Marquadt algorithm along with simultaneous keyword. Otherwise use gradient descent.");
304 18 : keys.add("compulsory", "LM_MIXING", "1", "Initial mixing parameter when using Levenberg-Marquadt minimization.");
305 18 : keys.add("optional", "RESTART_FMT", "the format that should be used to output real numbers in EDS restarts");
306 18 : keys.add("optional", "OUT_RESTART", "Output file for all information needed to continue EDS simulation. "
307 : "If you have the RESTART directive set (global or for EDS), this file will be appended to. "
308 : "Note that the header will be printed again if appending.");
309 18 : keys.add("optional", "IN_RESTART", "Read this file to continue an EDS simulation. "
310 : "If same as OUT_RESTART and you have not set the RESTART directive, the file will be backed-up and overwritten with new output. "
311 : "If you do have the RESTART flag set and it is the same name as OUT_RESTART, this file will be appended.");
312 :
313 18 : keys.addFlag("RAMP", false, "Slowly increase bias constant to a fixed value");
314 18 : keys.addFlag("COVAR", false, "Utilize the covariance matrix when updating the bias. Default Off, but may be enabled due to other options");
315 18 : keys.addFlag("FREEZE", false, "Fix bias at current level (only used for restarting).");
316 18 : keys.addFlag("MEAN", false, "Instead of using final bias level from restart, use average. Can only be used in conjunction with FREEZE");
317 :
318 9 : keys.use("RESTART");
319 :
320 18 : keys.addOutputComponent("force2", "default", "squared value of force from the bias");
321 18 : keys.addOutputComponent("pressure", "default", "If using virial keyword, this is the current sum of virial terms. It is in units of pressure (energy / vol^3)");
322 18 : keys.addOutputComponent("_coupling", "default", "For each named CV biased, there will be a corresponding output CV_coupling storing the current linear bias prefactor.");
323 9 : }
324 :
325 8 : EDS::EDS(const ActionOptions &ao) : PLUMED_BIAS_INIT(ao),
326 8 : ncvs_(getNumberOfArguments()),
327 8 : scale_(ncvs_, 0.0),
328 8 : current_coupling_(ncvs_, 0.0),
329 8 : set_coupling_(ncvs_, 0.0),
330 8 : target_coupling_(ncvs_, 0.0),
331 8 : max_coupling_range_(ncvs_, 25.0),
332 8 : max_coupling_grad_(ncvs_, 0.0),
333 8 : coupling_rate_(ncvs_, 1.0),
334 8 : coupling_accum_(ncvs_, 0.0),
335 8 : means_(ncvs_, 0.0),
336 8 : step_size_(ncvs_, 0.0),
337 8 : pseudo_virial_(ncvs_),
338 8 : out_coupling_(ncvs_, NULL),
339 8 : in_restart_name_(""),
340 8 : out_restart_name_(""),
341 8 : fmt_("%f"),
342 8 : b_adaptive_(true),
343 8 : b_freeze_(false),
344 8 : b_equil_(true),
345 8 : b_ramp_(false),
346 8 : b_covar_(false),
347 8 : b_restart_(false),
348 8 : b_write_restart_(false),
349 8 : b_lm_(false),
350 8 : b_virial_(false),
351 8 : b_weights_(false),
352 8 : seed_(0),
353 8 : update_period_(0),
354 8 : avg_coupling_count_(1),
355 8 : update_calls_(0),
356 8 : kbt_(0.0),
357 8 : multi_prop_(-1.0),
358 8 : lm_mixing_par_(0.1),
359 8 : virial_scaling_(0.),
360 8 : pseudo_virial_sum_(0.0),
361 8 : max_logweight_(0.0),
362 8 : wsum_(0.0),
363 32 : value_force2_(NULL)
364 : {
365 8 : double temp = -1.0;
366 8 : bool b_mean = false;
367 : std::vector<Value *> wvalues;
368 :
369 8 : addComponent("force2");
370 8 : componentIsNotPeriodic("force2");
371 8 : value_force2_ = getPntrToComponent("force2");
372 :
373 20 : for (unsigned int i = 0; i < ncvs_; ++i)
374 : {
375 12 : std::string comp = getPntrToArgument(i)->getName() + "_coupling";
376 12 : addComponent(comp);
377 12 : componentIsNotPeriodic(comp);
378 12 : out_coupling_[i] = getPntrToComponent(comp);
379 : }
380 :
381 8 : parseVector("CENTER", center_);
382 8 : parseArgumentList("CENTER_ARG", center_values_);
383 8 : parseArgumentList("LOGWEIGHTS", wvalues);
384 8 : parseVector("BIAS_SCALE", scale_);
385 8 : parseVector("RANGE", max_coupling_range_);
386 8 : parseVector("FIXED", target_coupling_);
387 8 : parseVector("INIT", set_coupling_);
388 8 : parse("PERIOD", update_period_);
389 8 : parse("TEMP", temp);
390 8 : parse("SEED", seed_);
391 8 : parse("MULTI_PROP", multi_prop_);
392 8 : parse("LM_MIXING", lm_mixing_par_);
393 8 : parse("RESTART_FMT", fmt_);
394 8 : parse("VIRIAL", virial_scaling_);
395 8 : fmt_ = " " + fmt_; // add space since parse strips them
396 8 : parse("OUT_RESTART", out_restart_name_);
397 8 : parseFlag("LM", b_lm_);
398 8 : parseFlag("RAMP", b_ramp_);
399 8 : parseFlag("FREEZE", b_freeze_);
400 8 : parseFlag("MEAN", b_mean);
401 8 : parseFlag("COVAR", b_covar_);
402 8 : parse("IN_RESTART", in_restart_name_);
403 8 : checkRead();
404 :
405 : /*
406 : * Things that are different when using changing centers:
407 : * 1. Scale
408 : * 2. The log file
409 : * 3. Reading Restarts
410 : */
411 :
412 8 : if (center_.size() == 0)
413 : {
414 1 : if (center_values_.size() == 0)
415 0 : error("Must set either CENTER or CENTER_ARG");
416 1 : else if (center_values_.size() != ncvs_)
417 0 : error("CENTER_ARG must contain the same number of variables as ARG");
418 1 : b_c_values_ = true;
419 1 : center_.resize(ncvs_);
420 1 : log.printf(" EDS will use possibly varying centers\n");
421 : }
422 : else
423 : {
424 7 : if (center_.size() != ncvs_)
425 0 : error("Must have same number of CENTER arguments as ARG arguments");
426 7 : else if (center_values_.size() != 0)
427 0 : error("You can only set CENTER or CENTER_ARG. Not both");
428 7 : b_c_values_ = false;
429 7 : log.printf(" EDS will use fixed centers\n");
430 : }
431 :
432 : // check for weights
433 8 : if (wvalues.size() > 1)
434 : {
435 0 : error("LOGWEIGHTS can only support one weight set. Please only pass one action");
436 : }
437 8 : else if (wvalues.size() == 1)
438 : {
439 1 : logweights_ = dynamic_cast<ReweightBase *>(wvalues[0]->getPntrToAction());
440 1 : b_weights_ = true;
441 : }
442 :
443 8 : log.printf(" setting scaling:");
444 8 : if (scale_.size() > 0 && scale_.size() < ncvs_)
445 : {
446 0 : error("the number of BIAS_SCALE values be the same as number of CVs");
447 : }
448 8 : else if (scale_.size() == 0 && b_c_values_)
449 : {
450 0 : log.printf(" Setting SCALE to be 1 for all CVs\n");
451 0 : scale_.resize(ncvs_);
452 0 : for (unsigned int i = 0; i < ncvs_; ++i)
453 0 : scale_[i] = 1;
454 : }
455 8 : else if (scale_.size() == 0 && !b_c_values_)
456 : {
457 2 : log.printf(" (default) ");
458 :
459 2 : scale_.resize(ncvs_);
460 6 : for (unsigned int i = 0; i < scale_.size(); ++i)
461 : {
462 4 : if (center_[i] == 0)
463 0 : error("BIAS_SCALE parameter has been set to CENTER value of 0 (as is default). This will divide by 0, so giving up. See doc for EDS bias");
464 4 : scale_[i] = center_[i];
465 : }
466 : }
467 : else
468 : {
469 14 : for (unsigned int i = 0; i < scale_.size(); ++i)
470 8 : log.printf(" %f", scale_[i]);
471 : }
472 8 : log.printf("\n");
473 :
474 8 : if (b_lm_)
475 : {
476 1 : log.printf(" EDS will perform Levenberg-Marquardt minimization with mixing parameter = %f\n", lm_mixing_par_);
477 1 : differences_.resize(ncvs_);
478 1 : alpha_vector_.resize(ncvs_);
479 1 : alpha_vector_2_.resize(ncvs_);
480 1 : covar_.resize(ncvs_, ncvs_);
481 1 : covar2_.resize(ncvs_, ncvs_);
482 1 : lm_inv_.resize(ncvs_, ncvs_);
483 1 : covar2_ *= 0;
484 1 : lm_inv_ *= 0;
485 1 : if (multi_prop_ != 1)
486 0 : log.printf(" WARNING - doing LM minimization but MULTI_PROP!=1\n");
487 : }
488 7 : else if (b_covar_)
489 : {
490 1 : log.printf(" EDS will utilize covariance matrix for update steps\n");
491 1 : covar_.resize(ncvs_, ncvs_);
492 : }
493 : else
494 : {
495 6 : log.printf(" EDS will utilize variance for update steps\n");
496 6 : ssds_.resize(ncvs_);
497 : }
498 :
499 8 : b_virial_ = virial_scaling_;
500 :
501 8 : if (b_virial_)
502 : {
503 1 : if (ncvs_ == 1)
504 0 : error("Minimizing the virial is only valid with multiply correlated collective variables.");
505 : // check that the CVs can be used to compute pseudo-virial
506 1 : log.printf(" EDS will compute virials of CVs and penalize with scale of %f. Checking CVs are valid...", virial_scaling_);
507 4 : for (unsigned int i = 0; i < ncvs_; ++i)
508 : {
509 3 : auto a = dynamic_cast<ActionAtomistic *>(getPntrToArgument(i)->getPntrToAction());
510 3 : if (!a)
511 0 : error("If using VIRIAL keyword, you must have normal CVs as arguments to EDS. Offending action: " + getPntrToArgument(i)->getPntrToAction()->getName());
512 : // cppcheck-suppress nullPointerRedundantCheck
513 3 : if (!(a->getPbc().isOrthorombic()))
514 3 : log.printf(" WARNING: EDS Virial should have a orthorombic cell\n");
515 : }
516 1 : log.printf("done.\n");
517 1 : addComponent("pressure");
518 1 : componentIsNotPeriodic("pressure");
519 1 : value_pressure_ = getPntrToComponent("pressure");
520 : }
521 :
522 8 : if (b_mean && !b_freeze_)
523 : {
524 0 : error("EDS keyword MEAN can only be used along with keyword FREEZE");
525 : }
526 :
527 8 : if (in_restart_name_ != "")
528 : {
529 2 : b_restart_ = true;
530 2 : log.printf(" reading simulation information from file: %s\n", in_restart_name_.c_str());
531 2 : readInRestart(b_mean);
532 : }
533 : else
534 : {
535 :
536 6 : if (temp >= 0.0)
537 6 : kbt_ = plumed.getAtoms().getKBoltzmann() * temp;
538 : else
539 0 : kbt_ = plumed.getAtoms().getKbT();
540 :
541 : // in driver, this results in kbt of 0
542 6 : if (kbt_ == 0)
543 : {
544 0 : error(" Unable to determine valid kBT. "
545 : "Could be because you are runnning from driver or MD didn't give temperature.\n"
546 : "Consider setting temperature manually with the TEMP keyword.");
547 : kbt_ = 1;
548 : }
549 :
550 6 : log.printf(" kBT = %f\n", kbt_);
551 6 : log.printf(" Updating every %i steps\n", update_period_);
552 :
553 6 : if (!b_c_values_)
554 : {
555 5 : log.printf(" with centers:");
556 14 : for (unsigned int i = 0; i < ncvs_; ++i)
557 : {
558 9 : log.printf(" %f ", center_[i]);
559 : }
560 : }
561 : else
562 : {
563 1 : log.printf(" with actions centers:");
564 2 : for (unsigned int i = 0; i < ncvs_; ++i)
565 : {
566 1 : log.printf(" %s ", center_values_[i]->getName().c_str());
567 : // add dependency on these actions
568 1 : addDependency(center_values_[i]->getPntrToAction());
569 : }
570 : }
571 :
572 6 : log.printf("\n with initial ranges / rates:\n");
573 16 : for (unsigned int i = 0; i < max_coupling_range_.size(); ++i)
574 : {
575 : // this is just an empirical guess. Bigger range, bigger grads. Less frequent updates, bigger changes
576 : //
577 : // using the current maxing out scheme, max_coupling_range is the biggest step that can be taken in any given interval
578 10 : max_coupling_range_[i] *= kbt_;
579 10 : max_coupling_grad_[i] = max_coupling_range_[i];
580 10 : log.printf(" %f / %f\n", max_coupling_range_[i], max_coupling_grad_[i]);
581 : }
582 :
583 6 : if (seed_ > 0)
584 : {
585 2 : log.printf(" setting random seed = %i", seed_);
586 2 : rand_.setSeed(seed_);
587 : }
588 :
589 16 : for (unsigned int i = 0; i < ncvs_; ++i)
590 10 : if (target_coupling_[i] != 0.0)
591 1 : b_adaptive_ = false;
592 :
593 6 : if (!b_adaptive_)
594 : {
595 1 : if (b_ramp_)
596 : {
597 1 : log.printf(" ramping up coupling constants over %i steps\n", update_period_);
598 : }
599 :
600 1 : log.printf(" with starting coupling constants");
601 2 : for (unsigned int i = 0; i < set_coupling_.size(); ++i)
602 1 : log.printf(" %f", set_coupling_[i]);
603 1 : log.printf("\n");
604 1 : log.printf(" and final coupling constants");
605 2 : for (unsigned int i = 0; i < target_coupling_.size(); ++i)
606 1 : log.printf(" %f", target_coupling_[i]);
607 1 : log.printf("\n");
608 : }
609 :
610 : // now do setup
611 6 : if (b_ramp_)
612 : {
613 1 : update_period_ *= -1;
614 : }
615 :
616 16 : for (unsigned int i = 0; i < set_coupling_.size(); ++i)
617 10 : current_coupling_[i] = set_coupling_[i];
618 :
619 : // if b_adaptive_, then first half will be used for equilibrating and second half for statistics
620 6 : if (update_period_ > 0)
621 : {
622 5 : update_period_ /= 2;
623 : }
624 : }
625 :
626 8 : if (b_freeze_)
627 : {
628 1 : b_adaptive_ = false;
629 1 : update_period_ = 0;
630 1 : if (b_mean)
631 : {
632 1 : log.printf(" freezing bias at the average level from the restart file\n");
633 : }
634 : else
635 : {
636 0 : log.printf(" freezing bias at current level\n");
637 : }
638 : }
639 :
640 8 : if (multi_prop_ == -1.0)
641 : {
642 5 : log.printf(" Will update each dimension stochastically with probability 1 / number of CVs\n");
643 5 : multi_prop_ = 1.0 / ncvs_;
644 : }
645 3 : else if (multi_prop_ > 0 && multi_prop_ <= 1.0)
646 : {
647 3 : log.printf(" Will update each dimension stochastically with probability %f\n", multi_prop_);
648 : }
649 : else
650 : {
651 0 : error(" MULTI_PROP must be between 0 and 1\n");
652 : }
653 :
654 8 : if (out_restart_name_.length() > 0)
655 : {
656 8 : log.printf(" writing restart information every %i steps to file %s with format %s\n", abs(update_period_), out_restart_name_.c_str(), fmt_.c_str());
657 8 : b_write_restart_ = true;
658 8 : setupOutRestart();
659 : }
660 :
661 16 : log << " Bibliography " << plumed.cite("White and Voth, J. Chem. Theory Comput. 10 (8), 3023-3030 (2014)") << "\n";
662 24 : log << " Bibliography " << plumed.cite("G. M. Hocky, T. Dannenhoffer-Lafage, G. A. Voth, J. Chem. Theory Comput. 13 (9), 4593-4603 (2017)") << "\n";
663 8 : }
664 :
665 2 : void EDS::readInRestart(const bool b_mean)
666 : {
667 2 : int adaptive_i = 0;
668 :
669 2 : in_restart_.open(in_restart_name_);
670 :
671 4 : if (in_restart_.FieldExist("kbt"))
672 : {
673 2 : in_restart_.scanField("kbt", kbt_);
674 : }
675 : else
676 : {
677 0 : error("No field 'kbt' in restart file");
678 : }
679 2 : log.printf(" with kBT = %f\n", kbt_);
680 :
681 4 : if (in_restart_.FieldExist("update_period"))
682 : {
683 2 : in_restart_.scanField("update_period", update_period_);
684 : }
685 : else
686 : {
687 0 : error("No field 'update_period' in restart file");
688 : }
689 2 : log.printf(" Updating every %i steps\n", update_period_);
690 :
691 4 : if (in_restart_.FieldExist("adaptive"))
692 : {
693 : // note, no version of scanField for boolean
694 2 : in_restart_.scanField("adaptive", adaptive_i);
695 : }
696 : else
697 : {
698 0 : error("No field 'adaptive' in restart file");
699 : }
700 2 : b_adaptive_ = bool(adaptive_i);
701 :
702 4 : if (in_restart_.FieldExist("seed"))
703 : {
704 2 : in_restart_.scanField("seed", seed_);
705 : }
706 : else
707 : {
708 0 : error("No field 'seed' in restart file");
709 : }
710 2 : if (seed_ > 0)
711 : {
712 0 : log.printf(" setting random seed = %i", seed_);
713 0 : rand_.setSeed(seed_);
714 : }
715 :
716 : double time, tmp;
717 2 : std::vector<double> avg_bias = std::vector<double>(center_.size());
718 : unsigned int N = 0;
719 : std::string cv_name;
720 :
721 24 : while (in_restart_.scanField("time", time))
722 : {
723 :
724 20 : for (unsigned int i = 0; i < ncvs_; ++i)
725 : {
726 : cv_name = getPntrToArgument(i)->getName();
727 20 : in_restart_.scanField(cv_name + "_center", set_coupling_[i]);
728 20 : in_restart_.scanField(cv_name + "_set", set_coupling_[i]);
729 20 : in_restart_.scanField(cv_name + "_target", target_coupling_[i]);
730 20 : in_restart_.scanField(cv_name + "_coupling", current_coupling_[i]);
731 20 : in_restart_.scanField(cv_name + "_maxrange", max_coupling_range_[i]);
732 20 : in_restart_.scanField(cv_name + "_maxgrad", max_coupling_grad_[i]);
733 20 : in_restart_.scanField(cv_name + "_accum", coupling_accum_[i]);
734 10 : in_restart_.scanField(cv_name + "_mean", means_[i]);
735 20 : if (in_restart_.FieldExist(cv_name + "_pseudovirial"))
736 : {
737 0 : if (b_virial_)
738 0 : in_restart_.scanField(cv_name + "_pseudovirial", pseudo_virial_[i]);
739 : else // discard the field
740 0 : in_restart_.scanField(cv_name + "_pseudovirial", tmp);
741 : }
742 : // unused due to difference between covar/nocovar
743 20 : in_restart_.scanField(cv_name + "_std", tmp);
744 :
745 10 : avg_bias[i] += current_coupling_[i];
746 : }
747 10 : N++;
748 :
749 10 : in_restart_.scanField();
750 : }
751 :
752 2 : log.printf(" with centers:");
753 4 : for (unsigned int i = 0; i < center_.size(); ++i)
754 : {
755 2 : log.printf(" %f", center_[i]);
756 : }
757 2 : log.printf("\n and scaling:");
758 4 : for (unsigned int i = 0; i < scale_.size(); ++i)
759 : {
760 2 : log.printf(" %f", scale_[i]);
761 : }
762 :
763 2 : log.printf("\n with initial ranges / rates:\n");
764 4 : for (unsigned int i = 0; i < max_coupling_range_.size(); ++i)
765 : {
766 2 : log.printf(" %f / %f\n", max_coupling_range_[i], max_coupling_grad_[i]);
767 : }
768 :
769 2 : if (!b_adaptive_ && update_period_ < 0)
770 : {
771 0 : log.printf(" ramping up coupling constants over %i steps\n", -update_period_);
772 : }
773 :
774 2 : if (b_mean)
775 : {
776 1 : log.printf("Loaded in averages for coupling constants...\n");
777 2 : for (unsigned int i = 0; i < current_coupling_.size(); ++i)
778 1 : current_coupling_[i] = avg_bias[i] / N;
779 2 : for (unsigned int i = 0; i < current_coupling_.size(); ++i)
780 1 : set_coupling_[i] = avg_bias[i] / N;
781 : }
782 :
783 2 : log.printf(" with current coupling constants:\n ");
784 4 : for (unsigned int i = 0; i < current_coupling_.size(); ++i)
785 2 : log.printf(" %f", current_coupling_[i]);
786 2 : log.printf("\n");
787 2 : log.printf(" with initial coupling constants:\n ");
788 4 : for (unsigned int i = 0; i < set_coupling_.size(); ++i)
789 2 : log.printf(" %f", set_coupling_[i]);
790 2 : log.printf("\n");
791 2 : log.printf(" and final coupling constants:\n ");
792 4 : for (unsigned int i = 0; i < target_coupling_.size(); ++i)
793 2 : log.printf(" %f", target_coupling_[i]);
794 2 : log.printf("\n");
795 :
796 2 : in_restart_.close();
797 2 : }
798 :
799 8 : void EDS::setupOutRestart()
800 : {
801 8 : out_restart_.link(*this);
802 8 : out_restart_.fmtField(fmt_);
803 8 : out_restart_.open(out_restart_name_);
804 : out_restart_.setHeavyFlush();
805 :
806 16 : out_restart_.addConstantField("adaptive").printField("adaptive", b_adaptive_);
807 16 : out_restart_.addConstantField("update_period").printField("update_period", update_period_);
808 16 : out_restart_.addConstantField("seed").printField("seed", seed_);
809 16 : out_restart_.addConstantField("kbt").printField("kbt", kbt_);
810 8 : }
811 :
812 27 : void EDS::writeOutRestart()
813 : {
814 : std::string cv_name;
815 27 : out_restart_.printField("time", getTimeStep() * getStep());
816 :
817 66 : for (unsigned int i = 0; i < ncvs_; ++i)
818 : {
819 : cv_name = getPntrToArgument(i)->getName();
820 78 : out_restart_.printField(cv_name + "_center", center_[i]);
821 78 : out_restart_.printField(cv_name + "_set", set_coupling_[i]);
822 78 : out_restart_.printField(cv_name + "_target", target_coupling_[i]);
823 78 : out_restart_.printField(cv_name + "_coupling", current_coupling_[i]);
824 78 : out_restart_.printField(cv_name + "_maxrange", max_coupling_range_[i]);
825 78 : out_restart_.printField(cv_name + "_maxgrad", max_coupling_grad_[i]);
826 78 : out_restart_.printField(cv_name + "_accum", coupling_accum_[i]);
827 39 : out_restart_.printField(cv_name + "_mean", means_[i]);
828 39 : if (b_virial_)
829 18 : out_restart_.printField(cv_name + "_pseudovirial", pseudo_virial_[i]);
830 39 : if (!b_covar_ && !b_lm_)
831 42 : out_restart_.printField(cv_name + "_std", ssds_[i] / (fmax(1, update_calls_ - 1)));
832 : else
833 36 : out_restart_.printField(cv_name + "_std", covar_(i, i) / (fmax(1, update_calls_ - 1)));
834 : }
835 27 : out_restart_.printField();
836 27 : }
837 :
838 40 : void EDS::calculate()
839 : {
840 :
841 : // get center values from action if necessary
842 40 : if (b_c_values_)
843 10 : for (unsigned int i = 0; i < ncvs_; ++i)
844 5 : center_[i] = center_values_[i]->get();
845 :
846 40 : apply_bias();
847 40 : }
848 :
849 40 : void EDS::apply_bias()
850 : {
851 : // Compute linear force as in "restraint"
852 : double ene = 0, totf2 = 0, cv, m, f;
853 :
854 100 : for (unsigned int i = 0; i < ncvs_; ++i)
855 : {
856 60 : cv = difference(i, center_[i], getArgument(i));
857 60 : m = current_coupling_[i];
858 60 : f = -m;
859 60 : ene += m * cv;
860 : setOutputForce(i, f);
861 60 : totf2 += f * f;
862 : }
863 :
864 : setBias(ene);
865 40 : value_force2_->set(totf2);
866 40 : }
867 :
868 12 : void EDS::update_statistics()
869 : {
870 : double s, N, w = 1.0;
871 12 : std::vector<double> deltas(ncvs_);
872 :
873 : // update weight max, if necessary
874 12 : if (b_weights_)
875 : {
876 2 : w = logweights_->getLogWeight();
877 2 : if (max_logweight_ < w)
878 : {
879 : // we have new max. Need to shift existing values
880 0 : wsum_ *= exp(max_logweight_ - w);
881 0 : max_logweight_ = w;
882 : }
883 : // convert to weight
884 2 : w = exp(w - max_logweight_);
885 2 : wsum_ += w;
886 : N = wsum_;
887 : }
888 : else
889 : {
890 10 : N = fmax(1, update_calls_);
891 : }
892 :
893 : // Welford, West, and Hanso online variance method
894 : // with weights (default = 1.0)
895 32 : for (unsigned int i = 0; i < ncvs_; ++i)
896 : {
897 20 : deltas[i] = difference(i, means_[i], getArgument(i)) * w;
898 20 : means_[i] += deltas[i] / N;
899 20 : if (!b_covar_ && !b_lm_)
900 8 : ssds_[i] += deltas[i] * difference(i, means_[i], getArgument(i));
901 : }
902 12 : if (b_covar_ || b_lm_)
903 : {
904 16 : for (unsigned int i = 0; i < ncvs_; ++i)
905 : {
906 36 : for (unsigned int j = i; j < ncvs_; ++j)
907 : {
908 24 : s = (N - 1) * deltas[i] * deltas[j] / N / N - covar_(i, j) / N;
909 24 : covar_(i, j) += s;
910 : // do this so we don't double count
911 24 : covar_(j, i) = covar_(i, j);
912 : }
913 : }
914 : }
915 12 : if (b_virial_)
916 2 : update_pseudo_virial();
917 12 : }
918 :
919 8 : void EDS::reset_statistics()
920 : {
921 20 : for (unsigned int i = 0; i < ncvs_; ++i)
922 : {
923 12 : means_[i] = 0;
924 12 : if (!b_covar_ && !b_lm_)
925 6 : ssds_[i] = 0;
926 : }
927 8 : if (b_covar_ || b_lm_)
928 8 : for (unsigned int i = 0; i < ncvs_; ++i)
929 24 : for (unsigned int j = 0; j < ncvs_; ++j)
930 18 : covar_(i, j) = 0;
931 8 : if (b_virial_)
932 : {
933 4 : for (unsigned int i = 0; i < ncvs_; ++i)
934 3 : pseudo_virial_[i] = 0;
935 1 : pseudo_virial_sum_ = 0;
936 : }
937 8 : if (b_weights_)
938 : {
939 2 : wsum_ = 0;
940 2 : max_logweight_ = 0;
941 : }
942 8 : }
943 :
944 1 : void EDS::calc_lm_step_size()
945 : {
946 : // calulcate step size
947 : // uses scale here, which by default is center
948 :
949 1 : mult(covar_, covar_, covar2_);
950 4 : for (unsigned int i = 0; i < ncvs_; ++i)
951 : {
952 3 : differences_[i] = difference(i, center_[i], means_[i]);
953 3 : covar2_[i][i] += lm_mixing_par_ * covar2_[i][i];
954 : }
955 :
956 : // "step_size_vec" = 2*inv(covar*covar+ lambda diag(covar*covar))*covar*(mean-center)
957 1 : mult(covar_, differences_, alpha_vector_);
958 1 : Invert(covar2_, lm_inv_);
959 1 : mult(lm_inv_, alpha_vector_, alpha_vector_2_);
960 :
961 4 : for (unsigned int i = 0; i < ncvs_; ++i)
962 : {
963 3 : step_size_[i] = 2 * alpha_vector_2_[i] / kbt_ / scale_[i];
964 : }
965 1 : }
966 :
967 1 : void EDS::calc_covar_step_size()
968 : {
969 : // calulcate step size
970 : // uses scale here, which by default is center
971 : double tmp;
972 4 : for (unsigned int i = 0; i < ncvs_; ++i)
973 : {
974 : tmp = 0;
975 12 : for (unsigned int j = 0; j < ncvs_; ++j)
976 9 : tmp += difference(i, center_[i], means_[i]) * covar_(i, j);
977 3 : step_size_[i] = 2 * tmp / kbt_ / scale_[i] * update_calls_ / fmax(1, update_calls_ - 1);
978 : }
979 1 : }
980 :
981 6 : void EDS::calc_ssd_step_size()
982 : {
983 : double tmp;
984 12 : for (unsigned int i = 0; i < ncvs_; ++i)
985 : {
986 6 : tmp = 2. * difference(i, center_[i], means_[i]) * ssds_[i] / fmax(1, update_calls_ - 1);
987 6 : step_size_[i] = tmp / kbt_ / scale_[i];
988 : }
989 6 : }
990 :
991 2 : void EDS::update_pseudo_virial()
992 : {
993 : // We want to compute the bias force on each atom times the position
994 : // of the atoms.
995 : double p, netp = 0, netpv = 0;
996 : double volume = 0;
997 8 : for (unsigned int i = 0; i < ncvs_; ++i)
998 : {
999 : // checked in setup to ensure this cast is valid.
1000 6 : ActionAtomistic *cv = dynamic_cast<ActionAtomistic *>(getPntrToArgument(i)->getPntrToAction());
1001 : Tensor &v(cv->modifyVirial());
1002 6 : Tensor box(cv->getBox());
1003 : const unsigned int natoms = cv->getNumberOfAtoms();
1004 6 : if (!volume)
1005 2 : volume = box.determinant();
1006 :
1007 : // pressure contribution is -dBias / dV
1008 : // dBias / dV = alpha / w * dCV / dV
1009 : // to get partial of CV wrt to volume
1010 : // dCV/dV = sum dCV/dvij * vij / V
1011 : // where vij is box element
1012 : // add diagonal of virial tensor to get net pressure
1013 : // TODO: replace this with adjugate (Jacobi's Formula) for non-orthorombic case(?)
1014 6 : p = v(0, 0) * box(0, 0) + v(1, 1) * box(1, 1) + v(2, 2) * box(2, 2);
1015 6 : p /= volume;
1016 :
1017 6 : netp += p;
1018 :
1019 : // now scale for correct units in EDS algorithm
1020 6 : p *= (volume) / (kbt_ * natoms);
1021 :
1022 : // compute running mean of scaled
1023 6 : if (set_coupling_[i] != 0)
1024 0 : pseudo_virial_[i] = (p - pseudo_virial_[i]) / (fmax(1, update_calls_));
1025 : else
1026 6 : pseudo_virial_[i] = 0;
1027 : // update net pressure
1028 6 : netpv += pseudo_virial_[i];
1029 : }
1030 : // update pressure
1031 2 : value_pressure_->set(netp);
1032 2 : pseudo_virial_sum_ = netpv;
1033 2 : }
1034 :
1035 8 : void EDS::update_bias()
1036 : {
1037 8 : log.flush();
1038 8 : if (b_lm_)
1039 1 : calc_lm_step_size();
1040 7 : else if (b_covar_)
1041 1 : calc_covar_step_size();
1042 : else
1043 6 : calc_ssd_step_size();
1044 :
1045 20 : for (unsigned int i = 0; i < ncvs_; ++i)
1046 : {
1047 :
1048 : // multidimesional stochastic step
1049 12 : if (ncvs_ == 1 || (rand_.RandU01() < (multi_prop_)))
1050 : {
1051 :
1052 12 : if (b_virial_)
1053 : {
1054 : // apply virial regularization
1055 : // P * dP/dcoupling
1056 : // coupling is already included in virial term due to plumed propogating from bias to forces
1057 : // thus we need to divide by it to get the derivative (since force is linear in coupling)
1058 3 : if (fabs(set_coupling_[i]) > 0.000000001) // my heuristic for if EDS has started to prevent / 0
1059 : // scale^2 here is to align units
1060 0 : step_size_[i] -= 2 * scale_[i] * scale_[i] * virial_scaling_ * pseudo_virial_sum_ * pseudo_virial_sum_ / set_coupling_[i];
1061 : }
1062 12 : if (step_size_[i] == 0)
1063 4 : continue;
1064 :
1065 : // clip gradient
1066 8 : step_size_[i] = copysign(fmin(fabs(step_size_[i]), max_coupling_grad_[i]), step_size_[i]);
1067 8 : coupling_accum_[i] += step_size_[i] * step_size_[i];
1068 :
1069 : // equation 5 in White and Voth, JCTC 2014
1070 : // no negative sign because it's in step_size
1071 8 : set_coupling_[i] += step_size_[i] * max_coupling_range_[i] / sqrt(coupling_accum_[i]);
1072 8 : coupling_rate_[i] = (set_coupling_[i] - current_coupling_[i]) / update_period_;
1073 : }
1074 : else
1075 : {
1076 : // do not change the bias
1077 0 : coupling_rate_[i] = 0;
1078 : }
1079 : }
1080 :
1081 : // reset means/vars
1082 8 : reset_statistics();
1083 8 : }
1084 :
1085 40 : void EDS::update()
1086 : {
1087 : // adjust parameters according to EDS recipe
1088 40 : update_calls_++;
1089 :
1090 : // if we aren't wating for the bias to equilibrate, set flag to collect data
1091 : // want statistics before writing restart
1092 40 : if (!b_equil_ && update_period_ > 0)
1093 12 : update_statistics();
1094 :
1095 : // write restart with correct statistics before bias update
1096 : // check if we're ramping or doing normal updates and then restart if needed. The ramping check
1097 : // is complicated because we could be frozen, finished ramping or not ramping.
1098 : // The + 2 is so we have an extra line showing that the bias isn't changing (for my sanity and yours)
1099 40 : if (b_write_restart_)
1100 : {
1101 40 : if (getStep() == 0 ||
1102 32 : ((update_period_ < 0 && !b_freeze_ && update_calls_ <= fabs(update_period_) + 2) ||
1103 24 : (update_period_ > 0 && update_calls_ % update_period_ == 0)))
1104 27 : writeOutRestart();
1105 : }
1106 :
1107 : int b_finished_equil_flag = 1;
1108 :
1109 : // assume forces already applied and saved
1110 : // are we ramping to a constant value and not done equilibrating?
1111 40 : if (update_period_ < 0)
1112 : {
1113 5 : if (update_calls_ <= fabs(update_period_) && !b_freeze_)
1114 : {
1115 4 : for (unsigned int i = 0; i < ncvs_; ++i)
1116 2 : current_coupling_[i] += (target_coupling_[i] - set_coupling_[i]) / fabs(update_period_);
1117 : }
1118 : // make sure we don't reset update calls
1119 : b_finished_equil_flag = 0;
1120 : }
1121 35 : else if (update_period_ == 0)
1122 : { // do we have a no-update case?
1123 : // not updating
1124 : // pass
1125 : }
1126 30 : else if (b_equil_)
1127 : {
1128 : // equilibrating
1129 : // check if we've reached the setpoint
1130 48 : for (unsigned int i = 0; i < ncvs_; ++i)
1131 : {
1132 30 : if (coupling_rate_[i] == 0 || pow(current_coupling_[i] - set_coupling_[i], 2) < pow(coupling_rate_[i], 2))
1133 : {
1134 14 : b_finished_equil_flag &= 1;
1135 : }
1136 : else
1137 : {
1138 16 : current_coupling_[i] += coupling_rate_[i];
1139 : b_finished_equil_flag = 0;
1140 : }
1141 : }
1142 : }
1143 :
1144 : // reduce all the flags
1145 40 : if (b_equil_ && b_finished_equil_flag)
1146 : {
1147 11 : b_equil_ = false;
1148 11 : update_calls_ = 0;
1149 : }
1150 :
1151 : // Now we update coupling constant, if necessary
1152 40 : if (!b_equil_ && update_period_ > 0 && update_calls_ == update_period_ && !b_freeze_)
1153 : {
1154 8 : update_bias();
1155 8 : update_calls_ = 0;
1156 8 : avg_coupling_count_++;
1157 8 : b_equil_ = true; // back to equilibration now
1158 : } // close update if
1159 :
1160 : // pass couplings out so they are accessible
1161 100 : for (unsigned int i = 0; i < ncvs_; ++i)
1162 : {
1163 60 : out_coupling_[i]->set(current_coupling_[i]);
1164 : }
1165 40 : }
1166 :
1167 16 : EDS::~EDS()
1168 : {
1169 8 : out_restart_.close();
1170 24 : }
1171 :
1172 0 : void EDS::turnOnDerivatives()
1173 : {
1174 : // do nothing
1175 : // this is to avoid errors triggered when a bias is used as a CV
1176 : // (This is done in ExtendedLagrangian.cpp)
1177 0 : }
1178 :
1179 : }
1180 : } // close the 2 namespaces
|
