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