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