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