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