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