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