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 : 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 mergeHistoryFiles;
221 : bool textoutput;
222 : bool withExternalForce;
223 : bool withExternalFict;
224 : vector<unsigned> reflectingWall;
225 : ABF ABFGrid;
226 : CZAR CZARestimator;
227 : double fullsamples;
228 : vector<double> maxFactors;
229 : UIestimator::UIestimator eabf_UI;
230 : Random rand;
231 :
232 : public:
233 : explicit DynamicReferenceRestraining(const ActionOptions &);
234 : void calculate();
235 : void update();
236 : void save(const string &filename, long long int step);
237 : void load(const string &filename);
238 : void backupFile(const string &filename);
239 : static void registerKeywords(Keywords &keys);
240 : bool is_file_exist(const char *fileName);
241 : };
242 :
243 10443 : PLUMED_REGISTER_ACTION(DynamicReferenceRestraining, "DRR")
244 :
245 13 : void DynamicReferenceRestraining::registerKeywords(Keywords &keys) {
246 13 : Bias::registerKeywords(keys);
247 13 : keys.use("ARG");
248 26 : keys.add("optional", "KAPPA", "specifies that the restraint is harmonic and "
249 : "what the values of the force constants on "
250 : "each of the variables are (default to "
251 : "\\f$k_BT\\f$/(GRID_SPACING)^2)");
252 26 : keys.add("compulsory", "TAU", "0.5", "specifies relaxation time on each of "
253 : "variables are, similar to "
254 : "extended Time Constant in Colvars");
255 26 : keys.add("compulsory", "FRICTION", "8.0",
256 : "add a friction to the variable, similar to extended Langevin Damping "
257 : "in Colvars");
258 26 : keys.add("compulsory", "GRID_MIN", "the lower bounds for the grid (GRID_BIN "
259 : "or GRID_SPACING should be specified)");
260 26 : keys.add("compulsory", "GRID_MAX", "the upper bounds for the grid (GRID_BIN "
261 : "or GRID_SPACING should be specified)");
262 26 : keys.add("compulsory", "REFLECTINGWALL", "0", "whether add reflecting walls "
263 : "for each CV at GRID_MIN and GRID_MAX. Setting non-zero values will "
264 : "enable this feature");
265 26 : keys.add("optional", "GRID_BIN", "the number of bins for the grid");
266 26 : keys.add("optional", "GRID_SPACING", "the approximate grid spacing (to be "
267 : "used as an alternative or together "
268 : "with GRID_BIN)");
269 26 : keys.add("optional", "ZGRID_MIN", "the lower bounds for the grid (ZGRID_BIN"
270 : " or ZGRID_SPACING should be specified)");
271 26 : keys.add("optional", "ZGRID_MAX", "the upper bounds for the grid (ZGRID_BIN"
272 : " or ZGRID_SPACING should be specified)");
273 26 : keys.add("optional", "ZGRID_BIN", "the number of bins for the grid");
274 26 : keys.add("optional", "ZGRID_SPACING", "the approximate grid spacing (to be "
275 : "used as an alternative or together "
276 : "with ZGRID_BIN)");
277 26 : keys.add("optional", "EXTERNAL_FORCE", "use forces from other action instead"
278 : " of internal spring force, this disable the extended system!");
279 26 : keys.add("optional", "EXTERNAL_FICT", "position of external fictitious "
280 : "particles, useful for UIESTIMATOR");
281 26 : keys.add("compulsory", "FULLSAMPLES", "500",
282 : "number of samples in a bin prior to application of the ABF");
283 26 : keys.add("compulsory", "MAXFACTOR", "1.0",
284 : "maximum scaling factor of biasing force");
285 26 : keys.add("compulsory", "OUTPUTFREQ", "write results to a file every N steps");
286 26 : keys.add("optional", "HISTORYFREQ", "save history to a file every N steps");
287 26 : keys.addFlag("NOCZAR", false, "disable the CZAR estimator");
288 26 : keys.addFlag("UI", false,
289 : "enable the umbrella integration estimator");
290 26 : keys.add("optional", "UIRESTARTPREFIX",
291 : "specify the restart files for umbrella integration");
292 26 : keys.add("optional", "OUTPUTPREFIX",
293 : "specify the output prefix (default to the label name)");
294 26 : keys.add("optional", "TEMP", "the system temperature - needed when FRICTION "
295 : "is present. If not provided will be taken from "
296 : "MD code (if available)");
297 26 : keys.add(
298 : "optional", "EXTTEMP",
299 : "the temperature of extended variables (default to system temperature)");
300 26 : keys.add("optional", "DRR_RFILE",
301 : "specifies the restart file (.drrstate file)");
302 26 : keys.addFlag("NOBIAS", false, "DO NOT apply bias forces.");
303 26 : keys.addFlag("TEXTOUTPUT", false, "use text output for grad and count files "
304 : "instead of boost::serialization binary "
305 : "output");
306 26 : keys.addFlag("MERGEHISTORYFILES", false, "output all historic results "
307 : "to a single file rather than multiple .drrstate files. "
308 : "This option is effective only when textOutput is on.");
309 13 : componentsAreNotOptional(keys);
310 26 : keys.addOutputComponent(
311 : "_fict", "default",
312 : "one or multiple instances of this quantity can be referenced "
313 : "elsewhere in the input file. "
314 : "These quantities will named with the arguments of the bias followed by "
315 : "the character string _tilde. It is possible to add forces on these "
316 : "variable.");
317 26 : keys.addOutputComponent(
318 : "_vfict", "default",
319 : "one or multiple instances of this quantity can be referenced "
320 : "elsewhere in the input file. "
321 : "These quantities will named with the arguments of the bias followed by "
322 : "the character string _tilde. It is NOT possible to add forces on these "
323 : "variable.");
324 26 : keys.addOutputComponent(
325 : "_biasforce", "default",
326 : "The bias force from eABF/DRR of the fictitious particle.");
327 26 : keys.addOutputComponent("_springforce", "default", "Spring force between real CVs and extended CVs");
328 26 : keys.addOutputComponent("_fictNoPBC", "default",
329 : "the positions of fictitious particles (without PBC).");
330 13 : }
331 :
332 12 : DynamicReferenceRestraining::DynamicReferenceRestraining(
333 12 : const ActionOptions &ao)
334 12 : : PLUMED_BIAS_INIT(ao), firsttime(true), nobias(false),
335 12 : fictNoPBC(getNumberOfArguments(), 0.0), real(getNumberOfArguments(), 0.0),
336 12 : springlength(getNumberOfArguments(), 0.0),
337 12 : fict(getNumberOfArguments(), 0.0), vfict(getNumberOfArguments(), 0.0),
338 12 : vfict_laststep(getNumberOfArguments(), 0.0),
339 12 : ffict(getNumberOfArguments(), 0.0), fbias(getNumberOfArguments(), 0.0),
340 12 : kappa(getNumberOfArguments(), 0.0), tau(getNumberOfArguments(), 0.0),
341 12 : friction(getNumberOfArguments(), 0.0), etemp(getNumberOfArguments(), 0.0),
342 12 : ffict_measured(getNumberOfArguments(), 0.0),
343 12 : biasforceValue(getNumberOfArguments(), NULL),
344 12 : springforceValue(getNumberOfArguments(), NULL),
345 12 : fictValue(getNumberOfArguments(), NULL),
346 12 : vfictValue(getNumberOfArguments(), NULL),
347 12 : fictNoPBCValue(getNumberOfArguments(), NULL),
348 12 : externalForceValue(getNumberOfArguments(), NULL),
349 12 : externalFictValue(getNumberOfArguments(), NULL),
350 12 : c1(getNumberOfArguments(), 0.0),
351 12 : c2(getNumberOfArguments(), 0.0), mass(getNumberOfArguments(), 0.0),
352 12 : delim(getNumberOfArguments()), outputname(""), cptname(""),
353 12 : outputprefix(""), ndims(getNumberOfArguments()), dt(0.0), kbt(0.0),
354 12 : outputfreq(0.0), historyfreq(-1.0), isRestart(false),
355 12 : useCZARestimator(true), useUIestimator(false), textoutput(false),
356 12 : withExternalForce(false), withExternalFict(false), mergeHistoryFiles(false),
357 12 : reflectingWall(getNumberOfArguments(), 0),
358 36 : maxFactors(getNumberOfArguments(), 1.0)
359 : {
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 : double temp = -1.0;
400 12 : parse("TEMP", temp);
401 12 : parse("FULLSAMPLES", fullsamples);
402 12 : parseVector("MAXFACTOR", maxFactors);
403 12 : parse("OUTPUTFREQ", outputfreq);
404 12 : parse("HISTORYFREQ", historyfreq);
405 24 : parse("OUTPUTPREFIX", outputprefix);
406 : string restart_prefix;
407 24 : parse("DRR_RFILE", restart_prefix);
408 : string uirprefix;
409 12 : parse("UIRESTARTPREFIX", uirprefix);
410 12 : parseArgumentList("EXTERNAL_FORCE", externalForceValue);
411 24 : parseArgumentList("EXTERNAL_FICT", externalFictValue);
412 12 : if (externalForceValue.empty()) {
413 11 : withExternalForce = false;
414 1 : } else if (externalForceValue.size() != ndims) {
415 0 : error("eABF/DRR: Number of forces doesn't match ARGS!");
416 : } else {
417 1 : withExternalForce = true;
418 : }
419 12 : if (withExternalForce && useUIestimator) {
420 1 : if (externalFictValue.empty()) {
421 0 : error("eABF/DRR: No external fictitious particles specified. UI estimator needs it.");
422 1 : } else if(externalFictValue.size() != ndims) {
423 0 : error("eABF/DRR: Number of fictitious particles doesn't match ARGS!");
424 : } else {
425 1 : withExternalFict = true;
426 : }
427 : }
428 12 : if (temp >= 0.0)
429 12 : kbt = plumed.getAtoms().getKBoltzmann() * temp;
430 : else {
431 0 : kbt = plumed.getAtoms().getKbT();
432 0 : if (kbt <= std::numeric_limits<double>::epsilon()) {
433 0 : error("eABF/DRR: It seems the MD engine does not setup the temperature correctly for PLUMED."
434 : "Please set it by the TEMP keyword manually.");
435 : }
436 : }
437 12 : if (fullsamples < 0.5) {
438 0 : fullsamples = 500.0;
439 0 : log << "eABF/DRR: The fullsamples parametre is not set. Set it to "
440 : "500(default)."
441 0 : << '\n';
442 : }
443 12 : if (getRestart()) {
444 1 : if (restart_prefix.length() != 0) {
445 1 : isRestart = true;
446 1 : firsttime = false;
447 1 : load(restart_prefix);
448 : } else {
449 0 : log << "eABF/DRR: You don't specify the file for restarting." << '\n';
450 0 : log << "eABF/DRR: So I assume you are splitting windows." << '\n';
451 0 : isRestart = false;
452 0 : firsttime = true;
453 : }
454 : }
455 :
456 12 : vector<string> gmin(ndims);
457 12 : vector<string> zgmin(ndims);
458 12 : parseVector("GRID_MIN", gmin);
459 24 : parseVector("ZGRID_MIN", zgmin);
460 12 : if (gmin.size() != ndims)
461 0 : error("eABF/DRR: not enough values for GRID_MIN");
462 12 : if (zgmin.size() != ndims) {
463 33 : log << "eABF/DRR: You didn't specify ZGRID_MIN. " << '\n'
464 11 : << "eABF/DRR: The GRID_MIN will be used instead.";
465 11 : zgmin = gmin;
466 : }
467 12 : vector<string> gmax(ndims);
468 12 : vector<string> zgmax(ndims);
469 12 : parseVector("GRID_MAX", gmax);
470 24 : parseVector("ZGRID_MAX", zgmax);
471 12 : if (gmax.size() != ndims)
472 0 : error("eABF/DRR: not enough values for GRID_MAX");
473 12 : if (zgmax.size() != ndims) {
474 33 : log << "eABF/DRR: You didn't specify ZGRID_MAX. " << '\n'
475 11 : << "eABF/DRR: The GRID_MAX will be used instead.";
476 11 : zgmax = gmax;
477 : }
478 12 : vector<unsigned> gbin(ndims);
479 12 : vector<unsigned> zgbin(ndims);
480 12 : vector<double> gspacing(ndims);
481 12 : vector<double> zgspacing(ndims);
482 12 : parseVector("GRID_BIN", gbin);
483 12 : parseVector("ZGRID_BIN", zgbin);
484 12 : parseVector("GRID_SPACING", gspacing);
485 24 : parseVector("ZGRID_SPACING", zgspacing);
486 12 : if (gbin.size() != ndims) {
487 3 : log << "eABF/DRR: You didn't specify GRID_BIN. Trying to use GRID_SPACING "
488 : "instead."
489 3 : << '\n';
490 3 : if (gspacing.size() != ndims) {
491 0 : error("eABF/DRR: not enough values for GRID_BIN");
492 : } else {
493 3 : gbin.resize(ndims);
494 6 : for (size_t i = 0; i < ndims; ++i) {
495 : double l, h;
496 3 : PLMD::Tools::convert(gmin[i], l);
497 3 : PLMD::Tools::convert(gmax[i], h);
498 3 : gbin[i] = std::nearbyint((h - l) / gspacing[i]);
499 3 : gspacing[i] = (h - l) / gbin[i];
500 3 : log << "GRID_BIN[" << i << "] is " << gbin[i] << '\n';
501 : }
502 : }
503 : }
504 12 : if (zgbin.size() != ndims) {
505 11 : log << "eABF/DRR: You didn't specify ZGRID_BIN. Trying to use ZGRID_SPACING instead." << '\n';
506 11 : if (zgspacing.size() != ndims) {
507 11 : log << "eABF/DRR: You didn't specify ZGRID_SPACING. Trying to use GRID_SPACING or GRID_BIN instead." << '\n';
508 11 : zgbin = gbin;
509 11 : zgspacing = gspacing;
510 : } else {
511 0 : zgbin.resize(ndims);
512 0 : for (size_t i = 0; i < ndims; ++i) {
513 : double l, h;
514 0 : PLMD::Tools::convert(zgmin[i], l);
515 0 : PLMD::Tools::convert(zgmax[i], h);
516 0 : zgbin[i] = std::nearbyint((h - l) / zgspacing[i]);
517 0 : zgspacing[i] = (h - l) / zgbin[i];
518 0 : log << "ZGRID_BIN[" << i << "] is " << zgbin[i] << '\n';
519 : }
520 : }
521 : }
522 12 : checkRead();
523 :
524 : // Set up kbt for extended system
525 12 : log << "eABF/DRR: The fullsamples is " << fullsamples << '\n';
526 12 : log << "eABF/DRR: The kbt(real system) is " << kbt << '\n';
527 12 : dt = getTimeStep() * getStride();
528 12 : log << "eABF/DRR: timestep = " << getTimeStep() << " ps with stride = " << getStride() << " steps\n";
529 12 : vector<double> ekbt(ndims, 0.0);
530 12 : if (etemp.size() != ndims) {
531 12 : etemp.assign(ndims, kbt / plumed.getAtoms().getKBoltzmann());
532 : }
533 12 : if (tau.size() != ndims) {
534 0 : tau.assign(ndims, 0.5);
535 : }
536 12 : if (friction.size() != ndims) {
537 0 : friction.assign(ndims, 8.0);
538 : }
539 12 : if (maxFactors.size() != ndims) {
540 0 : maxFactors.assign(ndims, 1.0);
541 : }
542 31 : for (size_t i = 0; i < ndims; ++i) {
543 19 : log << "eABF/DRR: The maximum scaling factor [" << i << "] is " << maxFactors[i] << '\n';
544 19 : if (maxFactors[i] > 1.0) {
545 0 : log << "eABF/DRR: Warning! The maximum scaling factor larger than 1.0 is not recommended!" << '\n';
546 : }
547 : }
548 31 : for (size_t i = 0; i < ndims; ++i) {
549 19 : ekbt[i] = etemp[i] * plumed.getAtoms().getKBoltzmann();
550 19 : log << "eABF/DRR: The kbt(extended system) of [" << i << "] is " << ekbt[i]
551 19 : << '\n';
552 19 : log << "eABF/DRR: relaxation time tau [" << i << "] is " << tau[i] << '\n';
553 19 : log << "eABF/DRR: Extended variable [" << i << "] has friction: " << friction[i] << '\n';
554 : }
555 :
556 : // Set up the force grid
557 12 : vector<DRRAxis> zdelim(ndims);
558 31 : for (size_t i = 0; i < ndims; ++i) {
559 19 : log << "eABF/DRR: The " << i << " dimensional grid minimum is " << gmin[i]
560 19 : << '\n';
561 19 : log << "eABF/DRR: The " << i << " dimensional grid maximum is " << gmax[i]
562 19 : << '\n';
563 19 : log << "eABF/DRR: The " << i << " dimensional grid has " << gbin[i]
564 19 : << " bins" << '\n';
565 19 : log << "eABF/DRR: The " << i << " dimensional zgrid minimum is " << zgmin[i]
566 19 : << '\n';
567 19 : log << "eABF/DRR: The " << i << " dimensional zgrid maximum is " << zgmax[i]
568 19 : << '\n';
569 19 : log << "eABF/DRR: The " << i << " dimensional zgrid has " << zgbin[i]
570 19 : << " bins" << '\n';
571 : double l, h;
572 19 : PLMD::Tools::convert(gmin[i], l);
573 19 : PLMD::Tools::convert(gmax[i], h);
574 19 : delim[i].set(l, h, gbin[i]);
575 : double zl,zh;
576 19 : PLMD::Tools::convert(zgmin[i], zl);
577 19 : PLMD::Tools::convert(zgmax[i], zh);
578 19 : zdelim[i].set(zl, zh, zgbin[i]);
579 : }
580 12 : if (kappa.size() != ndims) {
581 9 : kappa.resize(ndims, 0.0);
582 25 : for (size_t i = 0; i < ndims; ++i) {
583 16 : if (kappa[i] <= 0) {
584 16 : log << "eABF/DRR: The spring force constant kappa[" << i
585 16 : << "] is not set." << '\n';
586 16 : kappa[i] = ekbt[i] / (delim[i].getWidth() * delim[i].getWidth());
587 16 : log << "eABF/DRR: set kappa[" << i
588 16 : << "] according to bin width(ekbt/(binWidth^2))." << '\n';
589 : }
590 16 : log << "eABF/DRR: The spring force constant kappa[" << i << "] is "
591 16 : << std::fixed << std::setprecision(10) << kappa[i] << '\n';
592 : }
593 : } else {
594 3 : log << "eABF/DRR: The kappa have been set manually." << '\n';
595 6 : for (size_t i = 0; i < ndims; ++i) {
596 3 : log << "eABF/DRR: The spring force constant kappa[" << i << "] is "
597 3 : << std::fixed << std::setprecision(10) << kappa[i] << '\n';
598 : }
599 : }
600 :
601 31 : for (size_t i = 0; i < ndims; ++i) {
602 19 : mass[i] = kappa[i] * tau[i] * tau[i] / (4 * pi * pi);
603 19 : log << "eABF/DRR: Fictitious mass[" << i << "] is " << mass[i] << '\n';
604 19 : c1[i] = exp(-0.5 * friction[i] * dt);
605 19 : c2[i] = sqrt(ekbt[i] * (1.0 - c1[i] * c1[i]) / mass[i]);
606 : }
607 :
608 31 : for (size_t i = 0; i < ndims; ++i) {
609 : // Position output
610 19 : string comp = getPntrToArgument(i)->getName() + "_fict";
611 19 : addComponentWithDerivatives(comp);
612 19 : if (getPntrToArgument(i)->isPeriodic()) {
613 : string a, b;
614 : double c, d;
615 16 : getPntrToArgument(i)->getDomain(a, b);
616 16 : getPntrToArgument(i)->getDomain(c, d);
617 16 : componentIsPeriodic(comp, a, b);
618 16 : delim[i].setPeriodicity(c, d);
619 : zdelim[i].setPeriodicity(c, d);
620 : } else
621 3 : componentIsNotPeriodic(comp);
622 19 : fictValue[i] = getPntrToComponent(comp);
623 : // Velocity output
624 19 : comp = getPntrToArgument(i)->getName() + "_vfict";
625 19 : addComponent(comp);
626 19 : componentIsNotPeriodic(comp);
627 19 : vfictValue[i] = getPntrToComponent(comp);
628 : // Bias force from eABF/DRR output
629 19 : comp = getPntrToArgument(i)->getName() + "_biasforce";
630 19 : addComponent(comp);
631 19 : componentIsNotPeriodic(comp);
632 19 : biasforceValue[i] = getPntrToComponent(comp);
633 : // Spring force output, useful for perform egABF and other analysis
634 19 : comp = getPntrToArgument(i)->getName() + "_springforce";
635 19 : addComponent(comp);
636 19 : componentIsNotPeriodic(comp);
637 19 : springforceValue[i] = getPntrToComponent(comp);
638 : // Position output, no pbc-aware
639 19 : comp = getPntrToArgument(i)->getName() + "_fictNoPBC";
640 19 : addComponent(comp);
641 19 : componentIsNotPeriodic(comp);
642 19 : fictNoPBCValue[i] = getPntrToComponent(comp);
643 : }
644 :
645 12 : if (outputprefix.length() == 0) {
646 0 : outputprefix = getLabel();
647 : }
648 : // Support multiple replica
649 12 : string replica_suffix = plumed.getSuffix();
650 12 : if (replica_suffix.empty() == false) {
651 4 : outputprefix = outputprefix + replica_suffix;
652 : }
653 12 : outputname = outputprefix + ".drrstate";
654 12 : cptname = outputprefix + ".cpt.drrstate";
655 :
656 12 : if (!isRestart) {
657 : // If you want to use on-the-fly text output for CZAR and naive estimator,
658 : // you should turn it to true first!
659 11 : ABFGrid = ABF(delim, ".abf", fullsamples, maxFactors, textoutput);
660 : // Just initialize it even useCZARestimator is off.
661 22 : CZARestimator = CZAR(zdelim, ".czar", kbt, textoutput);
662 11 : log << "eABF/DRR: The init function of the grid is finished." << '\n';
663 : } else {
664 : // ABF Parametres are not saved in binary files
665 : // So manully set them up
666 1 : ABFGrid.setParameters(fullsamples, maxFactors);
667 : }
668 12 : if (useCZARestimator) {
669 12 : log << "eABF/DRR: Using corrected z-average restraint estimator of gradients" << '\n';
670 24 : log << " Bibliography " << plumed.cite("Lesage, Lelièvre, Stoltz and Hénin, "
671 24 : "J. Phys. Chem. B 3676, 121 (2017)");
672 12 : log << plumed.cite("Darve and Pohorille, J. Chem. Phys. 9169, 115 (2001)") << '\n';
673 : }
674 12 : if (useUIestimator) {
675 9 : log << "eABF/DRR: Using umbrella integration(Zheng and Yang's) estimator "
676 : "of gradients."
677 9 : << '\n';
678 9 : log << "eABF/DRR: The UI estimator code is contributed by Haohao Fu."
679 9 : << '\n';
680 18 : log << " Bibliography " << plumed.cite(
681 18 : "Fu, Shao, Chipot and Cai, J. Chem. Theory Comput. 3506, 12 (2016)");
682 18 : log << plumed.cite("Zheng and Yang, J. Chem. Theory Comput. 810, 8 (2012)");
683 9 : log << plumed.cite("Darve and Pohorille, J. Chem. Phys. 9169, 115 (2001)") << '\n';
684 9 : vector<double> lowerboundary(zdelim.size(), 0);
685 9 : vector<double> upperboundary(zdelim.size(), 0);
686 9 : vector<double> width(zdelim.size(), 0);
687 24 : for (size_t i = 0; i < zdelim.size(); ++i) {
688 15 : lowerboundary[i] = zdelim[i].getMin();
689 15 : upperboundary[i] = zdelim[i].getMax();
690 15 : width[i] = zdelim[i].getWidth();
691 : }
692 : vector<string> input_filename;
693 : bool uirestart = false;
694 9 : if (isRestart && (uirprefix.length() != 0)) {
695 1 : input_filename.push_back(uirprefix);
696 : uirestart = true;
697 : }
698 9 : if (isRestart && (uirprefix.length() == 0)) {
699 0 : input_filename.push_back(outputprefix);
700 : }
701 9 : eabf_UI = UIestimator::UIestimator(
702 9 : lowerboundary, upperboundary, width, kappa, outputprefix, int(outputfreq),
703 18 : uirestart, input_filename, kbt / plumed.getAtoms().getKBoltzmann());
704 9 : }
705 24 : }
706 :
707 123 : void DynamicReferenceRestraining::calculate() {
708 123 : long long int step_now = getStep();
709 123 : if (firsttime) {
710 28 : for (size_t i = 0; i < ndims; ++i) {
711 17 : fict[i] = getArgument(i);
712 17 : if(reflectingWall[i] && (fict[i]>=delim[i].getMax() || fict[i]<=delim[i].getMin())) {
713 0 : error("eABF/DRR: initial position not in the range of [gmin, gmax]!");
714 : }
715 : }
716 11 : firsttime = false;
717 : }
718 123 : if (step_now != 0) {
719 111 : if ((step_now % int(outputfreq)) == 0) {
720 15 : save(outputname, step_now);
721 15 : if (textoutput) {
722 10 : ABFGrid.writeAll(outputprefix);
723 10 : if (useCZARestimator) {
724 10 : CZARestimator.writeAll(outputprefix);
725 10 : CZARestimator.writeZCount(outputprefix);
726 : }
727 : }
728 : }
729 111 : if (historyfreq > 0 && (step_now % int(historyfreq)) == 0) {
730 4 : if (textoutput) {
731 : const string filename =
732 4 : outputprefix + ".another.drrstate";
733 4 : save(filename, step_now);
734 : const string textfilename =
735 4 : mergeHistoryFiles ? (outputprefix + ".hist") : (outputprefix + "." + std::to_string(step_now));
736 4 : ABFGrid.writeAll(textfilename, mergeHistoryFiles);
737 4 : if (useCZARestimator) {
738 4 : CZARestimator.writeAll(textfilename, mergeHistoryFiles);
739 4 : CZARestimator.writeZCount(textfilename, mergeHistoryFiles);
740 : }
741 : } else {
742 : const string filename =
743 0 : outputprefix + "." + std::to_string(step_now) + ".drrstate";
744 0 : save(filename, step_now);
745 : }
746 : }
747 111 : if (getCPT()) {
748 0 : log << "eABF/DRR: The MD engine is writing checkpoint so we also write a "
749 : "DRR state file at step: "
750 0 : << step_now << ".\n";
751 0 : save(cptname, step_now);
752 : }
753 : }
754 123 : if (withExternalForce == false) {
755 : double ene = 0.0;
756 294 : for (size_t i = 0; i < ndims; ++i) {
757 183 : real[i] = getArgument(i);
758 183 : springlength[i] = difference(i, fict[i], real[i]);
759 183 : fictNoPBC[i] = real[i] - springlength[i];
760 183 : double f = -kappa[i] * springlength[i];
761 183 : ffict_measured[i] = -f;
762 183 : ene += 0.5 * kappa[i] * springlength[i] * springlength[i];
763 : setOutputForce(i, f);
764 183 : ffict[i] = -f;
765 183 : fict[i] = fictValue[i]->bringBackInPbc(fict[i]);
766 183 : fictValue[i]->set(fict[i]);
767 183 : vfictValue[i]->set(vfict_laststep[i]);
768 183 : springforceValue[i]->set(ffict_measured[i]);
769 183 : fictNoPBCValue[i]->set(fictNoPBC[i]);
770 : }
771 : setBias(ene);
772 111 : ABFGrid.store_getbias(fict, ffict_measured, fbias);
773 : } else {
774 36 : for (size_t i = 0; i < ndims; ++i) {
775 24 : real[i] = getArgument(i);
776 24 : ffict_measured[i] = externalForceValue[i]->get();
777 24 : if (withExternalFict) {
778 24 : fictNoPBC[i] = externalFictValue[i]->get();
779 : }
780 24 : springforceValue[i]->set(ffict_measured[i]);
781 24 : fictNoPBCValue[i]->set(fictNoPBC[i]);
782 : }
783 12 : ABFGrid.store_getbias(real, ffict_measured, fbias);
784 12 : if (!nobias) {
785 0 : for (size_t i = 0; i < ndims; ++i) {
786 0 : setOutputForce(i, fbias[i]);
787 : }
788 : }
789 : }
790 123 : if (useCZARestimator) {
791 123 : CZARestimator.store(real, ffict_measured);
792 : }
793 123 : if (useUIestimator) {
794 87 : eabf_UI.update_output_filename(outputprefix);
795 261 : eabf_UI.update(int(step_now), real, fictNoPBC);
796 : }
797 123 : }
798 :
799 123 : void DynamicReferenceRestraining::update() {
800 123 : if (withExternalForce == false) {
801 294 : for (size_t i = 0; i < ndims; ++i) {
802 : // consider additional forces on the fictitious particle
803 : // (e.g. MetaD stuff)
804 183 : ffict[i] += fictValue[i]->getForce();
805 183 : if (!nobias) {
806 183 : ffict[i] += fbias[i];
807 : }
808 183 : biasforceValue[i]->set(fbias[i]);
809 : // update velocity (half step)
810 183 : vfict[i] += ffict[i] * 0.5 * dt / mass[i];
811 : // thermostat (half step)
812 183 : vfict[i] = c1[i] * vfict[i] + c2[i] * rand.Gaussian();
813 : // save full step velocity to be dumped at next step
814 183 : vfict_laststep[i] = vfict[i];
815 : // thermostat (half step)
816 183 : vfict[i] = c1[i] * vfict[i] + c2[i] * rand.Gaussian();
817 : // update velocity (half step)
818 183 : vfict[i] += ffict[i] * 0.5 * dt / mass[i];
819 : // update position (full step)
820 183 : fict[i] += vfict[i] * dt;
821 : // reflecting boundary
822 183 : if (reflectingWall[i]) {
823 24 : if (fict[i]<delim[i].getMin()) {
824 1 : fict[i]=delim[i].getMin()+(delim[i].getMin()-fict[i]);
825 1 : vfict[i]*=-1.0;
826 23 : } else if (fict[i]>delim[i].getMax()) {
827 0 : fict[i]=delim[i].getMax()-(fict[i]-delim[i].getMax());
828 0 : vfict[i]*=-1.0;
829 : }
830 : }
831 : }
832 : }
833 123 : }
834 :
835 19 : void DynamicReferenceRestraining::save(const string &filename,
836 : long long int step) {
837 19 : std::ofstream out;
838 19 : out.open(filename.c_str(), std::ios::binary);
839 19 : boost::archive::binary_oarchive oa(out);
840 19 : oa << step << fict << vfict << vfict_laststep << ffict << ABFGrid
841 19 : << CZARestimator;
842 19 : out.close();
843 19 : }
844 :
845 1 : void DynamicReferenceRestraining::load(const string &rfile_prefix) {
846 1 : string replica_suffix = plumed.getSuffix();
847 : string filename;
848 1 : if (replica_suffix.empty() == true) {
849 2 : filename = rfile_prefix + ".drrstate";
850 : } else {
851 0 : filename = rfile_prefix + "." + replica_suffix + ".drrstate";
852 : }
853 1 : std::ifstream in;
854 : long long int step;
855 1 : in.open(filename.c_str(), std::ios::binary);
856 1 : log << "eABF/DRR: Read restart file: " << filename << '\n';
857 1 : boost::archive::binary_iarchive ia(in);
858 1 : ia >> step >> fict >> vfict >> vfict_laststep >> ffict >> ABFGrid >>
859 1 : CZARestimator;
860 1 : in.close();
861 1 : log << "eABF/DRR: Restart at step: " << step << '\n';
862 1 : backupFile(filename);
863 2 : }
864 :
865 1 : void DynamicReferenceRestraining::backupFile(const string &filename) {
866 : bool isSuccess = false;
867 : long int i = 0;
868 2 : while (!isSuccess) {
869 : // If libstdc++ support C++17 we can simplify following code.
870 2 : const string bckname = "bck." + filename + "." + std::to_string(i);
871 1 : if (is_file_exist(bckname.c_str())) {
872 0 : ++i;
873 : } else {
874 1 : log << "eABF/DRR: Backup original restart file to " << bckname << '\n';
875 1 : std::ifstream src(filename.c_str(), std::ios::binary);
876 1 : std::ofstream dst(bckname.c_str(), std::ios::binary);
877 1 : dst << src.rdbuf();
878 1 : src.close();
879 1 : dst.close();
880 : isSuccess = true;
881 1 : }
882 : }
883 1 : }
884 :
885 : // Copy from
886 : // stackoverflow(https://stackoverflow.com/questions/12774207/fastest-way-to-check-if-a-file-exist-using-standard-c-c11-c)
887 1 : bool DynamicReferenceRestraining::is_file_exist(const char *fileName) {
888 1 : std::ifstream infile(fileName);
889 1 : return infile.good();
890 1 : }
891 : }
892 : }
893 :
894 : #endif
|
