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