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 just 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 : It's also possible to run extended generalized adaptive biasing force (egABF) described in \cite Zhao2017 .
125 : An egABF example:
126 : \plumedfile
127 : phi: TORSION ATOMS=5,7,9,15
128 : psi: TORSION ATOMS=7,9,15,17
129 :
130 : DRR ...
131 : LABEL=gabf_phi
132 : ARG=phi
133 : FULLSAMPLES=500
134 : GRID_MIN=-pi
135 : GRID_MAX=pi
136 : GRID_BIN=180
137 : FRICTION=8.0
138 : TAU=0.5
139 : OUTPUTFREQ=50000
140 : HISTORYFREQ=500000
141 : ... DRR
142 :
143 : DRR ...
144 : LABEL=gabf_psi
145 : ARG=psi
146 : FULLSAMPLES=500
147 : GRID_MIN=-pi
148 : GRID_MAX=pi
149 : GRID_BIN=180
150 : FRICTION=8.0
151 : TAU=0.5
152 : OUTPUTFREQ=50000
153 : HISTORYFREQ=500000
154 : ... DRR
155 :
156 : DRR ...
157 : LABEL=gabf_2d
158 : ARG=phi,psi
159 : EXTERNAL_FORCE=gabf_phi.phi_springforce,gabf_psi.psi_springforce
160 : EXTERNAL_FICT=gabf_phi.phi_fictNoPBC,gabf_psi.psi_fictNoPBC
161 : GRID_MIN=-pi,-pi
162 : GRID_MAX=pi,pi
163 : GRID_BIN=180,180
164 : NOBIAS
165 : OUTPUTFREQ=50000
166 : HISTORYFREQ=500000
167 : ... DRR
168 :
169 : PRINT STRIDE=10 ARG=phi,psi FILE=COLVAR
170 : \endplumedfile
171 :
172 : */
173 : //+ENDPLUMEDOC
174 :
175 : using std::vector;
176 : using std::string;
177 :
178 40 : class DynamicReferenceRestraining : public Bias {
179 : private:
180 : bool firsttime;
181 : bool nobias;
182 : vector<double> fictNoPBC;
183 : vector<double> real;
184 : vector<double> springlength; // spring lengths
185 : vector<double> fict; // coordinates of extended variables
186 : vector<double> vfict; // velocities of extended variables
187 : vector<double> vfict_laststep;
188 : vector<double> ffict; // forces exerted on extended variables
189 : vector<double> fbias; // bias forces from eABF
190 : vector<double> kappa;
191 : vector<double> tau;
192 : vector<double> friction;
193 : vector<double> etemp;
194 : vector<double> ffict_measured;
195 : vector<double> force_external;
196 : vector<double> fict_external;
197 : vector<Value *> biasforceValue;
198 : vector<Value *> springforceValue;
199 : vector<Value *> fictValue;
200 : vector<Value *> vfictValue;
201 : vector<Value *> fictNoPBCValue;
202 : vector<Value *> externalForceValue;
203 : vector<Value *> externalFictValue;
204 : vector<double> c1;
205 : vector<double> c2;
206 : vector<double> mass;
207 : vector<DRRAxis> delim;
208 : string outputname;
209 : string cptname;
210 : string outputprefix;
211 : const size_t ndims;
212 : double dt;
213 : double kbt;
214 : double outputfreq;
215 : double historyfreq;
216 : bool isRestart;
217 : bool useCZARestimator;
218 : bool useUIestimator;
219 : bool textoutput;
220 : bool withExternalForce;
221 : bool withExternalFict;
222 : ABF ABFGrid;
223 : CZAR CZARestimator;
224 : double fullsamples;
225 : double maxFactor;
226 : UIestimator::UIestimator eabf_UI;
227 : Random rand;
228 :
229 : public:
230 : explicit DynamicReferenceRestraining(const ActionOptions &);
231 : void calculate();
232 : void update();
233 : void save(const string &filename, long long int step);
234 : void load(const string &filename);
235 : void backupFile(const string &filename);
236 : static void registerKeywords(Keywords &keys);
237 : bool is_file_exist(const char *fileName);
238 : };
239 :
240 7376 : PLUMED_REGISTER_ACTION(DynamicReferenceRestraining, "DRR")
241 :
242 11 : void DynamicReferenceRestraining::registerKeywords(Keywords &keys) {
243 11 : Bias::registerKeywords(keys);
244 22 : keys.use("ARG");
245 44 : keys.add("optional", "KAPPA", "specifies that the restraint is harmonic and "
246 : "what the values of the force constants on "
247 : "each of the variables are (default to "
248 : "\\f$k_BT\\f$/(GRID_SPACING)^2)");
249 55 : keys.add("compulsory", "TAU", "0.5", "specifies relaxation time on each of "
250 : "variables are, similar to "
251 : "extended Time Constant in Colvars");
252 55 : keys.add("compulsory", "FRICTION", "8.0",
253 : "add a friction to the variable, similar to extended Langevin Damping "
254 : "in Colvars");
255 44 : keys.add("compulsory", "GRID_MIN", "the lower bounds for the grid (GRID_BIN "
256 : "or GRID_SPACING should be specified)");
257 44 : keys.add("compulsory", "GRID_MAX", "the upper bounds for the grid (GRID_BIN "
258 : "or GRID_SPACING should be specified)");
259 44 : keys.add("optional", "GRID_BIN", "the number of bins for the grid");
260 44 : keys.add("optional", "GRID_SPACING", "the approximate grid spacing (to be "
261 : "used as an alternative or together "
262 : "with GRID_BIN)");
263 44 : keys.add("optional", "ZGRID_MIN", "the lower bounds for the grid (ZGRID_BIN"
264 : " or ZGRID_SPACING should be specified)");
265 44 : keys.add("optional", "ZGRID_MAX", "the upper bounds for the grid (ZGRID_BIN"
266 : " or ZGRID_SPACING should be specified)");
267 44 : keys.add("optional", "ZGRID_BIN", "the number of bins for the grid");
268 44 : keys.add("optional", "ZGRID_SPACING", "the approximate grid spacing (to be "
269 : "used as an alternative or together "
270 : "with ZGRID_BIN)");
271 44 : keys.add("optional", "EXTERNAL_FORCE", "use forces from other action instead"
272 : " of internal spring force, this disable the extended system!");
273 44 : keys.add("optional", "EXTERNAL_FICT", "position of external fictitious "
274 : "particles, useful for UIESTIMATOR");
275 55 : keys.add("compulsory", "FULLSAMPLES", "500",
276 : "number of samples in a bin prior to application of the ABF");
277 55 : keys.add("compulsory", "MAXFACTOR", "1.0",
278 : "maximum scaling factor of biasing force");
279 44 : keys.add("compulsory", "OUTPUTFREQ", "write results to a file every N steps");
280 44 : keys.add("optional", "HISTORYFREQ", "save history to a file every N steps");
281 33 : keys.addFlag("NOCZAR", false, "disable the CZAR estimator");
282 33 : keys.addFlag("UI", false,
283 : "enable the umbrella integration estimator");
284 44 : keys.add("optional", "UIRESTARTPREFIX",
285 : "specify the restart files for umbrella integration");
286 44 : keys.add("optional", "OUTPUTPREFIX",
287 : "specify the output prefix (default to the label name)");
288 44 : keys.add("optional", "TEMP", "the system temperature - needed when FRICTION "
289 : "is present. If not provided will be taken from "
290 : "MD code (if available)");
291 44 : keys.add(
292 : "optional", "EXTTEMP",
293 : "the temperature of extended variables (default to system temperature)");
294 44 : keys.add("optional", "DRR_RFILE",
295 : "specifies the restart file (.drrstate file)");
296 33 : keys.addFlag("NOBIAS", false, "DO NOT apply bias forces.");
297 33 : keys.addFlag("TEXTOUTPUT", false, "use text output for grad and count files "
298 : "instead of boost::serialization binary "
299 : "output");
300 11 : componentsAreNotOptional(keys);
301 44 : keys.addOutputComponent(
302 : "_fict", "default",
303 : "one or multiple instances of this quantity can be referenced "
304 : "elsewhere in the input file. "
305 : "These quantities will named with the arguments of the bias followed by "
306 : "the character string _tilde. It is possible to add forces on these "
307 : "variable.");
308 44 : keys.addOutputComponent(
309 : "_vfict", "default",
310 : "one or multiple instances of this quantity can be referenced "
311 : "elsewhere in the input file. "
312 : "These quantities will named with the arguments of the bias followed by "
313 : "the character string _tilde. It is NOT possible to add forces on these "
314 : "variable.");
315 44 : keys.addOutputComponent(
316 : "_biasforce", "default",
317 : "The bias force from eABF/DRR of the fictitious particle.");
318 44 : keys.addOutputComponent("_springforce", "default", "Spring force between real CVs and extended CVs");
319 44 : keys.addOutputComponent("_fictNoPBC", "default",
320 : "the positions of fictitious particles (without PBC).");
321 11 : }
322 :
323 10 : DynamicReferenceRestraining::DynamicReferenceRestraining(
324 10 : const ActionOptions &ao)
325 : : PLUMED_BIAS_INIT(ao), firsttime(true), nobias(false),
326 : fictNoPBC(getNumberOfArguments(), 0.0), real(getNumberOfArguments(), 0.0),
327 : springlength(getNumberOfArguments(), 0.0),
328 : fict(getNumberOfArguments(), 0.0), vfict(getNumberOfArguments(), 0.0),
329 : vfict_laststep(getNumberOfArguments(), 0.0),
330 : ffict(getNumberOfArguments(), 0.0), fbias(getNumberOfArguments(), 0.0),
331 : kappa(getNumberOfArguments(), 0.0), tau(getNumberOfArguments(), 0.0),
332 : friction(getNumberOfArguments(), 0.0), etemp(getNumberOfArguments(), 0.0),
333 : ffict_measured(getNumberOfArguments(), 0.0),
334 : biasforceValue(getNumberOfArguments(), NULL),
335 : springforceValue(getNumberOfArguments(), NULL),
336 : fictValue(getNumberOfArguments(), NULL),
337 : vfictValue(getNumberOfArguments(), NULL),
338 : fictNoPBCValue(getNumberOfArguments(), NULL),
339 : externalForceValue(getNumberOfArguments(), NULL),
340 : externalFictValue(getNumberOfArguments(), NULL),
341 : c1(getNumberOfArguments(), 0.0),
342 : c2(getNumberOfArguments(), 0.0), mass(getNumberOfArguments(), 0.0),
343 : delim(getNumberOfArguments()), outputname(""), cptname(""),
344 10 : outputprefix(""), ndims(getNumberOfArguments()), dt(0.0), kbt(0.0),
345 : outputfreq(0.0), historyfreq(-1.0), isRestart(false),
346 : useCZARestimator(true), useUIestimator(false), textoutput(false),
347 280 : withExternalForce(false), withExternalFict(false)
348 : {
349 10 : log << "eABF/DRR: You now are using the extended adaptive biasing "
350 : "force(eABF) method."
351 10 : << '\n';
352 10 : log << "eABF/DRR: Some people also refer to it as dynamic reference "
353 : "restraining(DRR) method."
354 10 : << '\n';
355 10 : log << "eABF/DRR: Currently the CZAR and naive(ABF on extended variables) "
356 : "estimator is enabled by default."
357 10 : << '\n';
358 10 : log << "eABF/DRR: For reasons of performance, the umbrella integration "
359 : "estimator is not enabled by default."
360 10 : << '\n';
361 10 : log << "eABF/DRR: This method is originally implemented in "
362 : "colvars(https://github.com/colvars/colvars)."
363 10 : << '\n';
364 10 : log << "eABF/DRR: This code in plumed is heavily modified from "
365 : "ExtendedLagrangian.cpp and doesn't implemented all variants of "
366 : "eABF/DRR."
367 10 : << '\n';
368 10 : log << "eABF/DRR: The thermostat using here maybe different from colvars."
369 10 : << '\n';
370 10 : log << "eABF/DRR: To integrate the gradients file, you can use abf_integrate "
371 : "from https://github.com/colvars/colvars/tree/master/colvartools."
372 10 : << '\n';
373 10 : log << "eABF/DRR: Please reading relevant articles and using this bias "
374 : "method carefully!"
375 10 : << '\n';
376 20 : parseFlag("NOBIAS", nobias);
377 20 : parseFlag("UI", useUIestimator);
378 10 : bool noCZAR = false;
379 20 : parseFlag("NOCZAR", noCZAR);
380 : // noCZAR == false ? useCZARestimator = true : useCZARestimator = false;
381 20 : parseFlag("TEXTOUTPUT", textoutput);
382 20 : parseVector("TAU", tau);
383 20 : parseVector("FRICTION", friction);
384 20 : parseVector("EXTTEMP", etemp);
385 20 : parseVector("KAPPA", kappa);
386 10 : double temp = -1.0;
387 20 : parse("TEMP", temp);
388 20 : parse("FULLSAMPLES", fullsamples);
389 20 : parse("MAXFACTOR", maxFactor);
390 20 : parse("OUTPUTFREQ", outputfreq);
391 20 : parse("HISTORYFREQ", historyfreq);
392 20 : parse("OUTPUTPREFIX", outputprefix);
393 : string restart_prefix;
394 20 : parse("DRR_RFILE", restart_prefix);
395 : string uirprefix;
396 20 : parse("UIRESTARTPREFIX", uirprefix);
397 20 : parseArgumentList("EXTERNAL_FORCE", externalForceValue);
398 20 : parseArgumentList("EXTERNAL_FICT", externalFictValue);
399 10 : if (externalForceValue.empty()) {
400 9 : withExternalForce = false;
401 1 : } else if (externalForceValue.size() != ndims) {
402 0 : error("eABF/DRR: Number of forces doesn't match ARGS!");
403 : } else {
404 1 : withExternalForce = true;
405 : }
406 10 : if (withExternalForce && useUIestimator) {
407 1 : if (externalFictValue.empty()) {
408 0 : error("eABF/DRR: No external fictitious particles specified. UI estimator needs it.");
409 1 : } else if(externalFictValue.size() != ndims) {
410 0 : error("eABF/DRR: Number of fictitious particles doesn't match ARGS!");
411 : } else {
412 1 : withExternalFict = true;
413 : }
414 : }
415 10 : if (temp >= 0.0)
416 20 : kbt = plumed.getAtoms().getKBoltzmann() * temp;
417 : else
418 0 : kbt = plumed.getAtoms().getKbT();
419 10 : if (fullsamples < 0.5) {
420 0 : fullsamples = 500.0;
421 0 : log << "eABF/DRR: The fullsamples parametre is not set. Set it to "
422 : "500(default)."
423 0 : << '\n';
424 : }
425 10 : if (getRestart()) {
426 1 : if (restart_prefix.length() != 0) {
427 1 : isRestart = true;
428 1 : firsttime = false;
429 1 : load(restart_prefix);
430 : } else {
431 0 : log << "eABF/DRR: You don't specify the file for restarting." << '\n';
432 0 : log << "eABF/DRR: So I assume you are splitting windows." << '\n';
433 0 : isRestart = false;
434 0 : firsttime = true;
435 : }
436 : }
437 :
438 20 : vector<string> gmin(ndims);
439 20 : vector<string> zgmin(ndims);
440 20 : parseVector("GRID_MIN", gmin);
441 20 : parseVector("ZGRID_MIN", zgmin);
442 10 : if (gmin.size() != ndims)
443 0 : error("eABF/DRR: not enough values for GRID_MIN");
444 10 : if (zgmin.size() != ndims) {
445 9 : log << "eABF/DRR: You didn't specify ZGRID_MIN. " << '\n'
446 9 : << "eABF/DRR: The GRID_MIN will be used instead.";
447 9 : zgmin = gmin;
448 : }
449 20 : vector<string> gmax(ndims);
450 20 : vector<string> zgmax(ndims);
451 20 : parseVector("GRID_MAX", gmax);
452 20 : parseVector("ZGRID_MAX", zgmax);
453 10 : if (gmax.size() != ndims)
454 0 : error("eABF/DRR: not enough values for GRID_MAX");
455 10 : if (zgmax.size() != ndims) {
456 9 : log << "eABF/DRR: You didn't specify ZGRID_MAX. " << '\n'
457 9 : << "eABF/DRR: The GRID_MAX will be used instead.";
458 9 : zgmax = gmax;
459 : }
460 10 : vector<unsigned> gbin(ndims);
461 10 : vector<unsigned> zgbin(ndims);
462 10 : vector<double> gspacing(ndims);
463 10 : vector<double> zgspacing(ndims);
464 20 : parseVector("GRID_BIN", gbin);
465 20 : parseVector("ZGRID_BIN", zgbin);
466 20 : parseVector("GRID_SPACING", gspacing);
467 20 : parseVector("ZGRID_SPACING", zgspacing);
468 10 : if (gbin.size() != ndims) {
469 2 : log << "eABF/DRR: You didn't specify GRID_BIN. Trying to use GRID_SPACING "
470 : "instead."
471 2 : << '\n';
472 2 : if (gspacing.size() != ndims) {
473 0 : error("eABF/DRR: not enough values for GRID_BIN");
474 : } else {
475 2 : gbin.resize(ndims);
476 4 : for (size_t i = 0; i < ndims; ++i) {
477 : double l, h;
478 2 : PLMD::Tools::convert(gmin[i], l);
479 4 : PLMD::Tools::convert(gmax[i], h);
480 6 : gbin[i] = std::nearbyint((h - l) / gspacing[i]);
481 6 : gspacing[i] = (h - l) / gbin[i];
482 4 : log << "GRID_BIN[" << i << "] is " << gbin[i] << '\n';
483 : }
484 : }
485 : }
486 10 : if (zgbin.size() != ndims) {
487 9 : log << "eABF/DRR: You didn't specify ZGRID_BIN. Trying to use ZGRID_SPACING instead." << '\n';
488 9 : if (zgspacing.size() != ndims) {
489 9 : log << "eABF/DRR: You didn't specify ZGRID_SPACING. Trying to use GRID_SPACING or GRID_BIN instead." << '\n';
490 9 : zgbin = gbin;
491 9 : zgspacing = gspacing;
492 : } else {
493 0 : zgbin.resize(ndims);
494 0 : for (size_t i = 0; i < ndims; ++i) {
495 : double l, h;
496 0 : PLMD::Tools::convert(zgmin[i], l);
497 0 : PLMD::Tools::convert(zgmax[i], h);
498 0 : zgbin[i] = std::nearbyint((h - l) / zgspacing[i]);
499 0 : zgspacing[i] = (h - l) / zgbin[i];
500 0 : log << "ZGRID_BIN[" << i << "] is " << zgbin[i] << '\n';
501 : }
502 : }
503 : }
504 10 : checkRead();
505 :
506 : // Set up kbt for extended system
507 10 : log << "eABF/DRR: The fullsamples is " << fullsamples << '\n';
508 10 : log << "eABF/DRR: The maximum scaling factor is " << maxFactor << '\n';
509 10 : if (maxFactor > 1.0) {
510 0 : log << "eABF/DRR: Warning! The maximum scaling factor larger than 1.0 is not recommended!" << '\n';
511 : }
512 10 : log << "eABF/DRR: The kbt(real system) is " << kbt << '\n';
513 10 : dt = getTimeStep();
514 10 : vector<double> ekbt(ndims, 0.0);
515 10 : if (etemp.size() != ndims) {
516 30 : etemp.assign(ndims, kbt / plumed.getAtoms().getKBoltzmann());
517 : }
518 10 : if (tau.size() != ndims) {
519 0 : tau.assign(ndims, 0.5);
520 : }
521 10 : if (friction.size() != ndims) {
522 0 : friction.assign(ndims, 8.0);
523 : }
524 26 : for (size_t i = 0; i < ndims; ++i) {
525 48 : ekbt[i] = etemp[i] * plumed.getAtoms().getKBoltzmann();
526 16 : log << "eABF/DRR: The kbt(extended system) of [" << i << "] is " << ekbt[i]
527 32 : << '\n';
528 32 : log << "eABF/DRR: relaxation time tau [" << i << "] is " << tau[i] << '\n';
529 32 : log << "eABF/DRR: Extended variable [" << i << "] has friction: " << friction[i] << '\n';
530 : }
531 :
532 : // Set up the force grid
533 10 : vector<DRRAxis> zdelim(ndims);
534 26 : for (size_t i = 0; i < ndims; ++i) {
535 16 : log << "eABF/DRR: The " << i << " dimensional grid minimum is " << gmin[i]
536 32 : << '\n';
537 16 : log << "eABF/DRR: The " << i << " dimensional grid maximum is " << gmax[i]
538 32 : << '\n';
539 16 : log << "eABF/DRR: The " << i << " dimensional grid has " << gbin[i]
540 32 : << " bins" << '\n';
541 16 : log << "eABF/DRR: The " << i << " dimensional zgrid minimum is " << zgmin[i]
542 32 : << '\n';
543 16 : log << "eABF/DRR: The " << i << " dimensional zgrid maximum is " << zgmax[i]
544 32 : << '\n';
545 16 : log << "eABF/DRR: The " << i << " dimensional zgrid has " << zgbin[i]
546 32 : << " bins" << '\n';
547 : double l, h;
548 32 : PLMD::Tools::convert(gmin[i], l);
549 32 : PLMD::Tools::convert(gmax[i], h);
550 32 : delim[i].set(l, h, gbin[i]);
551 : double zl,zh;
552 32 : PLMD::Tools::convert(zgmin[i], zl);
553 32 : PLMD::Tools::convert(zgmax[i], zh);
554 32 : zdelim[i].set(zl, zh, zgbin[i]);
555 : }
556 10 : if (kappa.size() != ndims) {
557 8 : kappa.resize(ndims, 0.0);
558 22 : for (size_t i = 0; i < ndims; ++i) {
559 14 : if (kappa[i] <= 0) {
560 14 : log << "eABF/DRR: The spring force constant kappa[" << i
561 14 : << "] is not set." << '\n';
562 56 : kappa[i] = ekbt[i] / (delim[i].getWidth() * delim[i].getWidth());
563 14 : log << "eABF/DRR: set kappa[" << i
564 14 : << "] according to bin width(ekbt/(binWidth^2))." << '\n';
565 : }
566 14 : log << "eABF/DRR: The spring force constant kappa[" << i << "] is "
567 28 : << std::fixed << std::setprecision(10) << kappa[i] << '\n';
568 : }
569 : } else {
570 2 : log << "eABF/DRR: The kappa have been set manually." << '\n';
571 4 : for (size_t i = 0; i < ndims; ++i) {
572 2 : log << "eABF/DRR: The spring force constant kappa[" << i << "] is "
573 4 : << std::fixed << std::setprecision(10) << kappa[i] << '\n';
574 : }
575 : }
576 :
577 26 : for (size_t i = 0; i < ndims; ++i) {
578 48 : mass[i] = kappa[i] * tau[i] * tau[i] / (4 * pi * pi);
579 32 : log << "eABF/DRR: Fictitious mass[" << i << "] is " << mass[i] << '\n';
580 48 : c1[i] = exp(-0.5 * friction[i] * dt);
581 80 : c2[i] = sqrt(ekbt[i] * (1.0 - c1[i] * c1[i]) / mass[i]);
582 : }
583 :
584 42 : for (size_t i = 0; i < ndims; ++i) {
585 : // Position output
586 16 : string comp = getPntrToArgument(i)->getName() + "_fict";
587 16 : addComponentWithDerivatives(comp);
588 16 : if (getPntrToArgument(i)->isPeriodic()) {
589 : string a, b;
590 : double c, d;
591 14 : getPntrToArgument(i)->getDomain(a, b);
592 14 : getPntrToArgument(i)->getDomain(c, d);
593 14 : componentIsPeriodic(comp, a, b);
594 14 : delim[i].setPeriodicity(c, d);
595 : zdelim[i].setPeriodicity(c, d);
596 : } else
597 2 : componentIsNotPeriodic(comp);
598 16 : fictValue[i] = getPntrToComponent(comp);
599 : // Velocity output
600 32 : comp = getPntrToArgument(i)->getName() + "_vfict";
601 16 : addComponent(comp);
602 16 : componentIsNotPeriodic(comp);
603 16 : vfictValue[i] = getPntrToComponent(comp);
604 : // Bias force from eABF/DRR output
605 32 : comp = getPntrToArgument(i)->getName() + "_biasforce";
606 16 : addComponent(comp);
607 16 : componentIsNotPeriodic(comp);
608 16 : biasforceValue[i] = getPntrToComponent(comp);
609 : // Spring force output, useful for perform egABF and other analysis
610 32 : comp = getPntrToArgument(i)->getName() + "_springforce";
611 16 : addComponent(comp);
612 16 : componentIsNotPeriodic(comp);
613 16 : springforceValue[i] = getPntrToComponent(comp);
614 : // Position output, no pbc-aware
615 32 : comp = getPntrToArgument(i)->getName() + "_fictNoPBC";
616 16 : addComponent(comp);
617 16 : componentIsNotPeriodic(comp);
618 16 : fictNoPBCValue[i] = getPntrToComponent(comp);
619 : }
620 :
621 10 : if (outputprefix.length() == 0) {
622 0 : outputprefix = getLabel();
623 : }
624 : // Support multiple replica
625 10 : string replica_suffix = plumed.getSuffix();
626 10 : if (replica_suffix.empty() == false) {
627 4 : outputprefix = outputprefix + replica_suffix;
628 : }
629 20 : outputname = outputprefix + ".drrstate";
630 20 : cptname = outputprefix + ".cpt.drrstate";
631 :
632 10 : if (!isRestart) {
633 : // If you want to use on-the-fly text output for CZAR and naive estimator,
634 : // you should turn it to true first!
635 27 : ABFGrid = ABF(delim, ".abf", fullsamples, maxFactor, textoutput);
636 : // Just initialize it even useCZARestimator is off.
637 27 : CZARestimator = CZAR(zdelim, ".czar", kbt, textoutput);
638 9 : log << "eABF/DRR: The init function of the grid is finished." << '\n';
639 : } else {
640 : // ABF Parametres are not saved in binary files
641 : // So manully set them up
642 1 : ABFGrid.setParameters(fullsamples, maxFactor);
643 : }
644 10 : if (useCZARestimator) {
645 10 : log << "eABF/DRR: Using corrected z-average restraint estimator of gradients" << '\n';
646 40 : log << " Bibliography " << plumed.cite("Lesage, Lelièvre, Stoltz and Hénin, "
647 10 : "J. Phys. Chem. B 3676, 121 (2017)");
648 30 : log << plumed.cite("Darve and Pohorille, J. Chem. Phys. 9169, 115 (2001)") << '\n';
649 : }
650 10 : if (useUIestimator) {
651 8 : log << "eABF/DRR: Using umbrella integration(Zheng and Yang's) estimator "
652 : "of gradients."
653 8 : << '\n';
654 8 : log << "eABF/DRR: The UI estimator code is contributed by Haohao Fu."
655 8 : << '\n';
656 32 : log << " Bibliography " << plumed.cite(
657 8 : "Fu, Shao, Chipot and Cai, J. Chem. Theory Comput. 3506, 12 (2016)");
658 24 : log << plumed.cite("Zheng and Yang, J. Chem. Theory Comput. 810, 8 (2012)");
659 24 : log << plumed.cite("Darve and Pohorille, J. Chem. Phys. 9169, 115 (2001)") << '\n';
660 16 : vector<double> lowerboundary(zdelim.size(), 0);
661 16 : vector<double> upperboundary(zdelim.size(), 0);
662 16 : vector<double> width(zdelim.size(), 0);
663 36 : for (size_t i = 0; i < zdelim.size(); ++i) {
664 14 : lowerboundary[i] = zdelim[i].getMin();
665 14 : upperboundary[i] = zdelim[i].getMax();
666 14 : width[i] = zdelim[i].getWidth();
667 : }
668 8 : vector<string> input_filename;
669 : bool uirestart = false;
670 9 : if (isRestart && (uirprefix.length() != 0)) {
671 1 : input_filename.push_back(uirprefix);
672 : uirestart = true;
673 : }
674 9 : if (isRestart && (uirprefix.length() == 0)) {
675 0 : input_filename.push_back(outputprefix);
676 : }
677 16 : eabf_UI = UIestimator::UIestimator(
678 8 : lowerboundary, upperboundary, width, kappa, outputprefix, int(outputfreq),
679 16 : uirestart, input_filename, kbt / plumed.getAtoms().getKBoltzmann());
680 : }
681 10 : }
682 :
683 106 : void DynamicReferenceRestraining::calculate() {
684 106 : long long int step_now = getStep();
685 106 : if (firsttime) {
686 37 : for (size_t i = 0; i < ndims; ++i) {
687 14 : fict[i] = getArgument(i);
688 : }
689 9 : firsttime = false;
690 : }
691 106 : if (step_now != 0) {
692 96 : if ((step_now % int(outputfreq)) == 0) {
693 12 : save(outputname, step_now);
694 12 : if (textoutput) {
695 7 : ABFGrid.writeAll(outputprefix);
696 7 : if (useCZARestimator) {
697 7 : CZARestimator.writeAll(outputprefix);
698 : }
699 : }
700 : }
701 96 : if (historyfreq > 0 && (step_now % int(historyfreq)) == 0) {
702 : const string filename =
703 0 : outputprefix + "." + std::to_string(step_now) + ".drrstate";
704 0 : save(filename, step_now);
705 0 : if (textoutput) {
706 : const string textfilename =
707 0 : outputprefix + "." + std::to_string(step_now);
708 0 : ABFGrid.writeAll(textfilename);
709 0 : if (useCZARestimator) {
710 0 : CZARestimator.writeAll(textfilename);
711 : }
712 : }
713 : }
714 96 : if (getCPT()) {
715 0 : log << "eABF/DRR: The MD engine is writing checkpoint so we also write a "
716 : "DRR state file at step: "
717 0 : << step_now << ".\n";
718 0 : save(cptname, step_now);
719 : }
720 : }
721 106 : if (withExternalForce == false) {
722 : double ene = 0.0;
723 402 : for (size_t i = 0; i < ndims; ++i) {
724 154 : real[i] = getArgument(i);
725 462 : springlength[i] = difference(i, fict[i], real[i]);
726 462 : fictNoPBC[i] = real[i] - springlength[i];
727 308 : double f = -kappa[i] * springlength[i];
728 154 : ffict_measured[i] = -f;
729 308 : ene += 0.5 * kappa[i] * springlength[i] * springlength[i];
730 : setOutputForce(i, f);
731 154 : ffict[i] = -f;
732 462 : fict[i] = fictValue[i]->bringBackInPbc(fict[i]);
733 308 : fictValue[i]->set(fict[i]);
734 308 : vfictValue[i]->set(vfict_laststep[i]);
735 308 : springforceValue[i]->set(ffict_measured[i]);
736 308 : fictNoPBCValue[i]->set(fictNoPBC[i]);
737 : }
738 : setBias(ene);
739 94 : ABFGrid.store_getbias(fict, ffict_measured, fbias);
740 : } else {
741 60 : for (size_t i = 0; i < ndims; ++i) {
742 24 : real[i] = getArgument(i);
743 48 : ffict_measured[i] = externalForceValue[i]->get();
744 24 : if (withExternalFict) {
745 48 : fictNoPBC[i] = externalFictValue[i]->get();
746 : }
747 48 : springforceValue[i]->set(ffict_measured[i]);
748 48 : fictNoPBCValue[i]->set(fictNoPBC[i]);
749 : }
750 12 : ABFGrid.store_getbias(real, ffict_measured, fbias);
751 12 : if (!nobias) {
752 0 : for (size_t i = 0; i < ndims; ++i) {
753 0 : setOutputForce(i, fbias[i]);
754 : }
755 : }
756 : }
757 106 : if (useCZARestimator) {
758 106 : CZARestimator.store(real, ffict_measured);
759 : }
760 106 : if (useUIestimator) {
761 82 : eabf_UI.update_output_filename(outputprefix);
762 246 : eabf_UI.update(int(step_now), real, fictNoPBC);
763 : }
764 106 : }
765 :
766 106 : void DynamicReferenceRestraining::update() {
767 106 : if (withExternalForce == false) {
768 402 : for (size_t i = 0; i < ndims; ++i) {
769 : // consider additional forces on the fictitious particle
770 : // (e.g. MetaD stuff)
771 308 : ffict[i] += fictValue[i]->getForce();
772 154 : if (!nobias) {
773 308 : ffict[i] += fbias[i];
774 : }
775 308 : biasforceValue[i]->set(fbias[i]);
776 : // update velocity (half step)
777 462 : vfict[i] += ffict[i] * 0.5 * dt / mass[i];
778 : // thermostat (half step)
779 616 : vfict[i] = c1[i] * vfict[i] + c2[i] * rand.Gaussian();
780 : // save full step velocity to be dumped at next step
781 154 : vfict_laststep[i] = vfict[i];
782 : // thermostat (half step)
783 616 : vfict[i] = c1[i] * vfict[i] + c2[i] * rand.Gaussian();
784 : // update velocity (half step)
785 462 : vfict[i] += ffict[i] * 0.5 * dt / mass[i];
786 : // update position (full step)
787 308 : fict[i] += vfict[i] * dt;
788 : }
789 : }
790 106 : }
791 :
792 12 : void DynamicReferenceRestraining::save(const string &filename,
793 : long long int step) {
794 24 : std::ofstream out;
795 12 : out.open(filename.c_str(), std::ios::binary);
796 : boost::archive::binary_oarchive oa(out);
797 60 : oa << step << fict << vfict << vfict_laststep << ffict << ABFGrid
798 12 : << CZARestimator;
799 12 : out.close();
800 12 : }
801 :
802 1 : void DynamicReferenceRestraining::load(const string &rfile_prefix) {
803 1 : string replica_suffix = plumed.getSuffix();
804 : string filename;
805 1 : if (replica_suffix.empty() == true) {
806 2 : filename = rfile_prefix + ".drrstate";
807 : } else {
808 0 : filename = rfile_prefix + "." + replica_suffix + ".drrstate";
809 : }
810 2 : std::ifstream in;
811 : long long int step;
812 1 : in.open(filename.c_str(), std::ios::binary);
813 1 : log << "eABF/DRR: Read restart file: " << filename << '\n';
814 : boost::archive::binary_iarchive ia(in);
815 5 : ia >> step >> fict >> vfict >> vfict_laststep >> ffict >> ABFGrid >>
816 1 : CZARestimator;
817 1 : in.close();
818 1 : log << "eABF/DRR: Restart at step: " << step << '\n';
819 1 : backupFile(filename);
820 1 : }
821 :
822 1 : void DynamicReferenceRestraining::backupFile(const string &filename) {
823 : bool isSuccess = false;
824 : long int i = 0;
825 3 : while (!isSuccess) {
826 : // If libstdc++ support C++17 we can simplify following code.
827 4 : const string bckname = "bck." + filename + "." + std::to_string(i);
828 1 : if (is_file_exist(bckname.c_str())) {
829 0 : ++i;
830 : } else {
831 1 : log << "eABF/DRR: Backup original restart file to " << bckname << '\n';
832 2 : std::ifstream src(filename.c_str(), std::ios::binary);
833 2 : std::ofstream dst(bckname.c_str(), std::ios::binary);
834 1 : dst << src.rdbuf();
835 1 : src.close();
836 1 : dst.close();
837 : isSuccess = true;
838 : }
839 : }
840 1 : }
841 :
842 : // Copy from
843 : // stackoverflow(https://stackoverflow.com/questions/12774207/fastest-way-to-check-if-a-file-exist-using-standard-c-c11-c)
844 1 : bool DynamicReferenceRestraining::is_file_exist(const char *fileName) {
845 2 : std::ifstream infile(fileName);
846 1 : return infile.good();
847 : }
848 : }
849 5517 : }
850 :
851 : #endif
|
