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