Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2019
3 : National Institute of Advanced Industrial Science and Technology (AIST), Japan.
4 : This file contains module for LogMFD method proposed by Tetsuya Morishita(AIST).
5 :
6 : The LogMFD module is free software: you can redistribute it and/or modify
7 : it under the terms of the GNU Lesser General Public License as published by
8 : the Free Software Foundation, either version 3 of the License, or
9 : (at your option) any later version.
10 :
11 : The LogMFD module is distributed in the hope that it will be useful,
12 : but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : GNU Lesser General Public License for more details.
15 :
16 : You should have received a copy of the GNU Lesser General Public License
17 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
18 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
19 :
20 : //+PLUMEDOC BIAS LOGMFD
21 : /*
22 : Used to perform LogMFD, LogPD, and TAMD/d-AFED.
23 :
24 : ## LogMFD
25 :
26 : Consider a physical system of $N_q$ particles, for which the Hamiltonian is given as
27 :
28 : $$
29 : {H_{\rm MD}}\left( {\bf{\Gamma}} \right) = \sum\limits_{j = 1}^{{N_q}} {\frac{{{\bf{p}}_j^2}}{{2{m_j}}}} + \Phi \left( {\bf{q}} \right)
30 : $$
31 :
32 : where ${\bf q}_j$, ${\bf p}_j$ ($\bf{\Gamma}={\bf q},{\bf p}$), and $m_j$ are the position, momentum, and mass of particle $j$ respectively,
33 : and $\Phi$ is the potential energy function for ${\bf q}$.
34 : The free energy $F({\bf X})$ as a function of a set of $N$ collective variables (CVs) is given as
35 :
36 : $$
37 : \begin{aligned}
38 : F\left( {{\bf X}} \right) &= - {k_B}T\log \int {\exp \left[ { - \beta {H_{\rm MD}}} \right]\prod\limits_{i = 1}^N {\delta \left( {{s_i}\left( {{\bf q}} \right) - {X_i}} \right)} d{\bf{\Gamma }}} \\
39 : &\simeq - {k_B}T\log \int {\exp \left[ { - \beta \left\{ {{H_{\rm MD}} + \sum\limits_i^N {\frac{{{K_i}}}{2}{{\left( {{s_i}\left( {{\bf q}} \right) - {X_i}} \right)}^2}} } \right\}} \right]} d{\bf{\Gamma }}
40 : \end{aligned}
41 : $$
42 :
43 : where $\bf{s}$ are the CVs, $k_B$ is Boltzmann constant, $\beta=1/k_BT$,
44 : and $K_i$ is the spring constant which is large enough to invoke
45 :
46 : $$
47 : \delta \left( x \right) = \lim_{k \to \infty } \sqrt {\beta k/2\pi} \exp \left( -\beta kx^2/2 \right)
48 : $$
49 :
50 : In mean-force dynamics, ${\bf X}$ are treated as fictitious dynamical variables, which are associated with the following Hamiltonian,
51 :
52 : $$
53 : {H_{\rm log}} = \sum\limits_{i = 1}^N {\frac{{P_{{X_i}}^2}}{{2{M_i}}} + \Psi \left( {{\bf X}} \right)}
54 : $$
55 :
56 : where ${P_{{X_i}}}$ and $M_i$ are the momentum and mass of $X_i$, respectively, and $\Psi$ is the potential function for ${\bf X}$.
57 : We assume that $\Psi$ is a functional of $F({\bf X})$; $Ψ[F({\bf X})]$. The simplest form of $\Psi$ is $\Psi = F({\bf X})$,
58 : which corresponds to the TAMD/d-AFED methods that are discussed in the first two papers cited below (or the extended Lagrangian dynamics in the limit of the adiabatic decoupling between $\bf{q}$ and ${\bf X}$).
59 : In the case of LogMFD, the following form of $\Psi$ is introduced and discussed in the third paper cited below.
60 :
61 :
62 : $$
63 : {\Psi _{\rm log}}\left( {{\bf X}} \right) = \gamma \log \left[ {\alpha F\left( {{\bf X}} \right) + 1} \right]
64 : $$
65 :
66 : where $\alpha$ (ALPHA) and $\gamma$ (GAMMA) are positive parameters. The logarithmic form of $\Psi_{\rm log}$ ensures the dynamics of ${\bf X}$ on a much smoother energy surface [i.e., $\Psi_{\rm log}({\bf X})$] than $F({\bf X})$, thus enhancing the sampling in the ${\bf X}$-space. The parameters $\alpha$ and $\gamma$ determine the degree of flatness of $\Psi_{\rm log}$, but adjusting only $\alpha$ is normally sufficient to have a relatively flat surface (with keeping the relation $\gamma=1/\alpha$).
67 :
68 : The equation of motion for $X_i$ in LogMFD (no thermostat) is
69 :
70 : $$
71 : {M_i}{\ddot X_i} = - \left( {\frac{{\alpha \gamma }}{{\alpha F + 1}}} \right)\frac{{\partial F}}{{\partial {X_i}}}
72 : $$
73 :
74 : where $-\partial F/\partial X_i$ is evaluated as a canonical average under the condition that ${\bf X}$ is fixed;
75 :
76 : $$
77 : \begin{aligned}
78 : - \frac{{\partial F}}{{\partial {X_i}}} &\simeq \frac{1}{Z}\int {{K_i}\left( {{s_i}\left( {{\bf q}} \right) - {X_i}} \right)\exp \left[ { - \beta \left\{ {{H_{\rm MD}} + \sum\limits_i^N {\frac{{{K_i}}}{2}{{\left( {{s_i}\left( {{\bf q}} \right) - {X_i}} \right)}^2}} } \right\}} \right]} d{\bf{\Gamma }}\\
79 : &\equiv {\left\langle {{K_i}\left( {{s_i}\left( {{\bf q}} \right) - {X_i}} \right)} \right\rangle _{{\bf X}}}
80 : \end{aligned}
81 : $$
82 :
83 : where
84 :
85 : $$
86 : Z = \int {\exp \left[ { - \beta \left\{ {{H_{\rm MD}} + \sum\limits_i^N {\frac{{{K_i}}}{2}{{\left( {{s_i}\left( {{\bf q}} \right) - {X_i}} \right)}^2}} } \right\}} \right]} d{\bf{\Gamma }}
87 : $$
88 :
89 : The mean-force (MF) is practically evaluated by performing a shot-time canonical MD run of $N_m$ steps each time ${\bf X}$ is updated according to the equation of motion for ${\bf X}$.
90 :
91 : If the canonical average for the MF is effectively converged, the dynamical variables ${\bf q}$ and ${\bf X}$ are decoupled and they evolve adiabatically, which can be exploited for the on-the-fly evaluation of $F({\bf X})$. I.e., $H_{\rm log}$ should be a constant of motion in this case, thus $F({\bf X})$ can be evaluated each time ${\bf X}$ is updated as
92 :
93 :
94 : $$
95 : F\left( {{{\bf X}}\left( t \right)} \right) = \frac{1}{\alpha} \left[
96 : \exp \frac{1}{\gamma} \left\{ \left( H_{\rm log} - \sum_i \frac{P_{X_i}^2}{2M_i} \right) \right\} - 1 \right]
97 : $$
98 :
99 :
100 : This means that $F({\bf X})$ can be constructed without post processing (on-the-fly free energy reconstruction). Note that the on-the-fly free energy reconstruction is also possible in TAMD/d-AFED if the Hamiltonian-like conserved quantity is available (e.g., the Nose-Hoover type dynamics).
101 :
102 :
103 : ## LogPD
104 :
105 :
106 : The accuracy in the MF is critical to the on-the-fly free energy reconstruction. To improve the evaluation of the MF, parallel-dynamics (PD) is incorporated into LogMFD, leading to logarithmic parallel-dynamics (LogPD) that is discussed in the fourth paper in the reference section.
107 :
108 :
109 : In PD, the MF is evaluated by a non-equilibrium path-ensemble based on the Crooks-Jarzynski non-equilibrium work relation. To this end, multiple replicas of the MD system which run in parallel are introduced. The CVs [${\bf s}({\bf q})$] in each replica is restrained to the same value of ${\bf X}(t)$. A canonical MD run with $N_m$ steps is performed in each replica, then the MF on $X_i$ is evaluated using the MD trajectories from all replicas.
110 : The MF is practically calculated as
111 :
112 :
113 : $$
114 : - \frac{{\partial F}}{{\partial {X_i}}} = \sum\limits_{k = 1}^{{N_r}} {{W_k}} \sum\limits_{n = 1}^{{N_m}} {\frac{1}{{{N_m}}}{K_i}\left[ {{s_i}\left( {{{\bf q}}_n^k} \right) - {X_i}} \right]}
115 : $$
116 :
117 :
118 :
119 : where ${\bf q}^k_n$ indicates the ${\bf q}$-configuration at time step $n$ in the $N_m$-step shot-time MD run in the $k$th replica among $N_r$ replicas. $W_k$ is given as
120 :
121 : $$
122 : {W_k} = \frac{{{e^{ - \beta {w_k}\left( t \right)}}}}{{\sum\limits_{k=1}^{{N_r}} {{e^{ - \beta {w_k}\left( t \right)}}} }}
123 : $$
124 :
125 :
126 : where
127 :
128 :
129 : $$
130 : {w_k}\left( t \right) = \int\limits_{{X_0}}^{X\left( t \right)} {\sum\limits_{i=1}^N {\frac{{\partial H_{\rm MD}^k}}{{\partial {X_i}}}d{X_i}} }
131 : $$
132 :
133 : $$
134 : H_{\rm MD}^k\left( {{\bf{\Gamma }},{{\bf X}}} \right) = {H_{\rm MD}}\left( {{{\bf{\Gamma }}^k}} \right) + \sum\limits_{i = 1}^N {\frac{{{K_i}}}{2}{{\left( {s_i^k - {X_i}} \right)}^2}}
135 : $$
136 :
137 : and $s^k_i$ is the $i$th CV in the $k$th replica.
138 :
139 : $W_k$ comes from the Crooks-Jarzynski non-equilibrium work relation by which we can evaluate an equilibrium ensemble average from a set of non-equilibrium trajectories. Note that, to avoid possible numerical errors in the exponential function, the following form of $W_k$ is instead used in PLUMED,
140 :
141 : $$
142 : {W_k}\left( t \right) = \frac{{\exp \left[ { - \beta \left\{ {{w_k}\left( t \right) - {w_{\min }}\left( t \right)} \right\}} \right]}}{{\sum\nolimits_k {\exp \left[ { - \beta \left\{ {{w_k}\left( t \right) - {w_{\min }}\left( t \right)} \right\}} \right]} }}
143 : $$
144 :
145 : where
146 :
147 : $$
148 : {w_{\min }}\left( t \right) = {\rm Min}_k\left[ {{w_k}\left( t \right)} \right]
149 : $$
150 :
151 : With the MF evaluated using the PD approach, free energy profiles can be reconstructed more efficiently (requiring less elapsed computing time) in LogPD than with a single MD system in LogMFD. In the case that there exist more than one stable state separated by high energy barriers in the hidden subspace orthogonal to the CV-subspace, LogPD is particularly of use to incorporate all the contributions from such hidden states with appropriate weights (in the limit $N_r\to\infty$ ).
152 :
153 : Note that LogPD calculations should always be initiated with an equilibrium ${\bf q}$-configuration in each replica, because the Crooks-Jarzynski non-equilibrium work relation is invoked. Also note that LogPD is currently available only with Gromacs, while LogMFD can be performed with LAMMPS, Gromacs, Quantum Espresso, NAMD, and so on.
154 :
155 : ## Using LogMFD/PD with a thermostat
156 :
157 : Introducing a thermostat on ${\bf X}$ is often recommended in LogMFD/PD to maintain the adiabatic decoupling between ${\bf q}$ and ${\bf X}$. In the framework of the LogMFD approach, the Nose-Hoover type thermostat and the Gaussian isokinetic (velocity scaling) thermostat can be used to control the kinetic energy of ${\bf X}$.
158 :
159 : ## Nose-Hoover LogMFD/PD
160 :
161 : The equation of motion for $X_i$ coupled to a Nose-Hoover thermostat variable $\eta$ (single heat bath) is
162 :
163 : $$
164 : {M_i}{\ddot X_i} = - \left( {\frac{{\alpha \gamma }}{{\alpha F + 1}}} \right)\frac{{\partial F}}{{\partial {X_i}}} - {M_i}{\dot X_i}\dot \eta
165 : $$
166 :
167 : The equation of motion for $\eta$ is
168 :
169 : $$
170 : Q\ddot \eta = \sum\limits_{i = 1}^N {\frac{{P_{{X_i}}^2}}{{{M_i}}} - N{k_B}T}
171 : $$
172 :
173 : where $N$ is the number of the CVs. Since the following pseudo-Hamiltonian is a constant of motion in Nose-Hoover LogMFD/PD,
174 :
175 : $$
176 : H_{\rm log}^{\rm NH} = \sum\limits_{i = 1}^N {\frac{{P_{{X_i}}^2}}{{2{M_i}}} + \gamma \log \left[ {\alpha F\left( {{\bf X}} \right) + 1} \right] + \frac{1}{2}Q{{\dot \eta }^2} + N{k_B}T\eta}
177 : $$
178 :
179 : $F({\bf X}(t))$ is obtained at each MFD step as
180 :
181 : $$
182 : F\left( {{{\bf X}}\left( t \right)} \right) = \frac{1}{\alpha }\left[ {\exp \left\{ {{{ \frac{1}{\gamma} \left( {H_{\rm log}^{\rm NH} - \sum_i {\frac{{P_{{X_i}}^2}}{{2{M_i}}}} - \frac{1}{2}Q{{\dot \eta }^2} - N{k_B}T\eta} \right)} }} \right\} - 1} \right]
183 : $$
184 :
185 : ## Velocity scaling LogMFD/PD
186 :
187 : The velocity scaling algorithm (which is equivalent to the Gaussian isokinetic dynamics in the limit $\Delta t\to 0$) can also be employed to control the velocity of ${\bf X}$, $\bf{V}_x$.
188 :
189 : The following algorithm is introduced to perform the isokinetic LogMFD calculations that are discussed in the fifth paper cited belo;
190 :
191 : $$
192 : \begin{aligned}
193 : {V_{{X_i}}}\left( {{t_{n + 1}}} \right)
194 : &=
195 : V_{X_i}^\prime \left( {{t_n}} \right) + \Delta t \left[
196 : { - \left( {\frac{{\alpha \gamma }}{{\alpha F\left( {{t_n}} \right) + 1}}} \right)
197 : \frac{{\partial F\left( {{t_n}} \right)}}{{\partial {X_i}}}}
198 : \right]\\
199 : S\left( {{t_{n + 1}}} \right)
200 : &=
201 : \sqrt {\frac{{N{k_B}T}}{{\sum\limits_i {{M_i}V_{{X_i}}^2\left( {{t_{n + 1}}} \right)} }}} \\
202 : {V_{{X_i}}}^\prime \left( {{t_{n + 1}}} \right)
203 : &=
204 : S\left( {{t_{n + 1}}} \right){V_{{X_i}}}\left( {{t_{n + 1}}} \right)\\
205 : {X_i}\left( {{t_{n + 1}}} \right)
206 : &=
207 : {X_i}\left( {{t_n}} \right) + \Delta t V_{X_i}^\prime \left( {{t_{n + 1}}} \right)\\
208 : {\Psi_{\rm log}}\left( {{t_{n + 1}}} \right)
209 : &=
210 : N{k_B}T\log S\left( {{t_{n + 1}}} \right) + {\Psi_{\rm log}}\left( {{t_n}} \right)\\
211 : F\left( {{t_{n + 1}}} \right)
212 : &=
213 : \frac{1}{\alpha} \left[
214 : \exp \left\{ \Psi_{\rm log} \left( t_{n+1} \right) / \gamma \right\} - 1
215 : \right]
216 : \end{aligned}
217 : $$
218 :
219 : Note that $V_{X_i}^\prime\left( {{t_0}} \right)$ is assumed to be initially given, which meets the following relation,
220 :
221 : $$
222 : \sum\limits_{i = 1}^N M_i V_{X_i}^{\prime 2} \left( t_0 \right) = N{k_B}{T}
223 : $$
224 :
225 : It should be stressed that, in the same way as in the NVE and Nose-Hoover LogMFD/PD, $F({\bf X}(t))$ can be evaluated at each MFD step (on-the-fly free energy reconstruction) in Velocity Scaling LogMFD/PD.
226 :
227 :
228 : ## Examples
229 :
230 : ### Example-LoGMFD Example of LogMFD
231 :
232 : The following input file is called plumed.dat and tells plumed to restrain collective variables
233 : to the fictitious dynamical variables in LogMFD/PD.
234 :
235 : ```plumed
236 : UNITS TIME=fs LENGTH=1.0 ENERGY=kcal/mol MASS=1.0 CHARGE=1.0
237 : phi: TORSION ATOMS=5,7,9,15
238 : psi: TORSION ATOMS=7,9,15,17
239 :
240 : # LogMFD
241 : LOGMFD ...
242 : LABEL=logmfd
243 : ARG=phi,psi
244 : KAPPA=1000.0,1000.0
245 : DELTA_T=1.0
246 : INTERVAL=200
247 : TEMP=300.0
248 : FLOG=2.0
249 : MFICT=5000000,5000000
250 : VFICT=3.5e-4,3.5e-4
251 : ALPHA=4.0
252 : THERMOSTAT=NVE
253 : FICT_MAX=3.15,3.15
254 : FICT_MIN=-3.15,-3.15
255 : ... LOGMFD
256 : ```
257 :
258 : To submit this simulation with Gromacs, use the following command line
259 : to execute a LogMFD run with Gromacs-MD.
260 : Here topol.tpr is the input file
261 : which contains atomic coordinates and Gromacs parameters.
262 :
263 : ````
264 : gmx_mpi mdrun -s topol.tpr -plumed
265 : ````
266 :
267 : This command will output files named logmfd.out and replica.out.
268 :
269 : The output file logmfd.out records free energy and all fictitious dynamical variables at every MFD step.
270 :
271 :
272 : ````
273 : # LogMFD
274 : # CVs : phi psi
275 : # Mass for CV particles : 5000000.000000 5000000.000000
276 : # Mass for thermostat : 11923.224809
277 : # 1:iter_mfd, 2:Flog, 3:2*Ekin/gkb[K], 4:eta, 5:Veta,
278 : # 6:phi_fict(t), 7:phi_vfict(t), 8:phi_force(t),
279 : # 9:psi_fict(t), 10:psi_vfict(t), 11:psi_force(t),
280 : 0 2.000000 308.221983 0.000000 0.000000 -2.442363 0.000350 5.522717 2.426650 0.000350 7.443177
281 : 1 1.995466 308.475775 0.000000 0.000000 -2.442013 0.000350 -4.406246 2.427000 0.000350 11.531345
282 : 2 1.992970 308.615664 0.000000 0.000000 -2.441663 0.000350 -3.346513 2.427350 0.000350 15.763196
283 : 3 1.988619 308.859946 0.000000 0.000000 -2.441313 0.000350 6.463092 2.427701 0.000351 6.975422
284 : ...
285 : ````
286 :
287 : The output file replica.out records all collective variables at every MFD step.
288 :
289 : ````
290 : Replica No. 0 of 1.
291 : # 1:iter_mfd, 2:work, 3:weight,
292 : # 4:phi(q)
293 : # 5:psi(q)
294 : 0 0.000000e+00 1.000000e+00 -2.436841 2.434093
295 : 1 -4.539972e-03 1.000000e+00 -2.446420 2.438531
296 : 2 -7.038516e-03 1.000000e+00 -2.445010 2.443114
297 : 3 -1.139672e-02 1.000000e+00 -2.434850 2.434677
298 : ...
299 : ````
300 :
301 : ### Example of LogPD
302 :
303 : Use the following command line to execute a LogPD run using two MD replicas (note that only Gromacs is currently available for LogPD).
304 : Here 0/topol.tpr and 1/topol.tpr are the input files,
305 : each containing the atomic coordinates for the corresponding replica and Gromacs parameters. All the directories, 0/ and 1/, should contain the same plumed.dat.
306 :
307 : ````
308 : mpirun -np 2 gmx_mpi mdrun -s topol -plumed -multidir 0 1
309 : ````
310 :
311 : This command will output files named logmfd.out in 0/, and replica.out.0 and replica.out.1 in 0/ and 1/, respectively.
312 :
313 : The output file logmfd.out records free energy and all fictitious dynamical variables at every MFD step.
314 :
315 : ````
316 : # LogPD, replica parallel of LogMFD
317 : # number of replica : 2
318 : # CVs : phi psi
319 : # Mass for CV particles : 5000000.000000 5000000.000000
320 : # Mass for thermostat : 11923.224809
321 : # 1:iter_mfd, 2:Flog, 3:2*Ekin/gkb[K], 4:eta, 5:Veta,
322 : # 6:phi_fict(t), 7:phi_vfict(t), 8:phi_force(t),
323 : # 9:psi_fict(t), 10:psi_vfict(t), 11:psi_force(t),
324 : 0 2.000000 308.221983 0.000000 0.000000 -2.417903 0.000350 4.930899 2.451475 0.000350 -3.122024
325 : 1 1.999367 308.257404 0.000000 0.000000 -2.417552 0.000350 0.851133 2.451825 0.000350 -1.552718
326 : 2 1.999612 308.243683 0.000000 0.000000 -2.417202 0.000350 -1.588175 2.452175 0.000350 1.601274
327 : 3 1.999608 308.243922 0.000000 0.000000 -2.416852 0.000350 4.267745 2.452525 0.000350 -4.565621
328 : ...
329 : ````
330 :
331 :
332 : The output file replica.out.0 records all collective variables and the Jarzynski weight of replica No.0 at every MFD step.
333 :
334 : ````
335 : # Replica No. 0 of 2.
336 : # 1:iter_mfd, 2:work, 3:weight,
337 : # 4:phi(q)
338 : # 5:psi(q)
339 : 0 0.000000e+00 5.000000e-01 -2.412607 2.446191
340 : 1 -4.649101e-06 4.994723e-01 -2.421403 2.451318
341 : 2 1.520985e-03 4.983996e-01 -2.420269 2.455049
342 : 3 1.588855e-03 4.983392e-01 -2.411081 2.447394
343 : ...
344 : ````
345 :
346 : The output file replica.out.1 records all collective variables and the Jarzynski weight of replica No.1 at every MFD step.
347 :
348 : ````
349 : # Replica No. 1 of 2.
350 : # 1:iter_mfd, 2:work, 3:weight,
351 : # 4:phi(q)
352 : # 5:psi(q)
353 : 0 0.000000e+00 5.000000e-01 -2.413336 2.450516
354 : 1 -1.263077e-03 5.005277e-01 -2.412009 2.449229
355 : 2 -2.295444e-03 5.016004e-01 -2.417322 2.452512
356 : 3 -2.371507e-03 5.016608e-01 -2.414078 2.448521
357 : ...
358 : ````
359 :
360 : */
361 : //+ENDPLUMEDOC
362 :
363 : #include "bias/Bias.h"
364 : #include "core/ActionRegister.h"
365 : #include "tools/Communicator.h"
366 :
367 : #include <iostream>
368 :
369 : using namespace std;
370 : using namespace PLMD;
371 : using namespace bias;
372 :
373 : namespace PLMD {
374 : namespace logmfd {
375 : /**
376 : \brief class for LogMFD parameters, variables and subroutines.
377 : */
378 : class LogMFD : public Bias {
379 : bool firsttime; ///< flag that indicates first MFD step or not.
380 : int step_initial; ///< initial MD step.
381 : int interval; ///< input parameter, period of MD steps when fictitious dynamical variables are evolved.
382 : double delta_t; ///< input parameter, one time step of MFD when fictitious dynamical variables are evolved.
383 : string thermostat; ///< input parameter, type of thermostat for canonical dyanamics.
384 : double kbt; ///< k_B*temperature
385 : double kbtpd; ///< k_B*temperature for PD
386 :
387 : int TAMD; ///< input parameter, perform TAMD instead of LogMFD.
388 : double alpha; ///< input parameter, alpha parameter for LogMFD.
389 : double gamma; ///< input parameter, gamma parameter for LogMFD.
390 : std::vector<double> kappa; ///< input parameter, strength of the harmonic restraining potential.
391 :
392 : std::vector<double> fict_max; ///< input parameter, maximum of each fictitous dynamical variable.
393 : std::vector<double> fict_min; ///< input parameter, minimum of each fictitous dynamical variable.
394 :
395 :
396 : std::vector<double> fict; ///< current values of each fictitous dynamical variable.
397 : std::vector<double> vfict; ///< current velocity of each fictitous dynamical variable.
398 : std::vector<double> mfict; ///< mass of each fictitous dynamical variable.
399 :
400 : double xeta; ///< current eta variable of thermostat.
401 : double veta; ///< current velocity of eta variable of thermostat.
402 : double meta; ///< mass of eta variable of thermostat.
403 :
404 : double phivs; ///< potential used in VS method
405 : double work; ///< current works done by fictitious dynamical variables in this replica.
406 : double weight; ///< current weight of this replica.
407 : double flog; ///< current free energy
408 : double hlog; ///< value invariant
409 :
410 : struct {
411 : std::vector<double> fict;
412 : std::vector<double> vfict;
413 : std::vector<double> ffict;
414 : double xeta;
415 : double veta;
416 : double phivs;
417 : double work;
418 : double weight;
419 : } backup;
420 :
421 : std::vector<double> ffict; ///< current force of each fictitous dynamical variable.
422 : std::vector<double> fict_ave; ///< averaged values of each collective variable.
423 :
424 : std::vector<Value*> fictValue; ///< pointers to fictitious dynamical variables
425 : std::vector<Value*> vfictValue; ///< pointers to velocity of fictitious dynamical variables
426 :
427 : public:
428 : static void registerKeywords(Keywords& keys);
429 :
430 : explicit LogMFD(const ActionOptions&);
431 : void calculate();
432 : void update();
433 : void updateNVE();
434 : void updateNVT();
435 : void updateVS();
436 : void updateWork();
437 : void calcMeanForce();
438 : double calcEkin();
439 : double calcFlog();
440 : double calcClog();
441 :
442 : private:
443 : double sgn( double x ) {
444 185 : return x>0.0 ? 1.0 : x<0.0 ? -1.0 : 0.0;
445 : }
446 : };
447 :
448 : PLUMED_REGISTER_ACTION(LogMFD,"LOGMFD")
449 :
450 : /**
451 : \brief instruction of parameters for Plumed manual.
452 : */
453 17 : void LogMFD::registerKeywords(Keywords& keys) {
454 17 : Bias::registerKeywords(keys);
455 17 : keys.add("compulsory","INTERVAL",
456 : "Period of MD steps (N_m) to update fictitious dynamical variables." );
457 17 : keys.add("compulsory","DELTA_T",
458 : "Time step for the fictitious dynamical variables (DELTA_T=1 often works)." );
459 17 : keys.add("compulsory","THERMOSTAT",
460 : "Type of thermostat for the fictitious dynamical variables. NVE, NVT, VS are available." );
461 17 : keys.add("optional","TEMP",
462 : "Target temperature for the fictitious dynamical variables in LogMFD/PD. "
463 : "It is recommended to set TEMP to be the same as "
464 : "the temperature of the MD system in any thermostated LogMFD/PD run. "
465 : "If not provided, it will be taken from the temperature of the MD system (Gromacs only)." );
466 :
467 17 : keys.add("optional","TAMD",
468 : "When TAMD=1, TAMD/d-AFED calculations can be performed instead of LogMFD. Otherwise, the LogMFD protocol is switched on (default)." );
469 :
470 17 : keys.add("optional","ALPHA",
471 : "Alpha parameter for LogMFD. "
472 : "If not provided or provided as 0, it will be taken as 1/gamma. "
473 : "If gamma is also not provided, Alpha is set as 4, which is a sensible value when the unit of kcal/mol is used." );
474 17 : keys.add("optional","GAMMA",
475 : "Gamma parameter for LogMFD. "
476 : "If not provided or provided as 0, it will be taken as 1/alpha. "
477 : "If alpha is also not provided, Gamma is set as 0.25, which is a sensible value when the unit of kcal/mol is used." );
478 17 : keys.add("compulsory","KAPPA",
479 : "Spring constant of the harmonic restraining potential." );
480 :
481 17 : keys.add("compulsory","FICT_MAX",
482 : "Maximum values reachable for the fictitious dynamical variables. The variables will elastically bounce back at the boundary (mirror boundary)." );
483 17 : keys.add("compulsory","FICT_MIN",
484 : "Minimum values reachable for the fictitious dynamical variables. The variables will elastically bounce back at the boundary (mirror boundary)." );
485 :
486 17 : keys.add("optional","FICT",
487 : "The initial values of the fictitious dynamical variables. "
488 : "If not provided, they are set equal to their corresponding CVs for the initial atomic configuration." );
489 17 : keys.add("optional","VFICT",
490 : "The initial velocities of the fictitious dynamical variables. "
491 : "If not provided, they will be taken as 0. "
492 : "If THERMOSTAT=VS, they are instead automatically adjusted according to TEMP. " );
493 17 : keys.add("optional","MFICT",
494 : "Masses of each fictitious dynamical variable. "
495 : "If not provided, they will be taken as 10000." );
496 :
497 17 : keys.add("optional","XETA",
498 : "The initial eta variable of the Nose-Hoover thermostat "
499 : "for the fictitious dynamical variables. "
500 : "If not provided, it will be taken as 0." );
501 17 : keys.add("optional","VETA",
502 : "The initial velocity of eta variable. "
503 : "If not provided, it will be taken as 0." );
504 17 : keys.add("optional","META",
505 : "Mass of eta variable. "
506 : "If not provided, it will be taken as N*kb*T*100*100." );
507 :
508 17 : keys.add("compulsory","FLOG",
509 : "The initial free energy value in the LogMFD/PD run."
510 : "The origin of the free energy profile is adjusted by FLOG to "
511 : "realize F({X}(t)) > 0 at any X(t), "
512 : "resulting in enhanced barrier-crossing. "
513 : "(The value of H_log is automatically "
514 : "set according to FLOG). "
515 : "In fact, F({X}(t)) is correctly "
516 : "estimated in PLUMED even when F({X}(t)) < 0 in "
517 : "the LogMFD/PD run." );
518 :
519 17 : keys.add("optional","WORK",
520 : "The initial value of the work done by fictitious dynamical "
521 : "variables in each replica. "
522 : "If not provided, it will be taken as 0.");
523 :
524 17 : keys.add("optional","TEMPPD",
525 : "Temperature of the Boltzmann factor in the Jarzynski weight in LogPD (Gromacs only). "
526 : "TEMPPD should be the same as the "
527 : "temperature of the MD system, while TEMP may be (in principle) different from it. "
528 : "If not provided, TEMPPD is set to be the same value as TEMP." );
529 :
530 34 : keys.addOutputComponent("_fict","default","scalar",
531 : "For example, the fictitious collective variable for LogMFD is specified as "
532 : "ARG=dist12 and LABEL=logmfd in LOGMFD section in Plumed input file, "
533 : "the associated fictitious dynamical variable can be specified as "
534 : "PRINT ARG=dist12,logmfd.dist12_fict FILE=COLVAR");
535 34 : keys.addOutputComponent("_vfict","default","scalar",
536 : "For example, the fictitious collective variable for LogMFD is specified as "
537 : "ARG=dist12 and LABEL=logmfd in LOGMFD section in Plumed input file, the "
538 : "velocity of the associated fictitious dynamical variable can be specified as "
539 : "PRINT ARG=dist12,logmfd.dist12_vfict FILE=COLVAR");
540 17 : keys.addDOI("10.1021/jp805039u");
541 17 : keys.addDOI("10.1016/j.cplett.2006.05.062");
542 17 : keys.addDOI("10.1103/PhysRevE.85.066702");
543 17 : keys.addDOI("10.1021/acs.jctc.7b00252");
544 17 : keys.addDOI("10.1016/j.cplett.2018.07.011");
545 17 : }
546 :
547 :
548 : /**
549 : \brief constructor of LogMFD class
550 : \details This constructor initializes all parameters and variables,
551 : reads input file and set value of parameters and initial value of variables,
552 : and writes messages as Plumed log.
553 : */
554 13 : LogMFD::LogMFD( const ActionOptions& ao ):
555 : PLUMED_BIAS_INIT(ao),
556 13 : firsttime(true),
557 13 : step_initial(0),
558 13 : interval(10),
559 13 : delta_t(1.0),
560 26 : thermostat("NVE"),
561 13 : kbt(-1.0),
562 13 : kbtpd(-1.0),
563 13 : TAMD(0),
564 13 : alpha(0.0),
565 13 : gamma(0.0),
566 13 : kappa(getNumberOfArguments(),0.0),
567 13 : fict_max(getNumberOfArguments(),0.0),
568 13 : fict_min(getNumberOfArguments(),0.0),
569 13 : fict (getNumberOfArguments(),-999.0), // -999 means no initial values given in plumed.dat
570 13 : vfict(getNumberOfArguments(),0.0),
571 13 : mfict(getNumberOfArguments(),10000.0),
572 13 : xeta(0.0),
573 13 : veta(0.0),
574 13 : meta(0.0),
575 13 : phivs(0.0),
576 13 : work(0.0),
577 13 : weight(0.0),
578 13 : flog(0.0),
579 13 : hlog(0.0),
580 13 : ffict(getNumberOfArguments(),0.0),
581 13 : fict_ave(getNumberOfArguments(),0.0),
582 13 : fictValue(getNumberOfArguments(),NULL),
583 26 : vfictValue(getNumberOfArguments(),NULL) {
584 13 : backup.fict.resize(getNumberOfArguments(),0.0);
585 13 : backup.vfict.resize(getNumberOfArguments(),0.0);
586 13 : backup.ffict.resize(getNumberOfArguments(),0.0);
587 13 : backup.xeta = 0.0;
588 13 : backup.veta = 0.0;
589 13 : backup.phivs = 0.0;
590 13 : backup.work = 0.0;
591 13 : backup.weight = 0.0;
592 :
593 13 : parse("INTERVAL",interval);
594 13 : parse("DELTA_T",delta_t);
595 13 : parse("THERMOSTAT",thermostat);
596 13 : kbt = getkBT(); // read as temperature
597 13 : parse("TEMPPD",kbtpd); // read as temperature
598 :
599 13 : parse("TAMD",TAMD);
600 13 : parse("ALPHA",alpha);
601 13 : parse("GAMMA",gamma);
602 13 : parseVector("KAPPA",kappa);
603 :
604 13 : parseVector("FICT_MAX",fict_max);
605 13 : parseVector("FICT_MIN",fict_min);
606 :
607 13 : parseVector("FICT",fict);
608 13 : parseVector("VFICT",vfict);
609 26 : parseVector("MFICT",mfict);
610 : // FICT, VFICT and MFICT are optional, and parseVector empties the vector of an
611 : // optional keyword that is absent: put back the defaults.
612 13 : if( fict.empty() ) {
613 10 : fict.assign(getNumberOfArguments(),-999.0);
614 : }
615 13 : if( vfict.empty() ) {
616 8 : vfict.assign(getNumberOfArguments(),0.0);
617 : }
618 13 : if( mfict.empty() ) {
619 2 : mfict.assign(getNumberOfArguments(),10000.0);
620 : }
621 :
622 13 : parse("XETA",xeta);
623 13 : parse("VETA",veta);
624 13 : parse("META",meta);
625 :
626 13 : parse("FLOG",flog);
627 :
628 : // read initial value of work for each replica of LogPD.
629 : // multi_sim_comm is only valid on the first process of each replica (see Action::multi_sim_comm),
630 : // so the number of replicas and the replica index are taken from getNumberOfReplicas()
631 : // and getReplicaIndex(), which are correct on every process.
632 13 : if( getNumberOfReplicas()>1 ) {
633 10 : vector<double> vwork(getNumberOfReplicas(),0.0);
634 20 : parseVector("WORK",vwork);
635 : // WORK is optional and defaults to 0: parseVector empties vwork if it is absent
636 10 : if( vwork.empty() ) {
637 2 : vwork.assign(getNumberOfReplicas(),0.0);
638 : }
639 : // initial work of this replica
640 10 : work = vwork[getReplicaIndex()];
641 : } else {
642 3 : work = 0.0;
643 : }
644 :
645 13 : if( kbtpd>=0.0 ) {
646 0 : kbtpd *= getKBoltzmann();
647 : } else {
648 13 : kbtpd = kbt;
649 : }
650 :
651 13 : if( meta == 0.0 ) {
652 12 : const double nkt = getNumberOfArguments()*kbt;
653 12 : meta = nkt*100.0*100.0;
654 : }
655 :
656 13 : if(alpha == 0.0 && gamma == 0.0) {
657 10 : alpha = 4.0;
658 10 : gamma = 1 / alpha;
659 3 : } else if(alpha != 0.0 && gamma == 0.0) {
660 3 : gamma = 1 / alpha;
661 0 : } else if(alpha == 0.0 && gamma != 0.0) {
662 0 : alpha = 1 / gamma;
663 : }
664 :
665 13 : checkRead();
666 :
667 : // output messaages to Plumed's log file
668 13 : if( getNumberOfReplicas()>1 ) {
669 10 : if( TAMD ) {
670 0 : log.printf("TAMD-PD, replica parallel of TAMD, no logarithmic flattening.\n");
671 : } else {
672 10 : log.printf("LogPD, replica parallel of LogMFD.\n");
673 : }
674 10 : log.printf("number of replica : %d.\n", getNumberOfReplicas() );
675 : } else {
676 3 : if( TAMD ) {
677 0 : log.printf("TAMD, no logarithmic flattening.\n");
678 : } else {
679 3 : log.printf("LogMFD, logarithmic flattening.\n");
680 : }
681 : }
682 :
683 13 : log.printf(" with harmonic force constant ");
684 36 : for(unsigned i=0; i<kappa.size(); i++) {
685 23 : log.printf(" %f",kappa[i]);
686 : }
687 13 : log.printf("\n");
688 :
689 13 : log.printf(" with interval of cv(ideal) update ");
690 13 : log.printf(" %d", interval);
691 13 : log.printf("\n");
692 :
693 13 : log.printf(" with time step of cv(ideal) update ");
694 13 : log.printf(" %f", delta_t);
695 13 : log.printf("\n");
696 :
697 13 : if( !TAMD ) {
698 13 : log.printf(" with alpha, gamma ");
699 13 : log.printf(" %f %f", alpha, gamma);
700 13 : log.printf("\n");
701 : }
702 :
703 13 : log.printf(" with Thermostat for cv(ideal) ");
704 13 : log.printf(" %s", thermostat.c_str());
705 13 : log.printf("\n");
706 :
707 13 : log.printf(" with initial free energy ");
708 13 : log.printf(" %f", flog);
709 13 : log.printf("\n");
710 :
711 13 : log.printf(" with mass of cv(ideal)");
712 36 : for(unsigned i=0; i<mfict.size(); i++) {
713 23 : log.printf(" %f", mfict[i]);
714 : }
715 13 : log.printf("\n");
716 :
717 13 : log.printf(" with initial value of cv(ideal)");
718 36 : for(unsigned i=0; i<fict.size(); i++) {
719 23 : log.printf(" %f", fict[i]);
720 : }
721 13 : log.printf("\n");
722 :
723 13 : log.printf(" with initial velocity of cv(ideal)");
724 36 : for(unsigned i=0; i<vfict.size(); i++) {
725 23 : log.printf(" %f", vfict[i]);
726 : }
727 13 : log.printf("\n");
728 :
729 13 : log.printf(" with maximum value of cv(ideal) ");
730 36 : for(unsigned i=0; i<fict_max.size(); i++) {
731 23 : log.printf(" %f",fict_max[i]);
732 : }
733 13 : log.printf("\n");
734 :
735 13 : log.printf(" with minimum value of cv(ideal) ");
736 36 : for(unsigned i=0; i<fict_min.size(); i++) {
737 23 : log.printf(" %f",fict_min[i]);
738 : }
739 13 : log.printf("\n");
740 :
741 13 : log.printf(" and kbt ");
742 13 : log.printf(" %f\n",kbt);
743 13 : log.printf(" kbt for PD %f\n",kbtpd);
744 :
745 : // setup Value* variables
746 36 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
747 23 : std::string comp = getPntrToArgument(i)->getName()+"_fict";
748 23 : addComponentWithDerivatives(comp);
749 :
750 23 : if(getPntrToArgument(i)->isPeriodic()) {
751 : std::string a,b;
752 0 : getPntrToArgument(i)->getDomain(a,b);
753 0 : componentIsPeriodic(comp,a,b);
754 : } else {
755 23 : componentIsNotPeriodic(comp);
756 : }
757 23 : fictValue[i] = getPntrToComponent(comp);
758 :
759 46 : comp = getPntrToArgument(i)->getName()+"_vfict";
760 23 : addComponent(comp);
761 :
762 23 : componentIsNotPeriodic(comp);
763 23 : vfictValue[i] = getPntrToComponent(comp);
764 : }
765 13 : }
766 :
767 : /**
768 : \brief calculate forces for fictitious variables at every MD steps.
769 : \details This function calculates initial values of fictitious variables
770 : and write header messages to LogMFD log files at the first MFD step,
771 : calculates restraining fources comes from difference between the fictitious variable
772 : and collective variable at every MD steps.
773 : */
774 1740 : void LogMFD::calculate() {
775 1740 : if( firsttime ) {
776 13 : firsttime = false;
777 :
778 13 : step_initial = getStep();
779 :
780 : // set initial values of fictitious variables if they were not specified.
781 36 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
782 23 : if( fict[i] != -999.0 ) {
783 3 : continue; // -999 means no initial values given in plumed.dat
784 : }
785 :
786 : // use the collective variables as the initial of the fictitious variable.
787 20 : fict[i] = getArgument(i);
788 :
789 : // average values of fictitious variables by all replica.
790 : // multi_sim_comm is only valid on the first process of each replica:
791 : // sum there and broadcast the result to the other processes of the replica.
792 20 : if( getNumberOfReplicas()>1 ) {
793 20 : if( comm.Get_rank()==0 ) {
794 16 : multi_sim_comm.Sum(fict[i]);
795 : }
796 20 : comm.Bcast(fict[i],0);
797 20 : fict[i] /= getNumberOfReplicas();
798 : }
799 : }
800 :
801 : // initialize accumulation value to zero
802 36 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
803 23 : fict_ave[i] = 0.0;
804 : }
805 :
806 : // calculate invariant for NVE
807 13 : if(thermostat == "NVE") {
808 : // kinetic energy
809 11 : const double ekin = calcEkin();
810 : // potential energy
811 21 : const double pot = TAMD ? flog : sgn(flog)*gamma * std::log1p( alpha*fabs(flog) );
812 : // invariant
813 11 : hlog = pot + ekin;
814 2 : } else if(thermostat == "NVT") {
815 1 : const double nkt = getNumberOfArguments()*kbt;
816 : // kinetic energy
817 1 : const double ekin = calcEkin();
818 : // bath energy
819 1 : const double ekin_bath = 0.5*veta*veta*meta + xeta*nkt;
820 : // potential energy
821 1 : const double pot = TAMD ? flog : sgn(flog)*gamma * std::log1p( alpha*fabs(flog) );
822 : // invariant
823 1 : hlog = pot + ekin + ekin_bath;
824 1 : } else if(thermostat == "VS") {
825 : // kinetic energy
826 1 : const double ekin = calcEkin();
827 1 : if( ekin == 0.0 ) { // this means VFICT is not given.
828 : // initial velocities
829 2 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
830 1 : vfict[i] = sqrt(kbt/mfict[i]);
831 : }
832 : } else {
833 0 : const double nkt = getNumberOfArguments()*kbt;
834 0 : const double svs = sqrt(nkt/ekin/2);
835 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
836 0 : vfict[i] *= svs; // scale velocities
837 : }
838 : }
839 : // initial VS potential
840 1 : phivs = TAMD ? flog : sgn(flog)* gamma*std::log1p( alpha*fabs(flog) );
841 :
842 : // invariant
843 1 : hlog = 0.0;
844 : }
845 :
846 13 : weight = 1.0; // for replica parallel
847 :
848 : // open LogMFD's log file
849 13 : if( comm.Get_rank()==0 && getReplicaIndex()==0 ) {
850 7 : FILE *outlog = std::fopen("logmfd.out", "w");
851 :
852 : // output messages to LogMFD's log file
853 7 : if( getNumberOfReplicas()>1 ) {
854 : fprintf(outlog, "# LogPD, replica parallel of LogMFD\n");
855 4 : fprintf(outlog, "# number of replica : %d\n", getNumberOfReplicas() );
856 : } else {
857 : fprintf(outlog, "# LogMFD\n");
858 : }
859 :
860 : fprintf(outlog, "# CVs :");
861 18 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
862 : fprintf(outlog, " %s", getPntrToArgument(i)->getName().c_str() );
863 : }
864 : fprintf(outlog, "\n");
865 :
866 : fprintf(outlog, "# Mass for CV particles :");
867 18 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
868 11 : fprintf(outlog, "%15.6f", mfict[i]);
869 : }
870 : fprintf(outlog, "\n");
871 :
872 : fprintf(outlog, "# Mass for thermostat :");
873 7 : fprintf(outlog, "%15.6f", meta);
874 : fprintf(outlog, "\n");
875 : fprintf(outlog, "# 1:iter_mfd, 2:Flog, 3:2*Ekin/gkb[K], 4:eta, 5:Veta,\n");
876 :
877 18 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
878 11 : fprintf(outlog, "# %u:%s_fict(t), %u:%s_vfict(t), %u:%s_force(t),\n",
879 : 6+i*3, getPntrToArgument(i)->getName().c_str(),
880 : 7+i*3, getPntrToArgument(i)->getName().c_str(),
881 11 : 8+i*3, getPntrToArgument(i)->getName().c_str() );
882 : }
883 7 : fclose(outlog);
884 : }
885 :
886 13 : if( comm.Get_rank()==0 ) {
887 : // the number of replica is added to file name to distingwish replica.
888 11 : FILE *outlog2 = fopen("replica.out", "w");
889 11 : fprintf(outlog2, "# Replica No. %d of %d.\n",
890 : getReplicaIndex(), getNumberOfReplicas() );
891 :
892 : fprintf(outlog2, "# 1:iter_mfd, 2:work, 3:weight,\n");
893 30 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
894 19 : fprintf(outlog2, "# %u:%s(q)\n",
895 : 4+i, getPntrToArgument(i)->getName().c_str() );
896 : }
897 11 : fclose(outlog2);
898 : }
899 :
900 : // output messages to Plumed's log file
901 : // log.printf("LOGMFD thermostat parameters Xeta Veta Meta");
902 : // log.printf(" %f %f %f", xeta, veta, meta);
903 : // log.printf("\n");
904 : // log.printf("# 1:iter_mfd, 2:Flog, 3:2*Ekin/gkb[K], 4:eta, 5:Veta,");
905 : // log.printf("# 6:X1(t), 7:V1(t), 8:F1(t), 9:X2(t), 10:V2(t), 11:F2(t), ...");
906 :
907 : } // firsttime
908 :
909 : // calculate force for fictitious variable
910 : double ene=0.0;
911 3720 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
912 : // difference between fictitious variable and collective variable.
913 1980 : const double diff = difference(i,fict[i],getArgument(i));
914 : // restraining force.
915 1980 : const double f = -kappa[i]*diff;
916 1980 : setOutputForce(i,f);
917 :
918 : // restraining energy.
919 1980 : ene += 0.5*kappa[i]*diff*diff;
920 :
921 : // accumulate force, later it will be averaged.
922 1980 : ffict[i] += -f;
923 :
924 : // accumulate varience of collective variable, later it will be averaged.
925 1980 : fict_ave[i] += diff;
926 : }
927 :
928 1740 : setBias(ene);
929 3720 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
930 : // correct fict so that it is inside [min:max].
931 1980 : fict[i] = fictValue[i]->bringBackInPbc(fict[i]);
932 1980 : fictValue[i]->set(fict[i]);
933 1980 : vfictValue[i]->set(vfict[i]);
934 : }
935 1740 : } // calculate
936 :
937 : /**
938 : \brief update fictitious variables.
939 : \details This function manages evolution of fictitious variables.
940 : This function calculates mean force, updates fictitious variables by one MFD step,
941 : bounces back variables, updates free energy, and record logs.
942 : */
943 1740 : void LogMFD::update() {
944 1740 : if( (getStep()-step_initial)%interval != interval-1 ) {
945 : return;
946 : }
947 :
948 210 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
949 135 : backup.fict[i] = fict[i];
950 135 : backup.vfict[i] = vfict[i];
951 : }
952 75 : backup.xeta = xeta;
953 75 : backup.veta = veta;
954 75 : backup.phivs = phivs;
955 75 : backup.work = work;
956 75 : backup.weight = weight;
957 :
958 : // calc mean force for fictitious variables
959 75 : calcMeanForce();
960 :
961 : // record log for fictitious variables
962 75 : if( comm.Get_rank()==0 && getReplicaIndex()==0 ) {
963 39 : const double ekin = calcEkin();
964 39 : const double temp = 2.0*ekin/getNumberOfArguments()/getKBoltzmann();
965 :
966 39 : FILE *outlog = std::fopen("logmfd.out", "a");
967 39 : fprintf(outlog, "%*d", 8, (int)(getStep()-step_initial)/interval);
968 39 : fprintf(outlog, "%15.6f", flog);
969 : fprintf(outlog, "%15.6f", temp);
970 39 : fprintf(outlog, "%15.6f", xeta);
971 39 : fprintf(outlog, "%15.6f", veta);
972 102 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
973 63 : fprintf(outlog, "%15.6f", fict[i]);
974 63 : fprintf(outlog, "%15.6f", vfict[i]);
975 63 : fprintf(outlog, "%15.6f", ffict[i]);
976 : }
977 : fprintf(outlog," \n");
978 39 : fclose(outlog);
979 : }
980 :
981 : // record log for collective variables
982 75 : if( comm.Get_rank()==0 ) {
983 : // the number of replica is added to file name to distingwish replica.
984 63 : FILE *outlog2 = fopen("replica.out", "a");
985 63 : fprintf(outlog2, "%*d", 8, (int)(getStep()-step_initial)/interval);
986 63 : fprintf(outlog2, "%16.6e ", work);
987 63 : fprintf(outlog2, "%16.6e ", weight);
988 174 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
989 111 : fprintf(outlog2, "%15.6f", fict_ave[i]);
990 : }
991 : fprintf(outlog2," \n");
992 63 : fclose(outlog2);
993 : }
994 :
995 : // update fictitious variables
996 75 : if(thermostat == "NVE") {
997 65 : updateNVE();
998 10 : } else if(thermostat == "NVT") {
999 5 : updateNVT();
1000 5 : } else if(thermostat == "VS") {
1001 5 : updateVS();
1002 : }
1003 :
1004 : // update work done by fictitious dynamical variables
1005 75 : updateWork();
1006 :
1007 : // check boundary
1008 : bool reject = false;
1009 210 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1010 135 : if( fict[i] < fict_min[i] || fict_max[i] < fict[i] ) {
1011 : reject = true;
1012 0 : backup.vfict[i] *= -1.0;
1013 : }
1014 : }
1015 75 : if( reject ) {
1016 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1017 0 : fict[i] = backup.fict[i];
1018 0 : vfict[i] = backup.vfict[i];
1019 : }
1020 0 : xeta = backup.xeta;
1021 0 : veta = backup.veta;
1022 0 : phivs = backup.phivs;
1023 0 : work = backup.work;
1024 0 : weight = backup.weight;
1025 : }
1026 :
1027 : // update free energy
1028 75 : flog = calcFlog();
1029 :
1030 : // reset mean force
1031 210 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1032 135 : ffict[i] = 0.0;
1033 135 : fict_ave[i] = 0.0;
1034 : }
1035 :
1036 : } // update
1037 :
1038 : /**
1039 : \brief update fictitious variables by NVE mechanics.
1040 : \details This function updates ficitious variables by one NVE-MFD step using mean forces
1041 : and flattening coefficient and free energy.
1042 : */
1043 65 : void LogMFD::updateNVE() {
1044 65 : const double dt = delta_t;
1045 :
1046 : // get latest free energy and flattening coefficient
1047 65 : flog = calcFlog();
1048 65 : const double clog = calcClog();
1049 :
1050 : // update all ficitious variables by one MFD step
1051 190 : for( unsigned i=0; i<getNumberOfArguments(); ++i ) {
1052 : // update velocity (full step)
1053 125 : vfict[i]+=clog*ffict[i]*dt/mfict[i];
1054 : // update position (full step)
1055 125 : fict[i]+=vfict[i]*dt;
1056 : }
1057 65 : } // updateNVE
1058 :
1059 : /**
1060 : \brief update fictitious variables by NVT mechanics.
1061 : \details This function updates ficitious variables by one NVT-MFD step using mean forces
1062 : and flattening coefficient and free energy.
1063 : */
1064 5 : void LogMFD::updateNVT() {
1065 5 : const double dt = delta_t;
1066 5 : const double nkt = getNumberOfArguments()*kbt;
1067 :
1068 : // backup vfict
1069 10 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1070 5 : backup.vfict[i] = vfict[i];
1071 : }
1072 :
1073 : const int niter=5;
1074 30 : for(unsigned j=0; j<niter; ++j) {
1075 : // get latest free energy and flattening coefficient
1076 25 : flog = calcFlog();
1077 25 : const double clog = calcClog();
1078 :
1079 : // restore vfict from backup.vfict
1080 50 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1081 25 : vfict[i] = backup.vfict[i];
1082 : }
1083 :
1084 : // evolve vfict from backup.vfict by dt/2
1085 50 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1086 25 : vfict[i] *= exp(-0.25*dt*veta);
1087 25 : vfict[i] += 0.5*dt*clog*ffict[i]/mfict[i];
1088 25 : vfict[i] *= exp(-0.25*dt*veta);
1089 : }
1090 : }
1091 :
1092 : // get latest free energy and flattening coefficient
1093 5 : flog = calcFlog();
1094 5 : const double clog = calcClog();
1095 :
1096 : // evolve vfict by dt/2, and evolve fict by dt
1097 10 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1098 5 : vfict[i] *= exp(-0.25*dt*veta);
1099 5 : vfict[i] += 0.5*dt*clog*ffict[i]/mfict[i];
1100 5 : vfict[i] *= exp(-0.25*dt*veta);
1101 5 : fict[i] += dt*vfict[i];
1102 : }
1103 :
1104 : // evolve xeta and veta by dt
1105 5 : xeta += 0.5*dt*veta;
1106 5 : const double ekin = calcEkin();
1107 5 : veta += dt*(2.0*ekin-nkt)/meta;
1108 5 : xeta += 0.5*dt*veta;
1109 5 : } // updateNVT
1110 :
1111 : /**
1112 : \brief update fictitious variables by VS mechanics.
1113 : \details This function updates ficitious variables by one VS-MFD step using mean forces
1114 : and flattening coefficient and free energy.
1115 : */
1116 5 : void LogMFD::updateVS() {
1117 5 : const double dt = delta_t;
1118 5 : const double nkt = getNumberOfArguments()*kbt;
1119 :
1120 : // get latest free energy and flattening coefficient
1121 5 : flog = calcFlog();
1122 5 : const double clog = calcClog();
1123 :
1124 10 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1125 : // update velocity (full step)
1126 5 : vfict[i] += clog*ffict[i]*dt/mfict[i];
1127 : }
1128 :
1129 5 : const double ekin = calcEkin();
1130 5 : const double svs = sqrt(nkt/ekin/2);
1131 :
1132 10 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1133 : // update position (full step)
1134 5 : vfict[i] *= svs;
1135 5 : fict[i] += vfict[i]*dt;
1136 : }
1137 :
1138 : // evolve VS potential
1139 5 : phivs += nkt*std::log(svs);
1140 5 : } // updateVS
1141 :
1142 : /**
1143 : \brief update work done by fictious variables.
1144 : \details This function updates work done by ficitious variables.
1145 : */
1146 75 : void LogMFD::updateWork() {
1147 : // accumulate work, it was initialized as 0.0
1148 210 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1149 135 : work -= backup.ffict[i] * vfict[i] * delta_t;
1150 : }
1151 75 : } // updateWork
1152 :
1153 : /**
1154 : \brief calculate mean force for fictitious variables.
1155 : \details This function calculates mean forces by averaging forces accumulated during one MFD step,
1156 : update work variables done by fictitious variables by one MFD step,
1157 : calculate weight variable of this replica for LogPD.
1158 : */
1159 75 : void LogMFD::calcMeanForce() {
1160 : // cale temporal mean force for each CV
1161 210 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1162 135 : ffict[i] /= interval;
1163 : // backup force of replica
1164 135 : backup.ffict[i] = ffict[i];
1165 : // average of diff (getArgument(i)-fict[i])
1166 135 : fict_ave[i] /= interval;
1167 : // average of getArgument(i)
1168 135 : fict_ave[i] += fict[i];
1169 :
1170 : // correct fict_ave so that it is inside [min:max].
1171 135 : fict_ave[i] = fictValue[i]->bringBackInPbc(fict_ave[i]);
1172 : }
1173 :
1174 : // for replica parallel.
1175 : // multi_sim_comm is only valid on the first process of each replica:
1176 : // reduce there and broadcast the results to the other processes of the replica,
1177 : // so that all of them evolve the fictitious variables in the same way.
1178 75 : if( getNumberOfReplicas()>1 ) {
1179 : // find the minimum work among all replicas
1180 60 : double work_min = work;
1181 60 : if( comm.Get_rank()==0 ) {
1182 48 : multi_sim_comm.Min(work_min);
1183 : }
1184 60 : comm.Bcast(work_min,0);
1185 :
1186 : // weight of this replica.
1187 : // here, work is reduced by work_min to avoid all exp(-work/kbt)s disconverge
1188 60 : if( kbtpd == 0.0 ) {
1189 0 : weight = work==work_min ? 1.0 : 0.0;
1190 : } else {
1191 60 : weight = exp(-(work-work_min)/kbtpd);
1192 : }
1193 :
1194 : // normalize the weight
1195 60 : double sum_weight = weight;
1196 60 : if( comm.Get_rank()==0 ) {
1197 48 : multi_sim_comm.Sum(sum_weight);
1198 : }
1199 60 : comm.Bcast(sum_weight,0);
1200 60 : weight /= sum_weight;
1201 :
1202 : // weighting force of this replica
1203 180 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1204 120 : ffict[i] *= weight;
1205 : }
1206 :
1207 : // averaged mean forces of all replica.
1208 60 : if( comm.Get_rank()==0 ) {
1209 48 : multi_sim_comm.Sum(ffict);
1210 : }
1211 60 : comm.Bcast(ffict,0);
1212 : // now, mean force is obtained.
1213 : }
1214 75 : } // calcMeanForce
1215 :
1216 : /**
1217 : \brief calculate kinetic energy of fictitious variables.
1218 : \retval kinetic energy.
1219 : \details This function calculates sum of kinetic energy of all fictitious variables.
1220 : */
1221 237 : double LogMFD::calcEkin() {
1222 : double ekin=0.0;
1223 628 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1224 391 : ekin += mfict[i]*vfict[i]*vfict[i]*0.5;
1225 : }
1226 237 : return ekin;
1227 : } // calcEkin
1228 :
1229 : /**
1230 : \brief calculate free energy of fictitious variables.
1231 : \retval free energy.
1232 : \details This function calculates free energy by using invariant of canonical mechanics.
1233 : */
1234 175 : double LogMFD::calcFlog() {
1235 175 : const double nkt = getNumberOfArguments()*kbt;
1236 175 : const double ekin = calcEkin();
1237 : double pot;
1238 :
1239 175 : if (thermostat == "NVE") {
1240 130 : pot = hlog - ekin;
1241 45 : } else if (thermostat == "NVT") {
1242 35 : const double ekin_bath = 0.5*veta*veta*meta+xeta*nkt;
1243 35 : pot = hlog - ekin - ekin_bath;
1244 10 : } else if (thermostat == "VS") {
1245 10 : pot = phivs;
1246 : } else {
1247 : pot = 0.0; // never occurs
1248 : }
1249 :
1250 350 : return TAMD ? pot : sgn(pot)*expm1(fabs(pot)/gamma)/alpha;
1251 : } // calcFlog
1252 :
1253 : /**
1254 : \brief calculate coefficient for flattening.
1255 : \retval flattering coefficient.
1256 : \details This function returns 1.0 for TAMD, flattening coefficient for LogMFD.
1257 : */
1258 100 : double LogMFD::calcClog() {
1259 100 : return TAMD ? 1.0 : alpha*gamma/(alpha*fabs(flog)+1.0);
1260 : } // calcClog
1261 :
1262 : } // logmfd
1263 : } // PLMD
|
