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 16 : 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 6456 : PLUMED_REGISTER_ACTION(DynamicReferenceRestraining, "DRR")
184 :
185 5 : void DynamicReferenceRestraining::registerKeywords(Keywords &keys) {
186 5 : Bias::registerKeywords(keys);
187 10 : keys.use("ARG");
188 20 : 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 : "kbt/(GRID_SPACING)^2)");
192 25 : keys.add("compulsory", "TAU", "0.5", "specifies relaxation time on each of "
193 : "variables are, similar to "
194 : "extendedTimeConstant in Colvars");
195 25 : keys.add("compulsory", "FRICTION", "8.0",
196 : "add a friction to the variable, similar to extendedLangevinDamping "
197 : "in Colvars");
198 20 : keys.add("compulsory", "GRID_MIN", "the lower bounds for the grid (GRID_BIN "
199 : "or GRID_SPACING should be specified)");
200 20 : keys.add("compulsory", "GRID_MAX", "the upper bounds for the grid (GRID_BIN "
201 : "or GRID_SPACING should be specified)");
202 20 : keys.add("optional", "GRID_BIN", "the number of bins for the grid");
203 20 : keys.add("optional", "GRID_SPACING", "the approximate grid spacing (to be "
204 : "used as an alternative or together "
205 : "with GRID_BIN)");
206 25 : keys.add("compulsory", "FULLSAMPLES", "500",
207 : "number of samples in a bin prior to application of the ABF");
208 20 : keys.add("compulsory", "OUTPUTFREQ", "write results to a file every N steps");
209 20 : keys.add("optional", "HISTORYFREQ", "save history to a file every N steps");
210 15 : keys.addFlag("NOCZAR", false, "disable the CZAR estimator");
211 15 : keys.addFlag("UI", false,
212 : "enable the umbrella integration estimator");
213 20 : keys.add("optional", "UIRESTARTPREFIX",
214 : "specify the restart files for umbrella integration");
215 20 : keys.add("optional", "OUTPUTPREFIX",
216 : "specify the output prefix (default to the label name)");
217 20 : keys.add("optional", "TEMP", "the system temperature - needed when FRICTION "
218 : "is present. If not provided will be taken from "
219 : "MD code (if available)");
220 20 : keys.add(
221 : "optional", "EXTTEMP",
222 : "the temperature of extended variables (default to system temperature)");
223 20 : keys.add("optional", "DRR_RFILE",
224 : "specifies the restart file (.drrstate file)");
225 15 : keys.addFlag("NOBIAS", false, "DO NOT apply bias forces.");
226 15 : keys.addFlag("TEXTOUTPUT", false, "use text output for grad and count files "
227 : "instead of boost::serialization binary "
228 : "output");
229 5 : componentsAreNotOptional(keys);
230 20 : 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 : "variable.");
237 20 : 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 : "variable.");
244 20 : keys.addOutputComponent(
245 : "_biasforce", "default",
246 : "The bias force from eABF/DRR of the fictitious particle.");
247 5 : }
248 :
249 4 : DynamicReferenceRestraining::DynamicReferenceRestraining(
250 4 : const ActionOptions &ao)
251 : : PLUMED_BIAS_INIT(ao), firsttime(true), nobias(false),
252 : fictNoPBC(getNumberOfArguments(), 0.0), real(getNumberOfArguments(), 0.0),
253 : springlength(getNumberOfArguments(), 0.0),
254 : fict(getNumberOfArguments(), 0.0), vfict(getNumberOfArguments(), 0.0),
255 : vfict_laststep(getNumberOfArguments(), 0.0),
256 : ffict(getNumberOfArguments(), 0.0), fbias(getNumberOfArguments(), 0.0),
257 : kappa(getNumberOfArguments(), 0.0), tau(getNumberOfArguments(), 0.0),
258 : friction(getNumberOfArguments(), 0.0), etemp(getNumberOfArguments(), 0.0),
259 : ffict_measured(getNumberOfArguments(), 0.0),
260 : biasforceValue(getNumberOfArguments(), NULL),
261 : fictValue(getNumberOfArguments(), NULL),
262 : vfictValue(getNumberOfArguments(), NULL), c1(getNumberOfArguments(), 0.0),
263 : c2(getNumberOfArguments(), 0.0), mass(getNumberOfArguments(), 0.0),
264 : 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 92 : useCZARestimator(true), useUIestimator(false), textoutput(false)
268 : {
269 4 : log << "eABF/DRR: You now are using the extended adaptive biasing "
270 : "force(eABF) method."
271 4 : << '\n';
272 4 : log << "eABF/DRR: Some people also refer to it as dynamic reference "
273 : "restraining(DRR) method."
274 4 : << '\n';
275 4 : log << "eABF/DRR: Currently the CZAR and naive(ABF on extended variables) "
276 : "estimator is enabled by default."
277 4 : << '\n';
278 4 : log << "eABF/DRR: For reasons of performance, the umbrella integration "
279 : "estimator is not enabled by default."
280 4 : << '\n';
281 4 : log << "eABF/DRR: This method is originally implemented in "
282 : "colvars(https://github.com/colvars/colvars)."
283 4 : << '\n';
284 4 : log << "eABF/DRR: This code in plumed is heavily modified from "
285 : "ExtendedLagrangian.cpp and doesn't implemented all variants of "
286 : "eABF/DRR."
287 4 : << '\n';
288 4 : log << "eABF/DRR: The thermostat using here maybe different from colvars."
289 4 : << '\n';
290 4 : log << "eABF/DRR: To integrate the gradients file, you can use abf_integrate "
291 : "from https://github.com/colvars/colvars/tree/master/colvartools."
292 4 : << '\n';
293 4 : log << "eABF/DRR: Please reading relevant articles and using this bias "
294 : "method carefully!"
295 4 : << '\n';
296 8 : parseFlag("NOBIAS", nobias);
297 8 : parseFlag("UI", useUIestimator);
298 4 : bool noCZAR = false;
299 8 : parseFlag("NOCZAR", noCZAR);
300 4 : noCZAR == false ? useCZARestimator = true : useCZARestimator = false;
301 8 : parseFlag("TEXTOUTPUT", textoutput);
302 8 : parseVector("TAU", tau);
303 8 : parseVector("FRICTION", friction);
304 8 : parseVector("EXTTEMP", etemp);
305 8 : parseVector("KAPPA", kappa);
306 4 : double temp = -1.0;
307 8 : parse("TEMP", temp);
308 8 : parse("FULLSAMPLES", fullsamples);
309 8 : parse("OUTPUTFREQ", outputfreq);
310 8 : parse("HISTORYFREQ", historyfreq);
311 8 : parse("OUTPUTPREFIX", outputprefix);
312 : string restartfilename;
313 8 : parse("DRR_RFILE", restartfilename);
314 : string uirprefix;
315 8 : parse("UIRESTARTPREFIX", uirprefix);
316 4 : if (temp >= 0.0)
317 8 : 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 0 : log << "eABF/DRR: The fullsamples parametre is not set. Set it to "
323 : "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 8 : 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 8 : parseVector("GRID_MAX", gmax);
345 4 : if (gmax.size() != ndims)
346 0 : error("eABF/DRR: not enough values for GRID_MAX");
347 4 : vector<unsigned> gbin(ndims);
348 4 : vector<double> gspacing(ndims);
349 8 : parseVector("GRID_BIN", gbin);
350 8 : parseVector("GRID_SPACING", gspacing);
351 4 : if (gbin.size() != ndims) {
352 2 : log << "eABF/DRR: You didn't specify GRID_BIN. Trying to use GRID_SPACING "
353 : "instead."
354 2 : << '\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 4 : PLMD::Tools::convert(gmax[i], h);
363 6 : gbin[i] = std::nearbyint((h - l) / gspacing[i]);
364 6 : gspacing[i] = (h - l) / gbin[i];
365 4 : 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 4 : vector<double> ekbt(ndims, 0.0);
376 4 : if (etemp.size() != ndims) {
377 12 : 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 18 : 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 12 : log << "eABF/DRR: relaxation time tau [" << i << "] is " << tau[i] << '\n';
390 12 : 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 12 : PLMD::Tools::convert(gmin[i], l);
403 12 : PLMD::Tools::convert(gmax[i], h);
404 12 : 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 4 : << "] is not set." << '\n';
412 16 : kappa[i] = ekbt[i] / (delim[i].getWidth() * delim[i].getWidth());
413 4 : log << "eABF/DRR: set kappa[" << i
414 4 : << "] 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 18 : mass[i] = kappa[i] * tau[i] * tau[i] / (4 * pi * pi);
429 12 : log << "eABF/DRR: Fictitious mass[" << i << "] is " << mass[i] << '\n';
430 18 : c1[i] = exp(-0.5 * friction[i] * dt);
431 30 : c2[i] = sqrt(ekbt[i] * (1.0 - c1[i] * c1[i]) / mass[i]);
432 : }
433 :
434 16 : 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 : 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 4 : delim[i].setPeriodicity(c, d);
445 : } else
446 2 : componentIsNotPeriodic(comp);
447 6 : fictValue[i] = getPntrToComponent(comp);
448 : // Velocity output
449 12 : 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 12 : comp = getPntrToArgument(i)->getName() + "_biasforce";
455 6 : addComponent(comp);
456 6 : componentIsNotPeriodic(comp);
457 6 : biasforceValue[i] = getPntrToComponent(comp);
458 : }
459 :
460 4 : if (outputprefix.length() == 0)
461 0 : outputprefix = getLabel();
462 8 : outputname = outputprefix + ".drrstate";
463 8 : 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 9 : ABFGrid = ABF(delim, ".abf", textoutput);
469 : // Just initialize it even useCZARestimator is off.
470 9 : 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 12 : log << " Bibliography " << plumed.cite("Lesage, Lelièvre, Stoltz and Hénin, "
476 3 : "J. Phys. Chem. B 3676, 121 (2017)");
477 9 : log << plumed.cite("Darve and Pohorille, J. Chem. Phys. 9169, 115 (2001)") << '\n';
478 : }
479 4 : if (useUIestimator) {
480 4 : log << "eABF/DRR: Using umbrella integration(Zheng and Yang's) estimator "
481 : "of gradients."
482 4 : << '\n';
483 4 : log << "eABF/DRR: The UI estimator code is contributed by Haohao Fu."
484 4 : << '\n';
485 16 : log << " Bibliography " << plumed.cite(
486 4 : "Fu, Shao, Chipot and Cai, J. Chem. Theory Comput. 3506, 12 (2016)");
487 12 : log << plumed.cite("Zheng and Yang, J. Chem. Theory Comput. 810, 8 (2012)");
488 12 : log << plumed.cite("Darve and Pohorille, J. Chem. Phys. 9169, 115 (2001)") << '\n';
489 8 : vector<double> lowerboundary(delim.size(), 0);
490 8 : vector<double> upperboundary(delim.size(), 0);
491 8 : vector<double> width(delim.size(), 0);
492 16 : 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 4 : vector<string> input_filename;
498 : bool uirestart = false;
499 5 : if (isRestart && (uirprefix.length() != 0)) {
500 1 : input_filename.push_back(uirprefix);
501 : uirestart = true;
502 : }
503 5 : if (isRestart && (uirprefix.length() == 0)) {
504 0 : input_filename.push_back(getLabel());
505 : }
506 12 : eabf_UI = UIestimator::UIestimator(
507 4 : lowerboundary, upperboundary, width, kappa, getLabel(), int(outputfreq),
508 8 : uirestart, input_filename, kbt / plumed.getAtoms().getKBoltzmann());
509 : }
510 4 : }
511 :
512 34 : void DynamicReferenceRestraining::calculate() {
513 34 : long long int step_now = getStep();
514 34 : if (firsttime) {
515 11 : 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 : }
541 : }
542 : }
543 30 : if (getCPT()) {
544 0 : log << "eABF/DRR: The MD engine is writing checkpoint so we also write a "
545 : "DRR state file at step: "
546 0 : << step_now << ".\n";
547 0 : save(cptname, step_now);
548 : }
549 : }
550 : double ene = 0.0;
551 150 : for (size_t i = 0; i < ndims; ++i) {
552 58 : real[i] = getArgument(i);
553 174 : springlength[i] = difference(i, fict[i], real[i]);
554 174 : fictNoPBC[i] = real[i] - springlength[i];
555 116 : double f = -kappa[i] * springlength[i];
556 58 : ffict_measured[i] = -f;
557 116 : ene += 0.5 * kappa[i] * springlength[i] * springlength[i];
558 : setOutputForce(i, f);
559 58 : ffict[i] = -f;
560 174 : fict[i] = fictValue[i]->bringBackInPbc(fict[i]);
561 116 : fictValue[i]->set(fict[i]);
562 116 : vfictValue[i]->set(vfict_laststep[i]);
563 : }
564 : 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 102 : eabf_UI.update(int(step_now), real, fictNoPBC);
572 : }
573 34 : }
574 :
575 34 : void DynamicReferenceRestraining::update() {
576 150 : for (size_t i = 0; i < ndims; ++i) {
577 : // consider additional forces on the fictitious particle
578 : // (e.g. MetaD stuff)
579 116 : ffict[i] += fictValue[i]->getForce();
580 58 : if (!nobias) {
581 116 : ffict[i] += fbias[i];
582 : }
583 116 : biasforceValue[i]->set(fbias[i]);
584 : // update velocity (half step)
585 174 : vfict[i] += ffict[i] * 0.5 * dt / mass[i];
586 : // thermostat (half step)
587 232 : 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 232 : vfict[i] = c1[i] * vfict[i] + c2[i] * rand.Gaussian();
592 : // update velocity (half step)
593 174 : vfict[i] += ffict[i] * 0.5 * dt / mass[i];
594 : // update position (full step)
595 116 : fict[i] += vfict[i] * dt;
596 : }
597 34 : }
598 :
599 6 : void DynamicReferenceRestraining::save(const string &filename,
600 : long long int step) {
601 12 : std::ofstream out;
602 6 : out.open(filename.c_str(), std::ios::binary);
603 : boost::archive::binary_oarchive oa(out);
604 30 : oa << step << fict << vfict << vfict_laststep << ffict << ABFGrid
605 6 : << CZARestimator;
606 6 : out.close();
607 6 : }
608 :
609 1 : void DynamicReferenceRestraining::load(const string &filename) {
610 2 : 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 : boost::archive::binary_iarchive ia(in);
615 5 : ia >> step >> fict >> vfict >> vfict_laststep >> ffict >> ABFGrid >>
616 1 : CZARestimator;
617 1 : in.close();
618 1 : log << "eABF/DRR: Restart at step: " << step << '\n';
619 1 : backupFile(filename);
620 1 : }
621 :
622 1 : void DynamicReferenceRestraining::backupFile(const string &filename) {
623 : bool isSuccess = false;
624 : long int i = 0;
625 3 : while (!isSuccess) {
626 : // If libstdc++ support C++17 we can simplify following code.
627 4 : 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 2 : 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 : isSuccess = true;
638 : }
639 : }
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 2 : std::ifstream infile(fileName);
646 1 : return infile.good();
647 : }
648 : }
649 4839 : }
650 :
651 : #endif
|
