Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2020 of Glen Hocky
3 :
4 : The FISST 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 FISST 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 : #include "legendre_rule_fast.h"
25 :
26 : #include <iostream>
27 :
28 :
29 : using namespace PLMD;
30 : using namespace bias;
31 :
32 : //namespace is lowercase to match
33 : //module names being all lowercase
34 :
35 : namespace PLMD {
36 : namespace fisst {
37 :
38 : //+PLUMEDOC FISSTMOD_BIAS FISST
39 : /*
40 : Compute and apply the optimal linear force on an observable to enhance sampling of conformational distributions over a range of applied forces.
41 :
42 : This method is described in \cite Hartmann-FISST-2019
43 :
44 : If the system's Hamiltonian is given by:
45 : \f[
46 : H(\vec{p},\vec{q}) = \sum_{j} \frac{p_j^2}{2m_j} + U(\vec{q}),
47 : \f]
48 :
49 : This bias modifies the Hamiltonian to be:
50 : \f[
51 : H'(\vec{p},\vec{q}) = H(\vec{p},\vec{q}) - \bar{F} Q
52 : \f]
53 :
54 : where for CV \f$Q\f$, a coupling constant \f${\bar{F}}\f$ is determined
55 : adaptively according to the FISST algorithm.
56 :
57 : Specifically,
58 : \f[
59 : \bar{F}(Q)=\frac{ \int_{F_{min}}^{F_{max}} e^{\beta F Q(\vec{q})} \omega(F) F dF}{\int_{F_{min}}^{F_{max}} e^{\beta F Q(\vec{q})} \omega(F) dF},
60 : \f]
61 :
62 : where \f$\vec{q}\f$ are the molecular coordinates of the system, and \f$w(F)\f$ is a weighting function that is learned on the fly for each force by the FISST algorithm (starting from an initial weight distribution, uniform by default).
63 :
64 : The target for \f$w(F)=1/Z_q(F)\f$, where
65 : \f[
66 : Z_q(F) \equiv \int d\vec{q} e^{-\beta U(\vec{q}) + \beta F Q(\vec{q})}.
67 : \f]
68 :
69 : FISST also computes and writes Observable Weights \f$W_F(\vec{q}_t)\f$ for a molecular configuration at time \f$t\f$, so that averages of other quantities \f$A(\vec{q})\f$ can be reconstructed later at different force values (over a trajectory with \f$T\f$ samples):
70 : \f[
71 : \langle A \rangle_F = \frac{1}{T} \sum_t W_F(\vec{q}_t) A(\vec{q}_t).
72 : \f]
73 :
74 :
75 : \par Examples
76 :
77 : In the following example, an adaptive restraint is learned to bias the distance between two atoms in a system, for a force range of 0-100 pN.
78 :
79 : \plumedfile
80 : UNITS LENGTH=A TIME=fs ENERGY=kcal/mol
81 :
82 : b1: GROUP ATOMS=1
83 : b2: GROUP ATOMS=12
84 :
85 : dend: DISTANCE ATOMS=b1,b2
86 :
87 : #The conversion factor is 69.4786 pN = 1 kcal/mol/Angstrom
88 :
89 : #0 pN to 100 pN
90 : f: FISST MIN_FORCE=0 MAX_FORCE=1.44 PERIOD=100 NINTERPOLATE=31 ARG=dend OUT_RESTART=pull.restart.txt OUT_OBSERVABLE=pull.observable.txt OBSERVABLE_FREQ=1000
91 :
92 : PRINT ARG=dend,f.dend_fbar,f.bias,f.force2 FILE=pull.colvar.txt STRIDE=1000
93 : \endplumedfile
94 :
95 :
96 : */
97 : //+ENDPLUMEDOC
98 :
99 :
100 : class FISST : public Bias {
101 :
102 :
103 : private:
104 : /*We will get this and store it once, since on-the-fly changing number of CVs will be fatal*/
105 : const unsigned int ncvs_;
106 : std::vector<double> center_;
107 : std::vector<double> current_avg_force_;
108 :
109 : std::vector<double> forces_;
110 : std::vector<double> force_weight_;
111 : std::vector<double> gauss_weight_;
112 : std::vector<double> partition_estimate_;
113 : std::vector<double> observable_weight_;
114 :
115 : std::string in_restart_name_;
116 : std::string out_restart_name_;
117 : std::string out_observable_name_;
118 : std::string fmt_;
119 : std::string initial_weight_dist_;
120 : OFile out_restart_;
121 : OFile out_observable_;
122 : IFile in_restart_;
123 : bool b_freeze_;
124 : bool b_adaptive_;
125 : bool b_restart_;
126 : bool b_write_restart_;
127 : bool b_write_observable_;
128 : bool b_first_restart_sample_;
129 : int period_;
130 : int reset_period_;
131 : int observable_freq_;
132 : int n_interpolation_;
133 : int n_samples_;
134 : double kbt_;
135 : double beta_;
136 : //change min_force and max_force to vectors if going to do more than one cv
137 : double max_force_;
138 : double min_force_;
139 : double initial_weight_rate_;
140 : double threshold_;
141 : Random rand_;
142 :
143 :
144 : Value* value_force2_;
145 : void readInRestart();
146 : void NormalizeForceWeights();
147 : /*setup output restart*/
148 : void setupOutRestart();
149 : void setupOutObservable();
150 : /*write output restart*/
151 : void writeOutRestart();
152 : void writeOutObservable();
153 : void update_statistics();
154 : void update_bias();
155 : void apply_bias();
156 : void compute_observable_weight();
157 :
158 : public:
159 : explicit FISST(const ActionOptions&);
160 : void calculate();
161 : void update();
162 : void turnOnDerivatives();
163 : static void registerKeywords(Keywords& keys);
164 : ~FISST();
165 : };
166 :
167 13789 : PLUMED_REGISTER_ACTION(FISST,"FISST")
168 :
169 6 : void FISST::registerKeywords(Keywords& keys) {
170 6 : Bias::registerKeywords(keys);
171 6 : keys.use("ARG");
172 12 : keys.add("compulsory","PERIOD","Steps corresponding to the learning rate");
173 12 : keys.add("optional","RESET_PERIOD","Reset the learning statistics every time this number of steps comes around.");
174 12 : keys.add("compulsory","NINTERPOLATE","Number of grid points on which to do interpolation.");
175 12 : keys.add("compulsory","MIN_FORCE","Minimum force (per CV) to use for sampling. Units: [Energy]/[CV] (can be negative).");
176 12 : keys.add("compulsory","MAX_FORCE","Maximum force (per CV) to use for sampling.");
177 12 : keys.add("compulsory","CENTER","0","The CV value at which the applied bias energy will be zero");
178 12 : keys.add("optional","KBT","The system temperature in units of KB*T. If not provided will be taken from MD code (if available)");
179 :
180 12 : keys.add("optional","INITIAL_WEIGHT_DIST","Starting distribution for the force weights (options: UNIFORM, EXP, GAUSS).");
181 12 : keys.add("optional","INITIAL_WEIGHT_RATE","Rate of decay for exponential and gaussian distributions. W(F)~exp(-r |F|^d).");
182 :
183 12 : keys.add("optional","RESTART_FMT","the format that should be used to output real numbers in FISST restarts.");
184 12 : keys.add("optional","OUT_RESTART","Output file for all information needed to continue FISST simulation."
185 : "If you have the RESTART directive set (global or for FISST), this file will be appended to."
186 : "Note that the header will be printed again if appending.");
187 12 : keys.add("optional","IN_RESTART","Read this file to continue an FISST simulation. "
188 : "If same as OUT_RESTART and you have not set the RESTART directive, the file will be backed-up and overwritten with new output."
189 : "If you do have the RESTART flag set and it is the same name as OUT_RESTART, this file will be appended.");
190 12 : keys.add("optional","OUT_OBSERVABLE","Output file putting weights needed to compute observables at different force values."
191 : "If you have the RESTART directive set (global or for FISST), this file will be appended to. "
192 : "Note that the header will be printed again if appending.");
193 12 : keys.add("optional","OBSERVABLE_FREQ","How often to write out observable weights (default=period).");
194 12 : keys.addFlag("FREEZE",false,"Fix bias weights at current level (only used for restarting).");
195 6 : keys.use("RESTART");
196 12 : keys.addOutputComponent("force2","default","squared value of force from the bias.");
197 12 : keys.addOutputComponent("_fbar","default", "For each named CV biased, there will be a corresponding output CV_fbar storing the current linear bias prefactor.");
198 6 : }
199 :
200 2 : FISST::FISST(const ActionOptions&ao):
201 : PLUMED_BIAS_INIT(ao),
202 2 : ncvs_(getNumberOfArguments()),
203 0 : current_avg_force_(ncvs_,0.0),
204 2 : center_(ncvs_,0.0),
205 : //change min_force and max_force to vectors if going to do more than one cv
206 2 : min_force_(0.0),
207 2 : max_force_(0.0),
208 2 : in_restart_name_(""),
209 2 : out_restart_name_(""),
210 2 : out_observable_name_(""),
211 2 : fmt_("%e"),
212 2 : b_freeze_(false),
213 2 : b_restart_(false),
214 2 : b_write_restart_(false),
215 2 : b_write_observable_(false),
216 2 : b_first_restart_sample_(true),
217 2 : n_interpolation_(0),
218 2 : n_samples_(0),
219 2 : initial_weight_rate_(0),
220 2 : initial_weight_dist_("UNIFORM"),
221 2 : period_(0),
222 2 : reset_period_(0),
223 2 : observable_freq_(0),
224 2 : kbt_(0.0),
225 6 : value_force2_(NULL) {
226 2 : if(ncvs_==0) {
227 0 : error("Must specify at least one CV with ARG");
228 : }
229 :
230 : //temporary
231 2 : if(ncvs_>1) {
232 0 : error("FISST only supports using one CV right now");
233 : }
234 :
235 2 : addComponent("force2");
236 2 : componentIsNotPeriodic("force2");
237 2 : value_force2_ = getPntrToComponent("force2");
238 :
239 4 : for(unsigned int i = 0; i<ncvs_; i++) {
240 2 : std::string comp = getPntrToArgument(i)->getName() + "_fbar";
241 2 : addComponent(comp);
242 2 : componentIsNotPeriodic(comp);
243 : }
244 :
245 2 : parseVector("CENTER",center_);
246 : //change min_force and max_force to vectors if going to do more than one cv
247 2 : parse("MIN_FORCE",min_force_);
248 2 : parse("MAX_FORCE",max_force_);
249 2 : parse("PERIOD",period_);
250 2 : parse("RESET_PERIOD",reset_period_);
251 2 : parse("INITIAL_WEIGHT_DIST",initial_weight_dist_);
252 2 : parse("INITIAL_WEIGHT_RATE",initial_weight_rate_);
253 2 : parse("OBSERVABLE_FREQ",observable_freq_);
254 2 : parse("NINTERPOLATE",n_interpolation_);
255 2 : parseFlag("FREEZE",b_freeze_);
256 2 : parse("KBT",kbt_);
257 2 : parse("RESTART_FMT", fmt_);
258 2 : fmt_ = " " + fmt_;//add space since parse strips them
259 2 : parse("OUT_RESTART",out_restart_name_);
260 2 : parse("OUT_OBSERVABLE",out_observable_name_);
261 2 : parse("IN_RESTART",in_restart_name_);
262 2 : checkRead();
263 :
264 2 : if(center_.size() != ncvs_) {
265 0 : error("Must have same number of CENTER arguments as ARG arguments");
266 : }
267 :
268 2 : if(in_restart_name_ != "") {
269 0 : b_restart_ = true;
270 0 : log.printf(" reading simulation information from file: %s\n",in_restart_name_.c_str());
271 0 : readInRestart();
272 : } else {
273 :
274 2 : if(! kbt_ > 0.0) {
275 2 : kbt_ = plumed.getAtoms().getKbT();
276 : }
277 :
278 : //in driver, this results in kbt of 0
279 2 : if(kbt_ == 0) {
280 0 : error(" Unable to determine valid kBT. "
281 : "Could be because you are runnning from driver or MD didn't give temperature.\n"
282 : "Consider setting temperature manually with the KBT keyword.");
283 : }
284 :
285 2 : log.printf(" kBT = %f\n",kbt_);
286 2 : log.printf(" Updating with a time scale of %i steps\n",period_);
287 :
288 2 : log.printf(" Using centers for CVs of:");
289 4 : for(unsigned int i = 0; i< ncvs_; i++) {
290 2 : log.printf(" %f ",center_[i]);
291 : }
292 2 : log.printf("\n");
293 2 : observable_weight_.resize(n_interpolation_);
294 64 : for(unsigned int i = 0; i<n_interpolation_; i++) {
295 62 : observable_weight_[i] = 1.0;
296 : }
297 :
298 2 : forces_.resize(n_interpolation_);
299 2 : force_weight_.resize(n_interpolation_);
300 : //using code from the MIST project
301 2 : gauss_weight_.resize(n_interpolation_);
302 2 : legendre_compute_glr(n_interpolation_, &forces_[0], &gauss_weight_[0]);
303 2 : rescale(min_force_, max_force_, n_interpolation_, &forces_[0], &gauss_weight_[0]);
304 :
305 2 : log.printf("Using weight distribution %s with rate %f\n",initial_weight_dist_.c_str(),initial_weight_rate_);
306 2 : if(initial_weight_dist_ == "UNIFORM" ) {
307 32 : for(unsigned int i = 0; i<n_interpolation_; i++) {
308 31 : force_weight_[i] = 1.0;
309 : }
310 1 : } else if (initial_weight_dist_ == "EXP" ) {
311 32 : for(unsigned int i = 0; i<n_interpolation_; i++) {
312 31 : force_weight_[i] = exp(-fabs(forces_[i])*initial_weight_rate_);
313 : }
314 0 : } else if (initial_weight_dist_ == "GAUSS" ) {
315 0 : for(unsigned int i = 0; i<n_interpolation_; i++) {
316 0 : force_weight_[i] = exp(-pow(forces_[i],2)*initial_weight_rate_);
317 : }
318 : } else {
319 0 : error(" Specified weight distribution is not from the allowed list.");
320 :
321 : }
322 :
323 2 : partition_estimate_.resize(n_interpolation_);
324 2 : NormalizeForceWeights();
325 : double sum = 0.0;
326 64 : for(unsigned int i = 0; i<n_interpolation_; i++) {
327 : //setting partition estimate as 1/w_i
328 62 : partition_estimate_[i] = 1/force_weight_[i];
329 62 : log.printf("force/gauss weight/force_weight: %i %f %f %f\n",i,forces_[i],gauss_weight_[i],force_weight_[i]);
330 62 : sum+=gauss_weight_[i]*force_weight_[i];
331 : }
332 2 : log.printf("--Sum_i w_i g_i: %f\n",sum);
333 :
334 : }
335 :
336 : //set inverse temperature
337 2 : beta_ = 1/kbt_;
338 :
339 2 : if(b_freeze_ && b_restart_) {
340 0 : log.printf(" freezing weights read in from the restart file\n");
341 : }
342 :
343 2 : if(out_restart_name_.length()>0) {
344 2 : log.printf(" writing restart information every %i steps to file %s with format %s\n",abs(period_),out_restart_name_.c_str(), fmt_.c_str());
345 2 : b_write_restart_ = true;
346 2 : setupOutRestart();
347 : }
348 2 : if(out_observable_name_.length()>0) {
349 2 : if(observable_freq_==0) {
350 2 : observable_freq_ = period_;
351 : }
352 2 : log.printf(" writing observable information every %i steps to file %s with format %s\n",observable_freq_,out_observable_name_.c_str(), fmt_.c_str());
353 2 : b_write_observable_ = true;
354 2 : setupOutObservable();
355 : }
356 :
357 : //add citation later:
358 : //log<<" Bibliography "<<plumed.cite("")<<"\n";
359 2 : }
360 :
361 12 : void FISST::NormalizeForceWeights() {
362 : double denom = 0.0;
363 :
364 384 : for(unsigned i=0; i<n_interpolation_; i++) {
365 372 : denom += gauss_weight_[i] * force_weight_[i];
366 : }
367 :
368 384 : for(unsigned i=0; i<n_interpolation_; i++) {
369 372 : force_weight_[i] /= denom;
370 : }
371 12 : }
372 :
373 0 : void FISST::readInRestart() {
374 0 : in_restart_.open(in_restart_name_);
375 :
376 0 : if(in_restart_.FieldExist("kbt")) {
377 0 : in_restart_.scanField("kbt",kbt_);
378 : } else {
379 0 : error("No field 'kbt' in restart file");
380 : }
381 0 : log.printf(" with kBT = %f\n",kbt_);
382 :
383 0 : if(in_restart_.FieldExist("period")) {
384 0 : in_restart_.scanField("period",period_);
385 : } else {
386 0 : error("No field 'period' in restart file");
387 : }
388 0 : log.printf(" Updating every %i steps\n",period_);
389 :
390 : //this one can be optional
391 0 : if(in_restart_.FieldExist("reset_period")) {
392 0 : in_restart_.scanField("reset_period",reset_period_);
393 : }
394 0 : log.printf(" Resetting statistics every %i steps\n",reset_period_);
395 :
396 0 : if(in_restart_.FieldExist("n_interpolation")) {
397 0 : in_restart_.scanField("n_interpolation",n_interpolation_);
398 : } else {
399 0 : error("No field 'n_interpolation' in restart file");
400 : }
401 :
402 0 : if(in_restart_.FieldExist("min_force")) {
403 0 : in_restart_.scanField("min_force",min_force_);
404 : } else {
405 0 : error("No field 'min_force' in restart file");
406 : }
407 0 : if(in_restart_.FieldExist("max_force")) {
408 0 : in_restart_.scanField("max_force",max_force_);
409 : } else {
410 0 : error("No field 'max_force' in restart file");
411 : }
412 0 : log.printf(" with forces from min_force=%e to max_force=%e over %i bins\n",min_force_,max_force_,n_interpolation_);
413 :
414 : unsigned int N = 0;
415 : std::string cv_name;
416 : double tmp, time;
417 :
418 0 : while(in_restart_.scanField("time",time)) {
419 0 : in_restart_.scanField("nsamples",n_samples_);
420 :
421 0 : observable_weight_.resize(n_interpolation_);
422 0 : partition_estimate_.resize(n_interpolation_);
423 0 : force_weight_.resize(n_interpolation_);
424 0 : gauss_weight_.resize(n_interpolation_);
425 0 : forces_.resize(n_interpolation_);
426 :
427 0 : for(unsigned int i = 0; i<ncvs_; ++i) {
428 : cv_name = getPntrToArgument(i)->getName();
429 0 : in_restart_.scanField(cv_name,tmp);
430 0 : for(unsigned int j =0; j<n_interpolation_; ++j) {
431 0 : in_restart_.scanField(cv_name + "_f"+std::to_string(j),forces_[j]);
432 0 : in_restart_.scanField(cv_name + "_g"+std::to_string(j),gauss_weight_[j]);
433 0 : in_restart_.scanField(cv_name + "_w"+std::to_string(j),force_weight_[j]);
434 0 : in_restart_.scanField(cv_name + "_z"+std::to_string(j),partition_estimate_[j]);
435 : }
436 : }
437 : N++;
438 :
439 0 : in_restart_.scanField();
440 : }
441 :
442 : double sum = 0.0;
443 0 : for(unsigned int j =0; j<n_interpolation_; ++j) {
444 : //clear observable weight, which will be set later
445 0 : observable_weight_[j] = 1.0;
446 :
447 : //setting partition estimate as 1/w_i
448 0 : log.printf("force/gauss weight/force_weight: %i %e %e %e\n",j,forces_[j],gauss_weight_[j],force_weight_[j]);
449 0 : sum+=gauss_weight_[j]*force_weight_[j];
450 : }
451 0 : log.printf("--Sum_i w_i g_i: %f\n",sum);
452 :
453 0 : in_restart_.close();
454 0 : }
455 :
456 2 : void FISST::setupOutObservable() {
457 2 : out_observable_.link(*this);
458 2 : out_observable_.fmtField(fmt_);
459 2 : out_observable_.open(out_observable_name_);
460 : out_observable_.setHeavyFlush();
461 :
462 4 : out_observable_.addConstantField("kbt").printField("kbt",kbt_);
463 4 : out_observable_.addConstantField("n_interpolation").printField("n_interpolation",n_interpolation_);
464 4 : out_observable_.addConstantField("period").printField("period",period_);
465 4 : out_observable_.addConstantField("min_force").printField("min_force",min_force_);
466 4 : out_observable_.addConstantField("max_force").printField("max_force",max_force_);
467 2 : }
468 :
469 2 : void FISST::setupOutRestart() {
470 2 : out_restart_.link(*this);
471 2 : out_restart_.fmtField(fmt_);
472 2 : out_restart_.open(out_restart_name_);
473 : out_restart_.setHeavyFlush();
474 :
475 4 : out_restart_.addConstantField("kbt").printField("kbt",kbt_);
476 4 : out_restart_.addConstantField("n_interpolation").printField("n_interpolation",n_interpolation_);
477 4 : out_restart_.addConstantField("period").printField("period",period_);
478 2 : if(reset_period_>0) {
479 0 : out_restart_.addConstantField("reset_period").printField("reset_period",reset_period_);
480 : }
481 4 : out_restart_.addConstantField("min_force").printField("min_force",min_force_);
482 4 : out_restart_.addConstantField("max_force").printField("max_force",max_force_);
483 2 : }
484 :
485 10 : void FISST::writeOutRestart() {
486 : std::string cv_name;
487 10 : out_restart_.printField("time",getTimeStep()*getStep());
488 10 : out_restart_.printField("nsamples",n_samples_);
489 :
490 20 : for(unsigned int i = 0; i<ncvs_; ++i) {
491 : cv_name = getPntrToArgument(i)->getName();
492 10 : double Q_i = difference(i, center_[i], getArgument(i));
493 10 : out_restart_.printField(cv_name,Q_i);
494 320 : for(int j = 0; j < n_interpolation_; j++ ) {
495 : //have to update this for multiple cvs
496 620 : out_restart_.printField(cv_name + "_f"+std::to_string(j),forces_[j]);
497 620 : out_restart_.printField(cv_name + "_g"+std::to_string(j),gauss_weight_[j]);
498 620 : out_restart_.printField(cv_name + "_w"+std::to_string(j),force_weight_[j]);
499 620 : out_restart_.printField(cv_name + "_z"+std::to_string(j),partition_estimate_[j]);
500 : }
501 : }
502 10 : out_restart_.printField();
503 10 : }
504 :
505 10 : void FISST::writeOutObservable() {
506 : std::string cv_name;
507 10 : out_observable_.printField("time",getTimeStep()*getStep());
508 10 : out_observable_.printField("nsamples",n_samples_);
509 :
510 20 : for(unsigned int i = 0; i<ncvs_; ++i) {
511 : cv_name = getPntrToArgument(i)->getName();
512 10 : double Q_i = difference(i, center_[i], getArgument(i));
513 10 : out_observable_.printField(cv_name,Q_i);
514 10 : out_observable_.printField(cv_name + "_fbar",current_avg_force_[i]);
515 320 : for(int j = 0; j < n_interpolation_; j++ ) {
516 : //have to update this for multiple cvs
517 620 : out_observable_.printField(cv_name + "_f"+std::to_string(j),forces_[j]);
518 620 : out_observable_.printField(cv_name + "_ow"+std::to_string(j),observable_weight_[j]);
519 : }
520 : }
521 10 : out_observable_.printField();
522 10 : }
523 :
524 :
525 10 : void FISST::calculate() {
526 10 : if(getStep() == 0 ) {
527 2 : if(b_write_restart_) {
528 2 : writeOutRestart();
529 : }
530 2 : if(b_write_observable_) {
531 2 : writeOutObservable();
532 : }
533 : }
534 :
535 10 : if(! b_freeze_) {
536 10 : if(b_restart_ && b_first_restart_sample_) {
537 : //dont' update statistics if restarting and first sample
538 0 : b_first_restart_sample_ = false;
539 : } else {
540 10 : update_statistics();
541 : }
542 : }
543 10 : update_bias();
544 10 : apply_bias();
545 :
546 : //check about writing restart file
547 10 : if(getStep()>0 && getStep()%period_==0) {
548 8 : if(b_write_restart_) {
549 8 : writeOutRestart();
550 : }
551 : }
552 10 : if(getStep()>0 && getStep()%observable_freq_==0) {
553 8 : if(b_write_observable_) {
554 8 : compute_observable_weight();
555 8 : writeOutObservable();
556 : }
557 : }
558 10 : }
559 :
560 :
561 10 : void FISST::apply_bias() {
562 : //Compute linear force as in "restraint"
563 : double ene = 0, totf2 = 0, cv, m, f;
564 :
565 20 : for(unsigned int i = 0; i < ncvs_; ++i) {
566 10 : cv = difference(i, center_[i], getArgument(i));
567 10 : double fbar = current_avg_force_[i];
568 10 : ene -= fbar*cv;
569 : setOutputForce(i,fbar);
570 10 : totf2 += fbar*fbar;
571 :
572 10 : std::string fbar_name_ = getPntrToArgument(i)->getName() + "_fbar";
573 10 : Value* fbar_ = getPntrToComponent(fbar_name_);
574 : fbar_->set(fbar);
575 : };
576 :
577 : setBias(ene);
578 10 : value_force2_->set(totf2);
579 : //log.flush();
580 10 : }
581 :
582 10 : void FISST::update_statistics() {
583 : //get stride is for multiple time stepping
584 10 : double dt=getTimeStep()*getStride();
585 10 : double h = dt/(period_*getTimeStep());
586 : double fbar_denum_integral = 0.0;
587 :
588 10 : int step = getStep();
589 10 : if(reset_period_>0 && step>0 && step%reset_period_==0) {
590 0 : n_samples_=1;
591 : } else {
592 10 : n_samples_++;
593 : }
594 10 : double d_n_samples = (double)n_samples_;
595 :
596 20 : for(unsigned int i = 0; i < ncvs_; ++i) {
597 10 : double Q_i = difference(i, center_[i], getArgument(i));
598 320 : for(unsigned int j=0; j<n_interpolation_; j++) {
599 : //if multiple cvs, these need to be updated to have 2 columns
600 310 : double f_j = forces_[j];
601 310 : double w_j = force_weight_[j];
602 310 : double g_j = gauss_weight_[j];
603 :
604 310 : fbar_denum_integral += g_j * w_j * exp(beta_*f_j * Q_i);
605 : }
606 :
607 320 : for(unsigned int j=0; j<n_interpolation_; j++) {
608 310 : double f_j = forces_[j];
609 310 : double sample_weight = exp(beta_*f_j * Q_i) / fbar_denum_integral;
610 :
611 310 : partition_estimate_[j] = sample_weight/d_n_samples + partition_estimate_[j]*(d_n_samples-1)/(d_n_samples);
612 :
613 310 : double w_jn = force_weight_[j];
614 310 : double z_jn = partition_estimate_[j];
615 :
616 310 : double w_jp1 = (1.0 - h) * w_jn + h / z_jn;
617 310 : force_weight_[j] = w_jp1;
618 : }
619 : }
620 :
621 : // make sure that the weights are normalised
622 10 : NormalizeForceWeights();
623 10 : }
624 :
625 :
626 10 : void FISST::update_bias() {
627 20 : for(unsigned int i = 0; i < ncvs_; ++i) {
628 10 : double Q_i = difference(i, center_[i], getArgument(i));
629 : double fbar_num_integral = 0.0;
630 : double fbar_denum_integral = 0.0;
631 :
632 320 : for(unsigned int j=0; j<n_interpolation_; j++ ) {
633 310 : double f_j = forces_[j];
634 310 : double w_j = force_weight_[j];
635 310 : double g_j = gauss_weight_[j];
636 :
637 310 : fbar_num_integral += g_j * f_j * w_j * exp(beta_*f_j*Q_i);
638 310 : fbar_denum_integral += g_j * w_j * exp(beta_*f_j*Q_i);
639 : }
640 :
641 10 : current_avg_force_[i] = fbar_num_integral/fbar_denum_integral;
642 : }
643 10 : }
644 :
645 8 : void FISST::compute_observable_weight() {
646 8 : double obs_num = (max_force_ - min_force_);
647 :
648 16 : for(unsigned int i = 0; i < ncvs_; ++i) {
649 8 : double Q_i = difference(i, center_[i], getArgument(i));
650 :
651 256 : for(unsigned int j=0; j<n_interpolation_; j++ ) {
652 248 : double z_j = partition_estimate_[j];
653 248 : double f_j = forces_[j];
654 : double denum_integral = 0.0;
655 :
656 7936 : for( unsigned int k=0; k<n_interpolation_; k++ ) {
657 7688 : double f_k = forces_[k];
658 7688 : double w_k = force_weight_[k];
659 7688 : double g_k = gauss_weight_[k];
660 :
661 7688 : denum_integral += g_k * w_k * exp(beta_*(f_k-f_j)*Q_i);
662 : }
663 248 : observable_weight_[j] = obs_num/(denum_integral*z_j);
664 : }
665 : }
666 8 : }
667 :
668 :
669 :
670 10 : void FISST::update() {
671 : //pass
672 10 : }
673 :
674 4 : FISST::~FISST() {
675 2 : out_restart_.close();
676 2 : out_observable_.close();
677 6 : }
678 :
679 0 : void FISST::turnOnDerivatives() {
680 : // do nothing
681 : // this is to avoid errors triggered when a bias is used as a CV
682 : // (This is done in ExtendedLagrangian.cpp)
683 0 : }
684 :
685 :
686 : }
687 : }//close the 2 namespaces
|
