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