Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017 of Haochuan Chen (excluding colvar_UIestimator.h)
3 : Copyright (c) 2017 of Haohao Fu (colvar_UIestimator.h)
4 :
5 : This program is free software: you can redistribute it and/or modify
6 : it under the terms of the GNU Lesser General Public License as published
7 : by the Free Software Foundation, either version 3 of the License, or
8 : (at your option) any later version.
9 :
10 : This program is distributed in the hope that it will be useful,
11 : but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 : GNU Lesser General Public License for more details.
14 :
15 : You should have received a copy of the GNU Lesser General Public License
16 : along with this program. If not, see <http://www.gnu.org/licenses/>.
17 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
18 : #ifdef __PLUMED_HAS_BOOST_SERIALIZATION
19 : #include "core/ActionRegister.h"
20 : #include "bias/Bias.h"
21 : #include "core/Atoms.h"
22 : #include "core/PlumedMain.h"
23 : #include "DRR.h"
24 : #include "tools/Random.h"
25 : #include "tools/Tools.h"
26 : #include "colvar_UIestimator.h"
27 :
28 : #include <boost/archive/binary_iarchive.hpp>
29 : #include <boost/archive/binary_oarchive.hpp>
30 : #include <boost/serialization/vector.hpp>
31 : #include <cmath>
32 : #include <fstream>
33 : #include <iomanip>
34 : #include <iostream>
35 : #include <limits>
36 : #include <random>
37 : #include <string>
38 :
39 : using namespace PLMD;
40 : using namespace bias;
41 : using namespace std;
42 :
43 : namespace PLMD {
44 : namespace drr {
45 :
46 : //+PLUMEDOC EABFMOD_BIAS DRR
47 : /*
48 : Used to performed extended-system adaptive biasing force(eABF) \cite Lelievre2007 method
49 : on one or more collective variables. This method is also
50 : called dynamic reference restraining(DRR) \cite Zheng2012 .
51 :
52 : For each collective variable \f$\xi_i\f$, a fictitious variable \f$\lambda_i\f$
53 : is attached through a spring. The fictitious variable \f$\lambda_i\f$ undergoes
54 : overdamped langevin dynamics jusk like \ref EXTENDED_LAGRANGIAN. The ABF
55 : algorithm applies bias force on \f$\lambda_i\f$. The bias force acts on
56 : \f$\lambda_i\f$ is the negative average spring force on \f$\lambda_i\f$, which
57 : enhances the sampling of \f$\lambda_i\f$.
58 :
59 : \f[
60 : F_{bias}(\lambda_i)=k(\lambda_i-\langle\xi_i\rangle_{\lambda_i})
61 : \f]
62 :
63 : If spring force constant k is large enough, then \f$\xi_i\f$ synchronizes with
64 : \f$\lambda_i\f$. The naive(ABF) estimator is just the negative
65 : average spring force of \f$\lambda_i\f$.
66 :
67 : The naive(ABF) estimator is biased. There are unbiased estimators such as
68 : CZAR(Corrected z-averaged restraint) \cite Lesage2016 and UI(Umbrella
69 : Integration).
70 : The CZAR estimates the gradients as:
71 :
72 : \f[
73 : \frac{\partial{A}}{\partial{\xi_i}}\left({\xi}\right)=-\frac{1}{\beta}\frac{\partial\ln\tilde{\rho}\left(\xi\right)}{\partial{\xi_i}}+k\left(\langle\lambda_i\rangle_\xi-\xi_i\right)
74 : \f]
75 :
76 : The UI estimates the gradients as:
77 : \f[
78 : A'(\xi^*)=\frac{{\sum_\lambda}N\left(\xi^*,\lambda\right)\left[\frac{\xi^*-\langle\xi\rangle_\lambda}{\beta\sigma_\lambda^2}-k(\xi^*-\lambda)\right]}{{\sum_\lambda}N\left(\xi^*,\lambda\right)}
79 : \f]
80 :
81 : The code performing UI(colvar_UIestimator.h) is contributed by Haohao Fu \cite Fu2016 .
82 : It may be slow. I only change the boltzmann constant and output
83 : precision in it. For new version and issues, please see:
84 : https://github.com/fhh2626/colvars
85 :
86 : After running eABF/DRR, the \ref drr_tool utility can be used to extract the gradients and counts files from .drrstate. Naive(ABF) estimator's result is in .abf.grad and .abf.count files and CZAR estimator's result is in .czar.grad and .czar.count files. To get PMF, the abf_integrate(https://github.com/Colvars/colvars/tree/master/colvartools) is useful.
87 :
88 : \par Examples
89 :
90 : The following input tells plumed to perform a eABF/DRR simulation on two
91 : torsional angles.
92 : \plumedfile
93 : phi: TORSION ATOMS=5,7,9,15
94 : psi: TORSION ATOMS=7,9,15,17
95 :
96 : DRR ...
97 : LABEL=eabf
98 : ARG=phi,psi
99 : FULLSAMPLES=500
100 : GRID_MIN=-pi,-pi
101 : GRID_MAX=pi,pi
102 : GRID_BIN=180,180
103 : FRICTION=8.0,8.0
104 : TAU=0.5,0.5
105 : OUTPUTFREQ=50000
106 : HISTORYFREQ=500000
107 : ... DRR
108 :
109 : # monitor the two variables, their fictitious variables and applied forces.
110 : PRINT STRIDE=10 ARG=phi,psi,eabf.phi_fict,eabf.psi_fict,eabf.phi_biasforce,eabf.psi_biasforce FILE=COLVAR
111 : \endplumedfile
112 :
113 : The following input tells plumed to perform a eABF/DRR simulation on the
114 : distance of atom 10 and 92. The distance is restraint by \ref LOWER_WALLS and
115 : \ref UPPER_WALLS.
116 : \plumedfile
117 : dist1: DISTANCE ATOMS=10,92
118 : eabf_winall: DRR ARG=dist1 FULLSAMPLES=2000 GRID_MIN=1.20 GRID_MAX=3.20 GRID_BIN=200 FRICTION=8.0 TAU=0.5 OUTPUTFREQ=5000 HISTORYFREQ=500000
119 : uwall: UPPER_WALLS ARG=eabf_winall.dist1_fict AT=3.2 KAPPA=418.4
120 : lwall: LOWER_WALLS ARG=eabf_winall.dist1_fict AT=1.2 KAPPA=418.4
121 : PRINT STRIDE=10 ARG=dist1,eabf_winall.dist1_fict,eabf_winall.dist1_biasforce FILE=COLVAR
122 : \endplumedfile
123 :
124 : */
125 : //+ENDPLUMEDOC
126 :
127 : using std::vector;
128 : using std::string;
129 :
130 8 : class DynamicReferenceRestraining : public Bias {
131 : private:
132 : bool firsttime;
133 : bool nobias;
134 : vector<double> fictNoPBC;
135 : vector<double> real;
136 : vector<double> springlength; // spring lengths
137 : vector<double> fict; // coordinates of extended variables
138 : vector<double> vfict; // velocities of extended variables
139 : vector<double> vfict_laststep;
140 : vector<double> ffict; // forces exerted on extended variables
141 : vector<double> fbias; // bias forces from eABF
142 : vector<double> kappa;
143 : vector<double> tau;
144 : vector<double> friction;
145 : vector<double> etemp;
146 : vector<double> ffict_measured;
147 : vector<Value *> biasforceValue;
148 : vector<Value *> fictValue;
149 : vector<Value *> vfictValue;
150 : vector<double> c1;
151 : vector<double> c2;
152 : vector<double> mass;
153 : vector<DRRAxis> delim;
154 : string outputname;
155 : string cptname;
156 : string outputprefix;
157 : const size_t ndims;
158 : double dt;
159 : double kbt;
160 : double outputfreq;
161 : double historyfreq;
162 : bool isRestart;
163 : bool useCZARestimator;
164 : bool useUIestimator;
165 : bool textoutput;
166 : ABF ABFGrid;
167 : CZAR CZARestimator;
168 : double fullsamples;
169 : UIestimator::UIestimator eabf_UI;
170 : Random rand;
171 :
172 : public:
173 : explicit DynamicReferenceRestraining(const ActionOptions &);
174 : void calculate();
175 : void update();
176 : void save(const string &filename, long long int step);
177 : void load(const string &filename);
178 : void backupFile(const string &filename);
179 : static void registerKeywords(Keywords &keys);
180 : bool is_file_exist(const char *fileName);
181 : };
182 :
183 4825 : PLUMED_REGISTER_ACTION(DynamicReferenceRestraining, "DRR")
184 :
185 5 : void DynamicReferenceRestraining::registerKeywords(Keywords &keys) {
186 5 : Bias::registerKeywords(keys);
187 5 : keys.use("ARG");
188 : keys.add("optional", "KAPPA", "specifies that the restraint is harmonic and "
189 : "what the values of the force constants on "
190 : "each of the variables are (default to "
191 5 : "kbt/(GRID_SPACING)^2)");
192 : keys.add("compulsory", "TAU", "0.5", "specifies relaxation time on each of "
193 : "variables are, similar to "
194 5 : "extendedTimeConstant in Colvars");
195 : keys.add("compulsory", "FRICTION", "8.0",
196 : "add a friction to the variable, similar to extendedLangevinDamping "
197 5 : "in Colvars");
198 : keys.add("compulsory", "GRID_MIN", "the lower bounds for the grid (GRID_BIN "
199 5 : "or GRID_SPACING should be specified)");
200 : keys.add("compulsory", "GRID_MAX", "the upper bounds for the grid (GRID_BIN "
201 5 : "or GRID_SPACING should be specified)");
202 5 : keys.add("optional", "GRID_BIN", "the number of bins for the grid");
203 : keys.add("optional", "GRID_SPACING", "the approximate grid spacing (to be "
204 : "used as an alternative or together "
205 5 : "with GRID_BIN)");
206 : keys.add("compulsory", "FULLSAMPLES", "500",
207 5 : "number of samples in a bin prior to application of the ABF");
208 5 : keys.add("compulsory", "OUTPUTFREQ", "write results to a file every N steps");
209 5 : keys.add("optional", "HISTORYFREQ", "save history to a file every N steps");
210 5 : keys.addFlag("NOCZAR", false, "disable the CZAR estimator");
211 : keys.addFlag("UI", false,
212 5 : "enable the umbrella integration estimator");
213 : keys.add("optional", "UIRESTARTPREFIX",
214 5 : "specify the restart files for umbrella integration");
215 : keys.add("optional", "OUTPUTPREFIX",
216 5 : "specify the output prefix (default to the label name)");
217 : keys.add("optional", "TEMP", "the system temperature - needed when FRICTION "
218 : "is present. If not provided will be taken from "
219 5 : "MD code (if available)");
220 : keys.add(
221 : "optional", "EXTTEMP",
222 5 : "the temperature of extended variables (default to system temperature)");
223 : keys.add("optional", "DRR_RFILE",
224 5 : "specifies the restart file (.drrstate file)");
225 5 : keys.addFlag("NOBIAS", false, "DO NOT apply bias forces.");
226 : keys.addFlag("TEXTOUTPUT", false, "use text output for grad and count files "
227 : "instead of boost::serialization binary "
228 5 : "output");
229 5 : componentsAreNotOptional(keys);
230 : keys.addOutputComponent(
231 : "_fict", "default",
232 : "one or multiple instances of this quantity will be refereceable "
233 : "elsewhere in the input file. "
234 : "These quantities will named with the arguments of the bias followed by "
235 : "the character string _tilde. It is possible to add forces on these "
236 5 : "variable.");
237 : keys.addOutputComponent(
238 : "_vfict", "default",
239 : "one or multiple instances of this quantity will be refereceable "
240 : "elsewhere in the input file. "
241 : "These quantities will named with the arguments of the bias followed by "
242 : "the character string _tilde. It is NOT possible to add forces on these "
243 5 : "variable.");
244 : keys.addOutputComponent(
245 : "_biasforce", "default",
246 5 : "The bias force from eABF/DRR of the fictitious particle.");
247 5 : }
248 :
249 4 : DynamicReferenceRestraining::DynamicReferenceRestraining(
250 : const ActionOptions &ao)
251 : : PLUMED_BIAS_INIT(ao), firsttime(true), nobias(false),
252 8 : fictNoPBC(getNumberOfArguments(), 0.0), real(getNumberOfArguments(), 0.0),
253 4 : springlength(getNumberOfArguments(), 0.0),
254 8 : fict(getNumberOfArguments(), 0.0), vfict(getNumberOfArguments(), 0.0),
255 4 : vfict_laststep(getNumberOfArguments(), 0.0),
256 8 : ffict(getNumberOfArguments(), 0.0), fbias(getNumberOfArguments(), 0.0),
257 8 : kappa(getNumberOfArguments(), 0.0), tau(getNumberOfArguments(), 0.0),
258 8 : friction(getNumberOfArguments(), 0.0), etemp(getNumberOfArguments(), 0.0),
259 4 : ffict_measured(getNumberOfArguments(), 0.0),
260 4 : biasforceValue(getNumberOfArguments(), NULL),
261 4 : fictValue(getNumberOfArguments(), NULL),
262 8 : vfictValue(getNumberOfArguments(), NULL), c1(getNumberOfArguments(), 0.0),
263 8 : c2(getNumberOfArguments(), 0.0), mass(getNumberOfArguments(), 0.0),
264 4 : delim(getNumberOfArguments()), outputname(""), cptname(""),
265 4 : outputprefix(""), ndims(getNumberOfArguments()), dt(0.0), kbt(0.0),
266 : outputfreq(0.0), historyfreq(-1.0), isRestart(false),
267 88 : useCZARestimator(true), useUIestimator(false), textoutput(false)
268 : {
269 : log << "eABF/DRR: You now are using the extended adaptive biasing "
270 4 : "force(eABF) method."
271 8 : << '\n';
272 : log << "eABF/DRR: Some people also refer to it as dynamic reference "
273 4 : "restraining(DRR) method."
274 8 : << '\n';
275 : log << "eABF/DRR: Currently the CZAR and naive(ABF on extended variables) "
276 4 : "estimator is enabled by default."
277 8 : << '\n';
278 : log << "eABF/DRR: For reasons of performance, the umbrella integration "
279 4 : "estimator is not enabled by default."
280 8 : << '\n';
281 : log << "eABF/DRR: This method is originally implemented in "
282 4 : "colvars(https://github.com/colvars/colvars)."
283 8 : << '\n';
284 : log << "eABF/DRR: This code in plumed is heavily modified from "
285 : "ExtendedLagrangian.cpp and doesn't implemented all variants of "
286 4 : "eABF/DRR."
287 8 : << '\n';
288 4 : log << "eABF/DRR: The thermostat using here maybe different from colvars."
289 8 : << '\n';
290 : log << "eABF/DRR: To integrate the gradients file, you can use abf_integrate "
291 4 : "from https://github.com/colvars/colvars/tree/master/colvartools."
292 8 : << '\n';
293 : log << "eABF/DRR: Please reading relevant articles and using this bias "
294 4 : "method carefully!"
295 8 : << '\n';
296 4 : parseFlag("NOBIAS", nobias);
297 4 : parseFlag("UI", useUIestimator);
298 4 : bool noCZAR = false;
299 4 : parseFlag("NOCZAR", noCZAR);
300 4 : noCZAR == false ? useCZARestimator = true : useCZARestimator = false;
301 4 : parseFlag("TEXTOUTPUT", textoutput);
302 4 : parseVector("TAU", tau);
303 4 : parseVector("FRICTION", friction);
304 4 : parseVector("EXTTEMP", etemp);
305 4 : parseVector("KAPPA", kappa);
306 4 : double temp = -1.0;
307 4 : parse("TEMP", temp);
308 4 : parse("FULLSAMPLES", fullsamples);
309 4 : parse("OUTPUTFREQ", outputfreq);
310 4 : parse("HISTORYFREQ", historyfreq);
311 4 : parse("OUTPUTPREFIX", outputprefix);
312 4 : string restartfilename;
313 4 : parse("DRR_RFILE", restartfilename);
314 8 : string uirprefix;
315 4 : parse("UIRESTARTPREFIX", uirprefix);
316 4 : if (temp >= 0.0)
317 4 : kbt = plumed.getAtoms().getKBoltzmann() * temp;
318 : else
319 0 : kbt = plumed.getAtoms().getKbT();
320 4 : if (fullsamples < 0.5) {
321 0 : fullsamples = 500.0;
322 : log << "eABF/DRR: The fullsamples parametre is not set. Set it to "
323 0 : "500(default)."
324 0 : << '\n';
325 : }
326 4 : if (getRestart()) {
327 1 : if (restartfilename.length() != 0) {
328 1 : isRestart = true;
329 1 : firsttime = false;
330 1 : load(restartfilename);
331 : } else {
332 0 : log << "eABF/DRR: You don't specify the file for restarting." << '\n';
333 0 : log << "eABF/DRR: So I assume you are splitting windows." << '\n';
334 0 : isRestart = false;
335 0 : firsttime = true;
336 : }
337 : }
338 :
339 8 : vector<string> gmin(ndims);
340 4 : parseVector("GRID_MIN", gmin);
341 4 : if (gmin.size() != ndims)
342 0 : error("eABF/DRR: not enough values for GRID_MIN");
343 8 : vector<string> gmax(ndims);
344 4 : parseVector("GRID_MAX", gmax);
345 4 : if (gmax.size() != ndims)
346 0 : error("eABF/DRR: not enough values for GRID_MAX");
347 8 : vector<unsigned> gbin(ndims);
348 8 : vector<double> gspacing(ndims);
349 4 : parseVector("GRID_BIN", gbin);
350 4 : parseVector("GRID_SPACING", gspacing);
351 4 : if (gbin.size() != ndims) {
352 : log << "eABF/DRR: You didn't specify GRID_BIN. Trying to use GRID_SPACING "
353 2 : "instead."
354 4 : << '\n';
355 2 : if (gspacing.size() != ndims) {
356 0 : error("eABF/DRR: not enough values for GRID_BIN");
357 : } else {
358 2 : gbin.resize(ndims);
359 4 : for (size_t i = 0; i < ndims; ++i) {
360 : double l, h;
361 2 : PLMD::Tools::convert(gmin[i], l);
362 2 : PLMD::Tools::convert(gmax[i], h);
363 2 : gbin[i] = std::nearbyint((h - l) / gspacing[i]);
364 2 : gspacing[i] = (h - l) / gbin[i];
365 2 : log << "GRID_BIN[" << i << "] is " << gbin[i] << '\n';
366 : }
367 : }
368 : }
369 4 : checkRead();
370 :
371 : // Set up kbt for extended system
372 4 : log << "eABF/DRR: The fullsamples is " << fullsamples << '\n';
373 4 : log << "eABF/DRR: The kbt(real system) is " << kbt << '\n';
374 4 : dt = getTimeStep();
375 8 : vector<double> ekbt(ndims, 0.0);
376 4 : if (etemp.size() != ndims) {
377 4 : etemp.assign(ndims, kbt / plumed.getAtoms().getKBoltzmann());
378 : }
379 4 : if (tau.size() != ndims) {
380 0 : tau.assign(ndims, 0.5);
381 : }
382 4 : if (friction.size() != ndims) {
383 0 : friction.assign(ndims, 8.0);
384 : }
385 10 : for (size_t i = 0; i < ndims; ++i) {
386 6 : ekbt[i] = etemp[i] * plumed.getAtoms().getKBoltzmann();
387 6 : log << "eABF/DRR: The kbt(extended system) of [" << i << "] is " << ekbt[i]
388 12 : << '\n';
389 6 : log << "eABF/DRR: relaxation time tau [" << i << "] is " << tau[i] << '\n';
390 6 : log << "eABF/DRR: Extended variable [" << i << "] has friction: " << friction[i] << '\n';
391 : }
392 :
393 : // Set up the force grid
394 10 : for (size_t i = 0; i < ndims; ++i) {
395 6 : log << "eABF/DRR: The " << i << " dimensional grid minimum is " << gmin[i]
396 12 : << '\n';
397 6 : log << "eABF/DRR: The " << i << " dimensional grid maximum is " << gmax[i]
398 12 : << '\n';
399 6 : log << "eABF/DRR: The " << i << " dimensional grid has " << gbin[i]
400 12 : << " bins" << '\n';
401 : double l, h;
402 6 : PLMD::Tools::convert(gmin[i], l);
403 6 : PLMD::Tools::convert(gmax[i], h);
404 6 : delim[i].set(l, h, gbin[i]);
405 : }
406 4 : if (kappa.size() != ndims) {
407 2 : kappa.resize(ndims, 0.0);
408 6 : for (size_t i = 0; i < ndims; ++i) {
409 4 : if (kappa[i] <= 0) {
410 4 : log << "eABF/DRR: The spring force constant kappa[" << i
411 8 : << "] is not set." << '\n';
412 4 : kappa[i] = ekbt[i] / (delim[i].getWidth() * delim[i].getWidth());
413 4 : log << "eABF/DRR: set kappa[" << i
414 8 : << "] according to bin width(ekbt/(binWidth^2))." << '\n';
415 : }
416 4 : log << "eABF/DRR: The spring force constant kappa[" << i << "] is "
417 8 : << std::fixed << std::setprecision(10) << kappa[i] << '\n';
418 : }
419 : } else {
420 2 : log << "eABF/DRR: The kappa have been set manually." << '\n';
421 4 : for (size_t i = 0; i < ndims; ++i) {
422 2 : log << "eABF/DRR: The spring force constant kappa[" << i << "] is "
423 4 : << std::fixed << std::setprecision(10) << kappa[i] << '\n';
424 : }
425 : }
426 :
427 10 : for (size_t i = 0; i < ndims; ++i) {
428 6 : mass[i] = kappa[i] * tau[i] * tau[i] / (4 * pi * pi);
429 6 : log << "eABF/DRR: Fictitious mass[" << i << "] is " << mass[i] << '\n';
430 6 : c1[i] = exp(-0.5 * friction[i] * dt);
431 6 : c2[i] = sqrt(ekbt[i] * (1.0 - c1[i] * c1[i]) / mass[i]);
432 : }
433 :
434 10 : for (size_t i = 0; i < ndims; ++i) {
435 : // Position output
436 6 : string comp = getPntrToArgument(i)->getName() + "_fict";
437 6 : addComponentWithDerivatives(comp);
438 6 : if (getPntrToArgument(i)->isPeriodic()) {
439 8 : string a, b;
440 : double c, d;
441 4 : getPntrToArgument(i)->getDomain(a, b);
442 4 : getPntrToArgument(i)->getDomain(c, d);
443 4 : componentIsPeriodic(comp, a, b);
444 8 : delim[i].setPeriodicity(c, d);
445 : } else
446 2 : componentIsNotPeriodic(comp);
447 6 : fictValue[i] = getPntrToComponent(comp);
448 : // Velocity output
449 6 : comp = getPntrToArgument(i)->getName() + "_vfict";
450 6 : addComponent(comp);
451 6 : componentIsNotPeriodic(comp);
452 6 : vfictValue[i] = getPntrToComponent(comp);
453 : // Bias force from eABF/DRR output
454 6 : comp = getPntrToArgument(i)->getName() + "_biasforce";
455 6 : addComponent(comp);
456 6 : componentIsNotPeriodic(comp);
457 6 : biasforceValue[i] = getPntrToComponent(comp);
458 6 : }
459 :
460 4 : if (outputprefix.length() == 0)
461 4 : outputprefix = getLabel();
462 4 : outputname = outputprefix + ".drrstate";
463 4 : cptname = outputprefix + ".cpt.drrstate";
464 :
465 4 : if (!isRestart) {
466 : // If you want to use on-the-fly text output for CZAR and naive estimator,
467 : // you should turn it to true first!
468 3 : ABFGrid = ABF(delim, ".abf", textoutput);
469 : // Just initialize it even useCZARestimator is off.
470 3 : CZARestimator = CZAR(delim, ".czar", kbt, textoutput);
471 3 : log << "eABF/DRR: The init function of the grid is finished." << '\n';
472 : }
473 4 : if (useCZARestimator) {
474 3 : log << "eABF/DRR: Using corrected z-average restraint estimator of gradients" << '\n';
475 3 : log << " Bibliography " << plumed.cite("Lesage, Lelièvre, Stoltz and Hénin, "
476 6 : "J. Phys. Chem. B 3676, 121 (2017)");
477 3 : log << plumed.cite("Darve and Pohorille, J. Chem. Phys. 9169, 115 (2001)") << '\n';
478 : }
479 4 : if (useUIestimator) {
480 : log << "eABF/DRR: Using umbrella integration(Zheng and Yang's) estimator "
481 4 : "of gradients."
482 8 : << '\n';
483 4 : log << "eABF/DRR: The UI estimator code is contributed by Haohao Fu."
484 8 : << '\n';
485 4 : log << " Bibliography " << plumed.cite(
486 8 : "Fu, Shao, Chipot and Cai, J. Chem. Theory Comput. 3506, 12 (2016)");
487 4 : log << plumed.cite("Zheng and Yang, J. Chem. Theory Comput. 810, 8 (2012)");
488 4 : log << plumed.cite("Darve and Pohorille, J. Chem. Phys. 9169, 115 (2001)") << '\n';
489 4 : vector<double> lowerboundary(delim.size(), 0);
490 8 : vector<double> upperboundary(delim.size(), 0);
491 8 : vector<double> width(delim.size(), 0);
492 10 : for (size_t i = 0; i < delim.size(); ++i) {
493 6 : lowerboundary[i] = delim[i].getMin();
494 6 : upperboundary[i] = delim[i].getMax();
495 6 : width[i] = delim[i].getWidth();
496 : }
497 8 : vector<string> input_filename;
498 4 : bool uirestart = false;
499 4 : if (isRestart && (uirprefix.length() != 0)) {
500 1 : input_filename.push_back(uirprefix);
501 1 : uirestart = true;
502 : }
503 4 : if (isRestart && (uirprefix.length() == 0)) {
504 0 : input_filename.push_back(getLabel());
505 : }
506 16 : eabf_UI = UIestimator::UIestimator(
507 4 : lowerboundary, upperboundary, width, kappa, getLabel(), int(outputfreq),
508 12 : uirestart, input_filename, kbt / plumed.getAtoms().getKBoltzmann());
509 4 : }
510 4 : }
511 :
512 34 : void DynamicReferenceRestraining::calculate() {
513 34 : long long int step_now = getStep();
514 34 : if (firsttime) {
515 7 : for (size_t i = 0; i < ndims; ++i) {
516 4 : fict[i] = getArgument(i);
517 : }
518 3 : firsttime = false;
519 : }
520 34 : if (step_now != 0) {
521 30 : if ((step_now % int(outputfreq)) == 0) {
522 6 : save(outputname, step_now);
523 6 : if (textoutput) {
524 4 : ABFGrid.writeAll(outputprefix);
525 4 : if (useCZARestimator) {
526 2 : CZARestimator.writeAll(outputprefix);
527 : }
528 : }
529 : }
530 30 : if (historyfreq > 0 && (step_now % int(historyfreq)) == 0) {
531 : const string filename =
532 0 : outputprefix + "." + std::to_string(step_now) + ".drrstate";
533 0 : save(filename, step_now);
534 0 : if (textoutput) {
535 : const string textfilename =
536 0 : outputprefix + "." + std::to_string(step_now);
537 0 : ABFGrid.writeAll(textfilename);
538 0 : if (useCZARestimator) {
539 0 : CZARestimator.writeAll(textfilename);
540 0 : }
541 0 : }
542 : }
543 30 : if (getCPT()) {
544 : log << "eABF/DRR: The MD engine is writing checkpoint so we also write a "
545 0 : "DRR state file at step: "
546 0 : << step_now << ".\n";
547 0 : save(cptname, step_now);
548 : }
549 : }
550 34 : double ene = 0.0;
551 92 : for (size_t i = 0; i < ndims; ++i) {
552 58 : real[i] = getArgument(i);
553 58 : springlength[i] = difference(i, fict[i], real[i]);
554 58 : fictNoPBC[i] = real[i] - springlength[i];
555 58 : double f = -kappa[i] * springlength[i];
556 58 : ffict_measured[i] = -f;
557 58 : ene += 0.5 * kappa[i] * springlength[i] * springlength[i];
558 58 : setOutputForce(i, f);
559 58 : ffict[i] = -f;
560 58 : fict[i] = fictValue[i]->bringBackInPbc(fict[i]);
561 58 : fictValue[i]->set(fict[i]);
562 58 : vfictValue[i]->set(vfict_laststep[i]);
563 : }
564 34 : setBias(ene);
565 34 : ABFGrid.store_getbias(fict, ffict_measured, fbias, fullsamples);
566 34 : if (useCZARestimator) {
567 29 : CZARestimator.store(real, ffict_measured);
568 : }
569 34 : if (useUIestimator) {
570 34 : eabf_UI.update_output_filename(outputprefix);
571 34 : eabf_UI.update(int(step_now), real, fictNoPBC);
572 : }
573 34 : }
574 :
575 34 : void DynamicReferenceRestraining::update() {
576 92 : for (size_t i = 0; i < ndims; ++i) {
577 : // consider additional forces on the fictitious particle
578 : // (e.g. MetaD stuff)
579 58 : ffict[i] += fictValue[i]->getForce();
580 58 : if (!nobias) {
581 58 : ffict[i] += fbias[i];
582 : }
583 58 : biasforceValue[i]->set(fbias[i]);
584 : // update velocity (half step)
585 58 : vfict[i] += ffict[i] * 0.5 * dt / mass[i];
586 : // thermostat (half step)
587 58 : vfict[i] = c1[i] * vfict[i] + c2[i] * rand.Gaussian();
588 : // save full step velocity to be dumped at next step
589 58 : vfict_laststep[i] = vfict[i];
590 : // thermostat (half step)
591 58 : vfict[i] = c1[i] * vfict[i] + c2[i] * rand.Gaussian();
592 : // update velocity (half step)
593 58 : vfict[i] += ffict[i] * 0.5 * dt / mass[i];
594 : // update position (full step)
595 58 : fict[i] += vfict[i] * dt;
596 : }
597 34 : }
598 :
599 6 : void DynamicReferenceRestraining::save(const string &filename,
600 : long long int step) {
601 6 : std::ofstream out;
602 6 : out.open(filename.c_str(), std::ios::binary);
603 12 : boost::archive::binary_oarchive oa(out);
604 6 : oa << step << fict << vfict << vfict_laststep << ffict << ABFGrid
605 12 : << CZARestimator;
606 12 : out.close();
607 6 : }
608 :
609 1 : void DynamicReferenceRestraining::load(const string &filename) {
610 1 : std::ifstream in;
611 : long long int step;
612 1 : in.open(filename.c_str(), std::ios::binary);
613 1 : log << "eABF/DRR: Read restart file: " << filename << '\n';
614 2 : boost::archive::binary_iarchive ia(in);
615 1 : ia >> step >> fict >> vfict >> vfict_laststep >> ffict >> ABFGrid >>
616 2 : CZARestimator;
617 1 : in.close();
618 1 : log << "eABF/DRR: Restart at step: " << step << '\n';
619 2 : backupFile(filename);
620 1 : }
621 :
622 1 : void DynamicReferenceRestraining::backupFile(const string &filename) {
623 1 : bool isSuccess = false;
624 1 : long int i = 0;
625 3 : while (!isSuccess) {
626 : // If libstdc++ support C++17 we can simplify following code.
627 1 : const string bckname = "bck." + filename + "." + std::to_string(i);
628 1 : if (is_file_exist(bckname.c_str())) {
629 0 : ++i;
630 : } else {
631 1 : log << "eABF/DRR: Backup original restart file to " << bckname << '\n';
632 1 : std::ifstream src(filename.c_str(), std::ios::binary);
633 2 : std::ofstream dst(bckname.c_str(), std::ios::binary);
634 1 : dst << src.rdbuf();
635 1 : src.close();
636 1 : dst.close();
637 2 : isSuccess = true;
638 : }
639 1 : }
640 1 : }
641 :
642 : // Copy from
643 : // stackoverflow(https://stackoverflow.com/questions/12774207/fastest-way-to-check-if-a-file-exist-using-standard-c-c11-c)
644 1 : bool DynamicReferenceRestraining::is_file_exist(const char *fileName) {
645 1 : std::ifstream infile(fileName);
646 1 : return infile.good();
647 : }
648 : }
649 4821 : }
650 :
651 : #endif
|
