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