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=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 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, reconstructing free energy profiles can be performed more efficiently (requiring less elapsed computing time) in LogPD than with a single MD system in LogMFD. In the case that there exists 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, and NAMD.
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) &=& V_{X_i}^\prime \left( {{t_n}} \right) + \Delta t\left[ { - \left( {\frac{{\alpha \gamma }}{{\alpha F\left( {{t_n}} \right) + 1}}} \right)\frac{{\partial F\left( {{t_n}} \right)}}{{\partial {X_i}}}} \right]\\
193 : S\left( {{t_{n + 1}}} \right) &=& \sqrt {\frac{{N{k_B}T}}{{\sum\limits_i {{M_i}V_{{X_i}}^2\left( {{t_{n + 1}}} \right)} }}} \\
194 : {V_{{X_i}}}^\prime \left( {{t_{n + 1}}} \right) &=& S\left( {{t_{n + 1}}} \right){V_{{X_i}}}\left( {{t_{n + 1}}} \right)\\
195 : {X_i}\left( {{t_{n + 1}}} \right) &=& {X_i}\left( {{t_n}} \right) + \Delta t V_{X_i}^\prime \left( {{t_{n + 1}}} \right)\\
196 : {\Psi_{\rm log}}\left( {{t_{n + 1}}} \right) &=& N{k_B}T\log S\left( {{t_{n + 1}}} \right) + {\Psi_{\rm log}}\left( {{t_n}} \right)\\
197 : F\left( {{t_{n + 1}}} \right) &=& \frac{1}{\alpha} \left[
198 : \exp \left\{ \Psi_{\rm log} \left( t_{n+1} \right) / \gamma \right\} - 1 \right]
199 : \f}
200 :
201 : Note that \f$V_{X_i}^\prime\left( {{t_0}} \right)\f$ is assumed to be initially given, which meets the following relation,
202 :
203 : \f[
204 : \sum\limits_{i = 1}^N M_i V_{X_i}^{\prime 2} \left( t_0 \right) = N{k_B}{T}
205 : \f]
206 :
207 : 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.
208 :
209 :
210 : \par Examples
211 : \section Examples Examples
212 :
213 : \subsection Example-LoGMFD Example of LogMFD
214 :
215 : The following input file tells plumed to restrain collective variables
216 : to the fictitious dynamical variables in LogMFD/PD.
217 :
218 : plumed.dat
219 : \plumedfile
220 : UNITS TIME=fs LENGTH=1.0 ENERGY=kcal/mol MASS=1.0 CHARGE=1.0
221 : phi: TORSION ATOMS=5,7,9,15
222 : psi: TORSION ATOMS=7,9,15,17
223 :
224 : # LogMFD
225 : LOGMFD ...
226 : LABEL=logmfd
227 : ARG=phi,psi
228 : KAPPA=100.0,100.0
229 : DELTA_T=0.5
230 : INTERVAL=500
231 : TEMP=300.0
232 : FLOG=5.0
233 : MFICT=5000000.0,5000000.0
234 : # FICT=0.8,0.8
235 : VFICT=3.5e-4,3.5e-4
236 : ALPHA=4.0
237 : THERMOSTAT=NVE
238 : VETA=0.0
239 : META=20000.0
240 : FICT_MAX=3.1,3.1
241 : FICT_MIN=-3.1,-3.1
242 : ... LOGMFD
243 : \endplumedfile
244 :
245 : To submit this simulation with Gromacs, use the following command line
246 : to execute a LogMFD run with Gromacs-MD.
247 : Here TOPO/topol0.tpr is an input file
248 : which contains atomic coordinates and Gromacs parameters.
249 :
250 : \verbatim
251 : gmx_mpi mdrun -s TOPO/topol0.tpr -plumed
252 : \endverbatim
253 :
254 : This command will output files named logmfd.out and replica.out.
255 :
256 : The output file logmfd.out records free energy and all fictitious dynamical variables at every MFD step.
257 :
258 : logmfd.out
259 :
260 : \verbatim
261 : # LogMFD
262 : # CVs : phi psi
263 : # Mass for CV particles : 5000000.000000000 5000000.000000000
264 : # Mass for thermostat : 20000.000000000
265 : # 1:iter_md, 2:Flog, 3:2*Ekin/gkb[K], 4:eta, 5:Veta,
266 : # 6:phi_fict(t), 7:phi_vfict(t), 8:phi_force(t),
267 : # 9:psi_fict(t), 10:psi_vfict(t), 11:psi_force(t),
268 : 1 4.99918574 308.24149708 0.00000000 0.00000000 -2.85605938 0.00035002 5.19074544 2.79216364 0.00035000 -0.53762989
269 : 2 4.99836196 308.26124159 0.00000000 0.00000000 -2.85588436 0.00035005 4.71247605 2.79233863 0.00035000 -0.00532474
270 : 3 4.99743572 308.28344595 0.00000000 0.00000000 -2.85570932 0.00035007 5.34358230 2.79251363 0.00035000 -0.05119816
271 : ...
272 : \endverbatim
273 :
274 : The output file replica.out records all collective variables at every MFD step.
275 :
276 : replica.out
277 :
278 : \verbatim
279 : # Replica No. 0 of 1.
280 : # 1:iter_md, 2:work, 3:weight,
281 : # 4:phi(q)
282 : # 5:psi(q)
283 : 1 -8.142952e-04 1.000000e+00 -2.80432694 2.78661234
284 : 2 -1.638105e-03 1.000000e+00 -2.80893462 2.79211039
285 : 3 -2.564398e-03 1.000000e+00 -2.80244854 2.79182665
286 : ...
287 : \endverbatim
288 :
289 : \subsection Example-LogPD Example of LogPD
290 :
291 : Use the following command line to execute a LogPD run using two MD replicas (note that only Gromacs is currently available for LogPD).
292 : Here TOPO/topol0.tpr and TOPO/topol1.tpr are input files
293 : which contain atomic coordinates of each replica and Gromacs parameters.
294 :
295 : \verbatim
296 : mpirun -np 2 gmx_mpi mdrun -s TOPO/topol -plumed -multi 2
297 : \endverbatim
298 :
299 : This command will output files named logmfd.out, replica.out.0 and replica.out.1.
300 :
301 : The output file logmfd.out records free energy and all fictitious dynamical variables at every MFD step.
302 :
303 : logmfd.out
304 :
305 : \verbatim
306 : # LogPD, replica parallel of LogMFD
307 : # number of replica : 2
308 : # CVs : phi psi
309 : # Mass for CV particles : 5000000.000000000 5000000.000000000
310 : # Mass for thermostat : 20000.000000000
311 : # 1:iter_md, 2:Flog, 3:2*Ekin/gkb[K], 4:eta, 5:Veta,
312 : # 6:phi_fict(t), 7:phi_vfict(t), 8:phi_force(t),
313 : # 9:psi_fict(t), 10:psi_vfict(t), 11:psi_force(t),
314 : 1 5.00224715 308.16814691 0.00000000 0.00000000 -0.95937173 0.00034994 -12.91277494 0.78923967 0.00035000 0.07353010
315 : 2 5.00476934 308.10774854 0.00000000 0.00000000 -0.95919679 0.00034989 -11.20093553 0.78941467 0.00034999 -3.21098229
316 : 3 5.00702463 308.05376594 0.00000000 0.00000000 -0.95902187 0.00034983 -10.81712171 0.78958965 0.00034998 -2.07196718
317 : ...
318 : \endverbatim
319 :
320 :
321 : The output file replica.out.0 records all collective variables of replica No.0 at every MFD step.
322 :
323 : replica.out.0
324 :
325 : \verbatim
326 : # Replica No. 0 of 2.
327 : # 1:iter_md, 2:work, 3:weight,
328 : # 4:phi(q)
329 : # 5:psi(q)
330 : 1 1.843110e-03 5.003389e-01 -1.10929125 0.83348865
331 : 2 3.466179e-03 5.010942e-01 -1.05020764 0.78731283
332 : 3 4.927870e-03 5.017619e-01 -1.04968867 0.79635198
333 : ...
334 : \endverbatim
335 :
336 : The output file replica.out.1 records all collective variables of replica No.1 at every MFD step.
337 :
338 : replica.out.1
339 :
340 : \verbatim
341 : # Replica No. 1 of 2.
342 : # 1:iter_md, 2:work, 3:weight,
343 : # 4:phi(q)
344 : # 5:psi(q)
345 : 1 2.651173e-03 4.996611e-01 -1.06802968 0.74605205
346 : 2 6.075530e-03 4.989058e-01 -1.09264741 0.72681448
347 : 3 9.129358e-03 4.982381e-01 -1.08517238 0.74084241
348 : ...
349 : \endverbatim
350 :
351 : */
352 : //+ENDPLUMEDOC
353 :
354 : #include "bias/Bias.h"
355 : #include "core/ActionRegister.h"
356 : #include "core/Atoms.h"
357 : #include "core/PlumedMain.h"
358 :
359 : #include <iostream>
360 :
361 : using namespace std;
362 : using namespace PLMD;
363 : using namespace bias;
364 :
365 : namespace PLMD {
366 : namespace logmfd {
367 : /**
368 : \brief class for LogMFD parameters, variables and subroutines.
369 : */
370 9 : class LogMFD : public Bias {
371 : bool firsttime; ///< flag that indicates first MFD step or not.
372 : int interval; ///< input parameter, period of MD steps when fictitious dynamical variables are evolved.
373 : double delta_t; ///< input parameter, one time step of MFD when fictitious dynamical variables are evolved.
374 : string thermostat; ///< input parameter, type of thermostat for canonical dyanamics.
375 : double kbt; ///< k_B*temerature
376 :
377 : int TAMD; ///< input parameter, perform TAMD instead of LogMFD.
378 : double alpha; ///< input parameter, alpha parameter for LogMFD.
379 : double gamma; ///< input parameter, gamma parameter for LogMFD.
380 : std::vector<double> kappa; ///< input parameter, strength of the harmonic restraining potential.
381 :
382 : std::vector<double> fict_max; ///< input parameter, maximum of each fictitous dynamical variable.
383 : std::vector<double> fict_min; ///< input parameter, minimum of each fictitous dynamical variable.
384 :
385 : std::vector<double> fict; ///< current values of each fictitous dynamical variable.
386 : std::vector<double> vfict; ///< current velocity of each fictitous dynamical variable.
387 : std::vector<double> mfict; ///< mass of each fictitous dynamical variable.
388 :
389 : double xeta; ///< current eta variable of thermostat.
390 : double veta; ///< current velocity of eta variable of thermostat.
391 : double meta; ///< mass of eta variable of thermostat.
392 :
393 : double flog; ///< current free energy
394 :
395 : double hlog; ///< value invariant
396 : double phivs; ///< potential used in VS method
397 : double work; ///< current works done by fictitious dynamical variables in this replica.
398 : double weight; ///< current weight of this replica.
399 :
400 : std::vector<double> ffict; ///< current force of each fictitous dynamical variable.
401 : std::vector<double> fict_ave; ///< averaged values of each collective variable.
402 :
403 : std::vector<Value*> fictValue; ///< pointers to fictitious dynamical variables
404 : std::vector<Value*> vfictValue; ///< pointers to velocity of fictitious dynamical variables
405 :
406 : public:
407 : static void registerKeywords(Keywords& keys);
408 :
409 : explicit LogMFD(const ActionOptions&);
410 : void calculate();
411 : void update();
412 : void updateNVE();
413 : void updateNVT();
414 : void updateVS();
415 : void calcMeanForce();
416 : double calcEkin();
417 : double calcFlog();
418 : double calcClog();
419 :
420 : private:
421 : double sgn( double x ) {
422 58 : return x>0.0 ? 1.0 : x<0.0 ? -1.0 : 0.0;
423 : }
424 : };
425 :
426 7362 : PLUMED_REGISTER_ACTION(LogMFD,"LOGMFD")
427 :
428 : /**
429 : \brief instruction of parameters for Plumed manual.
430 : */
431 4 : void LogMFD::registerKeywords(Keywords& keys) {
432 4 : Bias::registerKeywords(keys);
433 8 : keys.use("ARG");
434 16 : keys.add("compulsory","INTERVAL",
435 : "Period of MD steps (\\f$N_m\\f$) to update fictitious dynamical variables." );
436 16 : keys.add("compulsory","DELTA_T",
437 : "Time step for the fictitious dynamical variables (MFD step)." );
438 16 : keys.add("compulsory","THERMOSTAT",
439 : "Type of thermostat for the fictitious dynamical variables. NVE, NVT, VS are available." );
440 16 : keys.add("optional","TEMP",
441 : "Temperature of the fictitious dynamical variables in LogMFD/PD thermostat. "
442 : "If not provided or provided as 0, it will be taken from the temperature of the MD system." );
443 :
444 16 : keys.add("optional","TAMD",
445 : "When TAMD=1, TAMD/d-AFED calculations can be performed instead of LogMFD. Otherwise, the LogMFD protocol is switched on (default)." );
446 :
447 16 : keys.add("optional","ALPHA",
448 : "Alpha parameter for LogMFD. "
449 : "If not provided or provided as 0, it will be taken as 1/gamma. "
450 : "If gamma is also not provided, Alpha is set as 4, which is a sensible value when the unit of kcal/mol is used." );
451 16 : keys.add("optional","GAMMA",
452 : "Gamma parameter for LogMFD. "
453 : "If not provided or provided as 0, it will be taken as 1/alpha. "
454 : "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." );
455 16 : keys.add("compulsory","KAPPA",
456 : "Spring constant of the harmonic restraining potential for the fictitious dynamical variables." );
457 :
458 16 : keys.add("compulsory","FICT_MAX",
459 : "Maximum values reachable for the fictitious dynamical variables. The variables will elastically bounce back at the boundary (mirror boundary)." );
460 16 : keys.add("compulsory","FICT_MIN",
461 : "Minimum values reachable for the fictitious dynamical variables. The variables will elastically bounce back at the boundary (mirror boundary)." );
462 :
463 16 : keys.add("optional","FICT",
464 : "The initial values of the fictitious dynamical variables. "
465 : "If not provided, they are set equal to their corresponding CVs for the initial atomic configuration." );
466 16 : keys.add("optional","VFICT",
467 : "The initial velocities of the fictitious dynamical variables. "
468 : "If not provided, they will be taken as 0." );
469 16 : keys.add("optional","MFICT",
470 : "Masses of each fictitious dynamical variable. "
471 : "If not provided, they will be taken as 10000." );
472 :
473 16 : keys.add("optional","XETA",
474 : "The initial eta variable of the Nose-Hoover thermostat "
475 : "for the fictitious dynamical variables. "
476 : "If not provided, it will be taken as 0." );
477 16 : keys.add("optional","VETA",
478 : "The initial velocity of eta variable. "
479 : "If not provided, it will be taken as 0." );
480 16 : keys.add("optional","META",
481 : "Mass of eta variable. "
482 : "If not provided, it will be taken as \\f$N*kb*T*100*100\\f$." );
483 :
484 16 : keys.add("compulsory","FLOG",
485 : "The initial free energy value in the LogMFD/PD run."
486 : "The origin of the free energy profile is adjusted by FLOG to "
487 : "realize \\f$F({\\bf X}(t)) > 0\\f$ at any \\f${\\bf X}(t)\\f$, "
488 : "resulting in enhanced barrier-crossing. "
489 : "(The value of \\f$H_{\\rm log}\\f$ is automatically "
490 : "set according to FLOG).");
491 :
492 16 : keys.add("optional","WORK",
493 : "The initial value of the work done by fictitious dynamical "
494 : "variables in each replica. "
495 : "If not provided, it will be taken as 0.");
496 :
497 4 : componentsAreNotOptional(keys);
498 16 : keys.addOutputComponent("_fict","default",
499 : "For example, the fictitious collective variable for LogMFD is specified as "
500 : "ARG=dist12 and LABEL=logmfd in LOGMFD section in Plumed input file, "
501 : "the associated fictitious dynamical variable can be specified as "
502 : "PRINT ARG=dist12,logmfd.dist12_fict FILE=COLVAR");
503 16 : keys.addOutputComponent("_vfict","default",
504 : "For example, the fictitious collective variable for LogMFD is specified as "
505 : "ARG=dist12 and LABEL=logmfd in LOGMFD section in Plumed input file, the "
506 : "velocity of the associated fictitious dynamical variable can be specified as "
507 : "PRINT ARG=dist12,logmfd.dist12_vfict FILE=COLVAR");
508 4 : }
509 :
510 :
511 : /**
512 : \brief constructor of LogMFD class
513 : \details This constructor initializes all parameters and variables,
514 : reads input file and set value of parameters and initial value of variables,
515 : and writes messages as Plumed log.
516 : */
517 3 : LogMFD::LogMFD( const ActionOptions& ao ):
518 : PLUMED_BIAS_INIT(ao),
519 : firsttime(true),
520 : interval(10),
521 : delta_t(1.0),
522 : thermostat("NVE"),
523 : kbt(-1.0),
524 : TAMD(0),
525 : alpha(0.0),
526 : gamma(0.0),
527 : kappa(getNumberOfArguments(),0.0),
528 : fict_max(getNumberOfArguments(),0.0),
529 : fict_min(getNumberOfArguments(),0.0),
530 : fict (getNumberOfArguments(),0.0),
531 : vfict(getNumberOfArguments(),0.0),
532 : mfict(getNumberOfArguments(),10000.0),
533 : xeta(0.0),
534 : veta(0.0),
535 : meta(0.0),
536 : flog(0.0),
537 : hlog(0.0),
538 : phivs(0.0),
539 : work(0.0),
540 : weight(0.0),
541 : ffict(getNumberOfArguments(),0.0),
542 : fict_ave(getNumberOfArguments(),0.0),
543 : fictValue(getNumberOfArguments(),NULL),
544 33 : vfictValue(getNumberOfArguments(),NULL)
545 : {
546 6 : parse("INTERVAL",interval);
547 6 : parse("DELTA_T",delta_t);
548 6 : parse("THERMOSTAT",thermostat);
549 6 : parse("TEMP",kbt); // read as temperature
550 :
551 6 : parse("TAMD",TAMD);
552 6 : parse("ALPHA",alpha);
553 6 : parse("GAMMA",gamma);
554 6 : parseVector("KAPPA",kappa);
555 :
556 6 : parseVector("FICT_MAX",fict_max);
557 6 : parseVector("FICT_MIN",fict_min);
558 :
559 6 : parseVector("FICT",fict);
560 6 : parseVector("VFICT",vfict);
561 6 : parseVector("MFICT",mfict);
562 :
563 6 : parse("XETA",xeta);
564 6 : parse("VETA",veta);
565 6 : parse("META",meta);
566 :
567 6 : parse("FLOG",flog);
568 :
569 : // read initial value of work for each replica of LogPD
570 3 : if( multi_sim_comm.Get_size()>1 ) {
571 0 : vector<double> vwork(multi_sim_comm.Get_size(),0.0);
572 0 : parseVector("WORK",vwork);
573 : // initial work of this replica
574 0 : work = vwork[multi_sim_comm.Get_rank()];
575 : }
576 : else {
577 3 : work = 0.0;
578 : }
579 :
580 3 : if( kbt>=0.0 ) {
581 6 : kbt *= plumed.getAtoms().getKBoltzmann();
582 : }
583 : else {
584 0 : kbt = plumed.getAtoms().getKbT();
585 : }
586 :
587 3 : if( meta == 0.0 ) {
588 2 : const double nkt = getNumberOfArguments()*kbt;
589 2 : meta = nkt*100.0*100.0;
590 : }
591 :
592 3 : if(alpha == 0.0 && gamma == 0.0) {
593 0 : alpha = 4.0;
594 0 : gamma = 1 / alpha;
595 : }
596 3 : else if(alpha != 0.0 && gamma == 0.0) {
597 3 : gamma = 1 / alpha;
598 : }
599 0 : else if(alpha == 0.0 && gamma != 0.0) {
600 0 : alpha = 1 / gamma;
601 : }
602 :
603 3 : checkRead();
604 :
605 : // output messaages to Plumed's log file
606 3 : if( multi_sim_comm.Get_size()>1 ) {
607 0 : if( TAMD ) {
608 0 : log.printf("TAMD-PD, replica parallel of TAMD, no logarithmic flattening.\n");
609 : }
610 : else {
611 0 : log.printf("LogPD, replica parallel of LogMFD.\n");
612 : }
613 0 : log.printf("number of replica : %d.\n", multi_sim_comm.Get_size() );
614 : }
615 : else {
616 3 : if( TAMD ) {
617 0 : log.printf("TAMD, no logarithmic flattening.\n");
618 : }
619 : else {
620 3 : log.printf("LogMFD, logarithmic flattening.\n");
621 : }
622 : }
623 :
624 3 : log.printf(" with harmonic force constant ");
625 15 : for(unsigned i=0; i<kappa.size(); i++) log.printf(" %f",kappa[i]);
626 3 : log.printf("\n");
627 :
628 3 : log.printf(" with interval of cv(ideal) update ");
629 3 : log.printf(" %d", interval);
630 3 : log.printf("\n");
631 :
632 3 : log.printf(" with time step of cv(ideal) update ");
633 3 : log.printf(" %f", delta_t);
634 3 : log.printf("\n");
635 :
636 3 : if( !TAMD ) {
637 3 : log.printf(" with alpha, gamma ");
638 3 : log.printf(" %f %f", alpha, gamma);
639 3 : log.printf("\n");
640 : }
641 :
642 3 : log.printf(" with Thermostat for cv(ideal) ");
643 6 : log.printf(" %s", thermostat.c_str());
644 3 : log.printf("\n");
645 :
646 3 : log.printf(" with initial free energy ");
647 3 : log.printf(" %f", flog);
648 3 : log.printf("\n");
649 :
650 3 : log.printf(" with mass of cv(ideal)");
651 15 : for(unsigned i=0; i<mfict.size(); i++) log.printf(" %f", mfict[i]);
652 3 : log.printf("\n");
653 :
654 3 : log.printf(" with initial value of cv(ideal)");
655 15 : for(unsigned i=0; i<fict.size(); i++) log.printf(" %f", fict[i]);
656 3 : log.printf("\n");
657 :
658 3 : log.printf(" with initial velocity of cv(ideal)");
659 15 : for(unsigned i=0; i<vfict.size(); i++) log.printf(" %f", vfict[i]);
660 3 : log.printf("\n");
661 :
662 3 : log.printf(" with maximum value of cv(ideal) ");
663 15 : for(unsigned i=0; i<fict_max.size(); i++) log.printf(" %f",fict_max[i]);
664 3 : log.printf("\n");
665 :
666 3 : log.printf(" with minimum value of cv(ideal) ");
667 15 : for(unsigned i=0; i<fict_min.size(); i++) log.printf(" %f",fict_min[i]);
668 3 : log.printf("\n");
669 :
670 3 : log.printf(" and kbt ");
671 3 : log.printf(" %f",kbt);
672 3 : log.printf("\n");
673 :
674 : // setup Value* variables
675 9 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
676 3 : std::string comp = getPntrToArgument(i)->getName()+"_fict";
677 3 : addComponentWithDerivatives(comp);
678 :
679 3 : if(getPntrToArgument(i)->isPeriodic()) {
680 : std::string a,b;
681 0 : getPntrToArgument(i)->getDomain(a,b);
682 0 : componentIsPeriodic(comp,a,b);
683 : }
684 : else {
685 3 : componentIsNotPeriodic(comp);
686 : }
687 3 : fictValue[i] = getPntrToComponent(comp);
688 :
689 6 : comp = getPntrToArgument(i)->getName()+"_vfict";
690 3 : addComponent(comp);
691 :
692 3 : componentIsNotPeriodic(comp);
693 3 : vfictValue[i] = getPntrToComponent(comp);
694 : }
695 3 : }
696 :
697 : /**
698 : \brief calculate forces for fictitious variables at every MD steps.
699 : \details This function calculates initial values of fictitious variables
700 : and write header messages to LogMFD log files at the first MFD step,
701 : calculates restraining fources comes from difference between the fictious variable
702 : and collective variable at every MD steps.
703 : */
704 1500 : void LogMFD::calculate() {
705 1500 : if( firsttime ) {
706 3 : firsttime = false;
707 :
708 : // set initial values of fictitious variables if they were not specified.
709 9 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
710 6 : if( fict[i] != 0.0 ) continue;
711 :
712 : // use the collective variables as the initial of the fictitious variable.
713 0 : fict[i] = getArgument(i);
714 :
715 : // average values of fictitious variables by all replica.
716 0 : if( multi_sim_comm.Get_size()>1 ) {
717 0 : multi_sim_comm.Sum(fict[i]);
718 0 : fict[i] /= multi_sim_comm.Get_size();
719 : }
720 : }
721 :
722 : // initialize accumulation value to zero
723 9 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
724 6 : fict_ave[i] = 0.0;
725 : }
726 :
727 : // calculate invariant for NVE
728 6 : if(thermostat == "NVE") {
729 : // kinetic energy
730 1 : const double ekin = calcEkin();
731 : // potential energy
732 2 : const double pot = TAMD ? flog : sgn(flog)*gamma * std::log1p( alpha*fabs(flog) );
733 : // invariant
734 1 : hlog = pot + ekin;
735 : }
736 2 : else if(thermostat == "NVT") {
737 1 : const double nkt = getNumberOfArguments()*kbt;
738 : // kinetic energy
739 1 : const double ekin = calcEkin();
740 : // bath energy
741 1 : const double ekin_bath = 0.5*veta*veta*meta + xeta*nkt;
742 : // potential energy
743 2 : const double pot = TAMD ? flog : sgn(flog)*gamma * std::log1p( alpha*fabs(flog) );
744 : // invariant
745 1 : hlog = pot + ekin + ekin_bath;
746 : }
747 1 : else if(thermostat == "VS") {
748 : // initial velocities
749 3 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
750 3 : vfict[i] = sqrt(kbt/mfict[i]);
751 : }
752 : // initial VS potential
753 2 : phivs = TAMD ? flog : sgn(flog)* gamma*std::log1p( alpha*fabs(flog) );
754 :
755 : // invariant
756 1 : hlog = 0.0;
757 : }
758 :
759 3 : weight = 1.0; // for replica parallel
760 :
761 : // open LogMFD's log file
762 3 : if( multi_sim_comm.Get_rank()==0 && comm.Get_rank()==0 ) {
763 3 : FILE *outlog = std::fopen("logmfd.out", "w");
764 :
765 : // output messages to LogMFD's log file
766 3 : if( multi_sim_comm.Get_size()>1 ) {
767 : fprintf(outlog, "# LogPD, replica parallel of LogMFD\n");
768 0 : fprintf(outlog, "# number of replica : %d\n", multi_sim_comm.Get_size() );
769 : }
770 : else {
771 : fprintf(outlog, "# LogMFD\n");
772 : }
773 :
774 : fprintf(outlog, "# CVs :");
775 9 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
776 : fprintf(outlog, " %s", getPntrToArgument(i)->getName().c_str() );
777 : }
778 : fprintf(outlog, "\n");
779 :
780 : fprintf(outlog, "# Mass for CV particles :");
781 9 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
782 6 : fprintf(outlog, "%18.9f", mfict[i]);
783 : }
784 : fprintf(outlog, "\n");
785 :
786 : fprintf(outlog, "# Mass for thermostat :");
787 3 : fprintf(outlog, "%18.9f", meta);
788 : fprintf(outlog, "\n");
789 : fprintf(outlog, "# 1:iter_md, 2:Flog, 3:2*Ekin/gkb[K], 4:eta, 5:Veta,\n");
790 :
791 9 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
792 3 : fprintf(outlog, "# %u:%s_fict(t), %u:%s_vfict(t), %u:%s_force(t),\n",
793 : 6+i*3, getPntrToArgument(i)->getName().c_str(),
794 : 7+i*3, getPntrToArgument(i)->getName().c_str(),
795 3 : 8+i*3, getPntrToArgument(i)->getName().c_str() );
796 : }
797 :
798 3 : fclose(outlog);
799 : }
800 :
801 3 : if( comm.Get_rank()==0 ) {
802 : // the number of replica is added to file name to distingwish replica.
803 3 : FILE *outlog2 = fopen("replica.out", "w");
804 9 : fprintf(outlog2, "# Replica No. %d of %d.\n",
805 6 : multi_sim_comm.Get_rank(), multi_sim_comm.Get_size() );
806 :
807 : fprintf(outlog2, "# 1:iter_md, 2:work, 3:weight,\n");
808 9 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
809 3 : fprintf(outlog2, "# %u:%s(q)\n",
810 : 4+i, getPntrToArgument(i)->getName().c_str() );
811 : }
812 3 : fclose(outlog2);
813 : }
814 :
815 : // output messages to Plumed's log file
816 : // log.printf("LOGMFD thermostat parameters Xeta Veta Meta");
817 : // log.printf(" %f %f %f", xeta, veta, meta);
818 : // log.printf("\n");
819 : // log.printf("# 1:iter_md, 2:Flog, 3:2*Ekin/gkb[K], 4:eta, 5:Veta,");
820 : // log.printf("# 6:X1(t), 7:V1(t), 8:F1(t), 9:X2(t), 10:V2(t), 11:F2(t), ...");
821 :
822 : } // firsttime
823 :
824 : // calculate force for fictitious variable
825 : double ene=0.0;
826 4500 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
827 : // difference between fictitious variable and collective variable.
828 1500 : const double diff = difference(i,fict[i],getArgument(i));
829 : // restraining force.
830 1500 : const double f = -kappa[i]*diff;
831 : setOutputForce(i,f);
832 :
833 : // restraining energy.
834 1500 : ene += 0.5*kappa[i]*diff*diff;
835 :
836 : // accumulate force, later it will be averaged.
837 1500 : ffict[i] += -f;
838 :
839 : // accumulate varience of collective variable, later it will be averaged.
840 1500 : fict_ave[i] += diff;
841 : }
842 :
843 : setBias(ene);
844 4500 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
845 : // correct fict so that it is inside [min:max].
846 6000 : fict[i] = fictValue[i]->bringBackInPbc(fict[i]);
847 3000 : fictValue[i]->set(fict[i]);
848 3000 : vfictValue[i]->set(vfict[i]);
849 : }
850 1500 : } // calculate
851 :
852 : /**
853 : \brief update fictitious variables.
854 : \details This function manages evolution of fictitious variables.
855 : This function calculates mean force, updates fictitious variables by one MFD step,
856 : bounces back variables, updates free energy, and record logs.
857 : */
858 1500 : void LogMFD::update() {
859 1500 : if(getStep() == 0 || getStep()%interval != 0 ) return;
860 :
861 : // calc mean force for fictitious variables
862 15 : calcMeanForce();
863 :
864 : // update fictitious variables
865 30 : if(thermostat == "NVE") {
866 5 : updateNVE();
867 : }
868 10 : else if(thermostat == "NVT") {
869 5 : updateNVT();
870 : }
871 5 : else if(thermostat == "VS") {
872 5 : updateVS();
873 : }
874 :
875 : // bounce back variables if they are beyond their min and max.
876 45 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
877 45 : if(fict[i] > fict_max[i]) {
878 0 : fict[i] = fict_max[i] - (fict[i] - fict_max[i]);
879 0 : vfict[i] *= -1.0;
880 : }
881 30 : if(fict[i] < fict_min[i]) {
882 0 : fict[i] = fict_min[i] + (fict_min[i] - fict[i]);
883 0 : vfict[i] *= -1.0;
884 : }
885 : }
886 :
887 : // update free energy
888 15 : flog = calcFlog();
889 :
890 : // record log for fictitious variables
891 15 : if( multi_sim_comm.Get_rank()==0 && comm.Get_rank()==0 ) {
892 15 : FILE *outlog = std::fopen("logmfd.out", "a");
893 :
894 15 : const double ekin = calcEkin();
895 30 : const double temp = 2.0*ekin/getNumberOfArguments()/plumed.getAtoms().getKBoltzmann();
896 :
897 15 : fprintf(outlog, "%*d", 8, (int)getStep()/interval);
898 15 : fprintf(outlog, "%17.8f", flog);
899 : fprintf(outlog, "%17.8f", temp);
900 15 : fprintf(outlog, "%17.8f", xeta);
901 15 : fprintf(outlog, "%17.8f", veta);
902 45 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
903 30 : fprintf(outlog, "%17.8f", fict[i]);
904 15 : fprintf(outlog, "%17.8f", vfict[i]);
905 15 : fprintf(outlog, "%17.8f", ffict[i]);
906 : }
907 : fprintf(outlog," \n");
908 15 : fclose(outlog);
909 : }
910 :
911 : // record log for collective variables
912 15 : if( comm.Get_rank()==0 ) {
913 : // the number of replica is added to file name to distingwish replica.
914 15 : FILE *outlog2 = fopen("replica.out", "a");
915 15 : fprintf(outlog2, "%*d", 8, (int)getStep()/interval);
916 15 : fprintf(outlog2, "%16.6e ", work);
917 15 : fprintf(outlog2, "%16.6e ", weight);
918 45 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
919 30 : fprintf(outlog2, "%17.8f", fict_ave[i]);
920 : }
921 : fprintf(outlog2," \n");
922 15 : fclose(outlog2);
923 : }
924 :
925 : // reset mean force
926 45 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
927 30 : ffict[i] = 0.0;
928 15 : fict_ave[i] = 0.0;
929 : }
930 :
931 : } // update
932 :
933 : /**
934 : \brief update fictitious variables by NVE mechanics.
935 : \details This function updates ficitious variables by one NVE-MFD step using mean forces
936 : and flattening coefficient and free energy.
937 : */
938 5 : void LogMFD::updateNVE() {
939 5 : const double dt = delta_t;
940 :
941 : // get latest free energy and flattening coefficient
942 5 : flog = calcFlog();
943 5 : const double clog = calcClog();
944 :
945 : // update all ficitious variables by one MFD step
946 15 : for( unsigned i=0; i<getNumberOfArguments(); ++i ) {
947 : // update velocity (full step)
948 20 : vfict[i]+=clog*ffict[i]*dt/mfict[i];
949 : // update position (full step)
950 10 : fict[i]+=vfict[i]*dt;
951 : }
952 5 : } // updateNVE
953 :
954 : /**
955 : \brief update fictitious variables by NVT mechanics.
956 : \details This function updates ficitious variables by one NVT-MFD step using mean forces
957 : and flattening coefficient and free energy.
958 : */
959 5 : void LogMFD::updateNVT() {
960 5 : const double dt = delta_t;
961 5 : const double nkt = getNumberOfArguments()*kbt;
962 :
963 : // backup vfict
964 5 : std::vector<double> vfict_backup(getNumberOfArguments());
965 15 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
966 10 : vfict_backup[i] = vfict[i];
967 : }
968 :
969 : const int niter=5;
970 55 : for(unsigned j=0; j<niter; ++j) {
971 : // get latest free energy and flattening coefficient
972 25 : flog = calcFlog();
973 25 : const double clog = calcClog();
974 :
975 : // restore vfict from vfict_backup
976 75 : for(unsigned l=0; l<getNumberOfArguments(); ++l) {
977 50 : vfict[l] = vfict_backup[l];
978 : }
979 :
980 : // evolve vfict from vfict_backup by dt/2
981 75 : for(unsigned m=0; m<getNumberOfArguments(); ++m) {
982 50 : vfict[m] *= exp(-0.25*dt*veta);
983 100 : vfict[m] += 0.5*dt*clog*ffict[m]/mfict[m];
984 50 : vfict[m] *= exp(-0.25*dt*veta);
985 : }
986 : }
987 :
988 : // get latest free energy and flattening coefficient
989 5 : flog = calcFlog();
990 5 : const double clog = calcClog();
991 :
992 : // evolve vfict by dt/2, and evolve fict by dt
993 15 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
994 10 : vfict[i] *= exp(-0.25*dt*veta);
995 20 : vfict[i] += 0.5*dt*clog*ffict[i]/mfict[i];
996 10 : vfict[i] *= exp(-0.25*dt*veta);
997 10 : fict[i] += dt*vfict[i];
998 : }
999 :
1000 : // evolve xeta and veta by dt
1001 5 : xeta += 0.5*dt*veta;
1002 5 : const double ekin = calcEkin();
1003 5 : veta += dt*(2.0*ekin-nkt)/meta;
1004 5 : xeta += 0.5*dt*veta;
1005 5 : } // updateNVT
1006 :
1007 : /**
1008 : \brief update fictitious variables by VS mechanics.
1009 : \details This function updates ficitious variables by one VS-MFD step using mean forces
1010 : and flattening coefficient and free energy.
1011 : */
1012 5 : void LogMFD::updateVS() {
1013 5 : const double dt = delta_t;
1014 5 : const double nkt = getNumberOfArguments()*kbt;
1015 :
1016 : // get latest free energy and flattening coefficient
1017 5 : flog = calcFlog();
1018 5 : const double clog = calcClog();
1019 :
1020 15 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1021 : // update velocity (full step)
1022 20 : vfict[i] += clog*ffict[i]*dt/mfict[i];
1023 : }
1024 :
1025 5 : const double ekin = calcEkin();
1026 5 : const double svs = sqrt(nkt/ekin/2);
1027 :
1028 15 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1029 : // update position (full step)
1030 10 : vfict[i] *= svs;
1031 10 : fict[i] += vfict[i]*dt;
1032 : }
1033 :
1034 : // evolve VS potential
1035 5 : phivs += nkt*std::log(svs);
1036 5 : } // updateVS
1037 :
1038 : /**
1039 : \brief calculate mean force for fictitious variables.
1040 : \details This function calculates mean forces by averaging forces accumulated during one MFD step,
1041 : update work variables done by fictitious variables by one MFD step,
1042 : calculate weight variable of this replica for LogPD.
1043 : */
1044 15 : void LogMFD::calcMeanForce() {
1045 : // cale temporal mean force for each CV
1046 45 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1047 30 : ffict[i] /= interval;
1048 : // average of diff (getArgument(i)-fict[i])
1049 15 : fict_ave[i] /= interval;
1050 : // average of getArgument(i)
1051 30 : fict_ave[i] += fict[i];
1052 :
1053 : // correct fict_ave so that it is inside [min:max].
1054 45 : fict_ave[i] = fictValue[i]->bringBackInPbc(fict_ave[i]);
1055 : }
1056 :
1057 : // accumulate work, it was initialized as 0.0
1058 45 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1059 45 : work -= ffict[i] * vfict[i] * delta_t; // modified sign
1060 : }
1061 :
1062 : // for replica parallel
1063 15 : if( multi_sim_comm.Get_size()>1 ) {
1064 : // find the minimum work among all replicas
1065 0 : double work_min = work;
1066 0 : multi_sim_comm.Min(work_min);
1067 :
1068 : // weight of this replica.
1069 : // here, work is reduced by work_min to avoid all exp(-work/kbt)s disconverge
1070 0 : if( kbt == 0.0 ) {
1071 0 : weight = work==work_min ? 1.0 : 0.0;
1072 : }
1073 : else {
1074 0 : weight = exp(-(work-work_min)/kbt);
1075 : }
1076 :
1077 : // normalize the weight
1078 0 : double sum_weight = weight;
1079 0 : multi_sim_comm.Sum(sum_weight);
1080 0 : weight /= sum_weight;
1081 :
1082 : // weighting force of this replica
1083 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1084 0 : ffict[i] *= weight;
1085 : }
1086 :
1087 : // mean forces of all replica.
1088 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1089 0 : multi_sim_comm.Sum(ffict[i]);
1090 : }
1091 : // now, mean force is obtained.
1092 : }
1093 15 : } // calcMeanForce
1094 :
1095 : /**
1096 : \brief calculate kinetic energy of fictitious variables.
1097 : \retval kinetic energy.
1098 : \details This function calculates sum of kinetic energy of all fictitious variables.
1099 : */
1100 82 : double LogMFD::calcEkin() {
1101 : double ekin=0.0;
1102 246 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
1103 246 : ekin += mfict[i]*vfict[i]*vfict[i]*0.5;
1104 : }
1105 82 : return ekin;
1106 : } // calcEkin
1107 :
1108 : /**
1109 : \brief calculate free energy of fictitious variables.
1110 : \retval free energy.
1111 : \details This function calculates free energy by using invariant of canonical mechanics.
1112 : */
1113 55 : double LogMFD::calcFlog() {
1114 55 : const double nkt = getNumberOfArguments()*kbt;
1115 55 : const double ekin = calcEkin();
1116 : double pot;
1117 :
1118 110 : if (thermostat == "NVE") {
1119 10 : pot = hlog - ekin;
1120 : }
1121 45 : else if (thermostat == "NVT") {
1122 35 : const double ekin_bath = 0.5*veta*veta*meta+xeta*nkt;
1123 35 : pot = hlog - ekin - ekin_bath;
1124 : }
1125 10 : else if (thermostat == "VS") {
1126 10 : pot = phivs;
1127 : }
1128 : else {
1129 : pot = 0.0; // never occurs
1130 : }
1131 :
1132 110 : return TAMD ? pot : sgn(pot)*expm1(fabs(pot)/gamma)/alpha;
1133 : } // calcFlog
1134 :
1135 : /**
1136 : \brief calculate coefficient for flattening.
1137 : \retval flattering coefficient.
1138 : \details This function returns 1.0 for TAMD, flattening coefficient for LogMFD.
1139 : */
1140 40 : double LogMFD::calcClog() {
1141 40 : return TAMD ? 1.0 : alpha*gamma/(alpha*fabs(flog)+1.0);
1142 : } // calcClog
1143 :
1144 : } // logmfd
1145 5517 : } // PLMD
|
