 |
 |
 |
| |
|
|
| This is part of the logmfd module | |
| It is only available if you configure PLUMED with ./configure –enable-modules=logmfd . Furthermore, this feature is still being developed so take care when using it and report any problems on the mailing list. |
Used to perform LogMFD, LogPD, and TAMD/d-AFED.
Consider a physical system of \(N_q\) particles, for which the Hamiltonian is given as
\[ {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) \]
where \({\bf q}_j\), \({\bf p}_j\) ( \(\bf{\Gamma}={\bf q},{\bf p}\)), and \(m_j\) are the position, momentum, and mass of particle \(j\) respectively, and \(\Phi\) is the potential energy function for \({\bf q}\). The free energy \(F({\bf X})\) as a function of a set of \(N\) collective variables (CVs) is given as
\begin{eqnarray*} 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 }}} \\ &\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 }} \end{eqnarray*}
where \(\bf{s}\) are the CVs, \(k_B\) is Boltzmann constant, \(\beta=k_BT\), and \(K_i\) is the spring constant which is large enough to invoke
\[ \delta \left( x \right) = \lim_{k \to \infty } \sqrt {\beta k/2\pi} \exp \left( -\beta kx^2/2 \right) \]
In mean-force dynamics, \({\bf X}\) are treated as fictitious dynamical variables, which are associated with the following Hamiltonian,
\[ {H_{\rm log}} = \sum\limits_{i = 1}^N {\frac{{P_{{X_i}}^2}}{{2{M_i}}} + \Psi \left( {{\bf X}} \right)} \]
where \({P_{{X_i}}}\) and \(M_i\) are the momentum and mass of \(X_i\), respectively, and \(\Psi\) is the potential function for \({\bf X}\). We assume that \(\Psi\) is a functional of \(F({\bf X})\); \(Ψ[F({\bf X})]\). The simplest form of \(\Psi\) is \(\Psi = F({\bf X})\), which corresponds to TAMD/d-AFED [1], [66] (or the extended Lagrangian dynamics in the limit of the adiabatic decoupling between \(\bf{q}\) and \({\bf X}\)). In the case of LogMFD, the following form of \(\Psi\) is introduced [71];
\[ {\Psi _{\rm log}}\left( {{\bf X}} \right) = \gamma \log \left[ {\alpha F\left( {{\bf X}} \right) + 1} \right] \]
where \(\alpha\) (ALPHA) and \(\gamma\) (GAMMA) are positive parameters. The logarithmic form of \(\Psi_{\rm log}\) ensures the dynamics of \({\bf X}\) on a much smoother energy surface [i.e., \(\Psi_{\rm log}({\bf X})\)] than \(F({\bf X})\), thus enhancing the sampling in the \({\bf X}\)-space. The parameters \(\alpha\) and \(\gamma\) determine the degree of flatness of \(\Psi_{\rm log}\), but adjusting only \(\alpha\) is normally sufficient to have a relatively flat surface (with keeping the relation \(\gamma=1/\alpha\)).
The equation of motion for \(X_i\) in LogMFD (no thermostat) is
\[ {M_i}{\ddot X_i} = - \left( {\frac{{\alpha \gamma }}{{\alpha F + 1}}} \right)\frac{{\partial F}}{{\partial {X_i}}} \]
where \(-\partial F/\partial X_i\) is evaluated as a canonical average under the condition that \({\bf X}\) is fixed;
\begin{eqnarray*} - \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 }}\\ &\equiv& {\left\langle {{K_i}\left( {{s_i}\left( {{\bf q}} \right) - {X_i}} \right)} \right\rangle _{{\bf X}}} \end{eqnarray*}
where
\[ 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 }} \]
The mean-force (MF) is practically evaluated by performing a shot-time canonical MD run each time \({\bf X}\) is updated according to the equation of motion for \({\bf X}\).
If the canonical average for the MF is effectively converged, the dynamical variables \({\bf q}\) and \({\bf X}\) are decoupled and they evolve adiabatically, which can be exploited for the on-the-fly evaluation of \(F({\bf X})\). I.e., \(H_{\rm log}\) should be a constant of motion in this case, thus \(F({\bf X})\) can be evaluated each time \({\bf X}\) is updated as
\[ F\left( {{{\bf X}}\left( t \right)} \right) = \frac{1}{\alpha} \left[ \exp \frac{1}{\gamma} \left\{ \left( H_{\rm log} - \sum_i \frac{P_{X_i}^2}{2M_i} \right) \right\} - 1 \right] \]
This means that \(F({\bf X})\) can be constructed without post processing (on-the-fly free energy reconstruction). Note that the on-the-fly free energy reconstruction is also possible in TAMD/d-AFED if the Hamiltonian-like conserved quantity is available (e.g., the Nose-Hoover type dynamics).
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) [72].
In PD, the MF is evaluated by a non-equilibrium path-ensemble based on the Crooks-Jarzynski non-equilibrium work relation. To this end, multiple replicas of the MD system which run in parallel are introduced. The CVs [ \({\bf s}({\bf q})\)] in each replica is restrained to the same value of \({\bf X}(t)\). A canonical MD run with \(N_m\) steps is performed in each replica, then the MF on \(X_i\) is evaluated using the MD trajectories from all replicas. The MF is practically calculated as
\[ - \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]} \]
where \({\bf q}^k_n\) indicates the \({\bf q}\)-configuration at time step \(n\) in the \(N_m\)-step shot-time MD run in the \(k\)th replica among \(N_r\) replicas. \(W_k\) is given as
\[ {W_k} = \frac{{{e^{ - \beta {w_k}\left( t \right)}}}}{{\sum\limits_{k=1}^{{N_r}} {{e^{ - \beta {w_k}\left( t \right)}}} }} \]
where
\[ {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}} } \]
\[ 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}} \]
and \(s^k_i\) is the \(i\)th CV in the \(k\)th replica.
\(W_k\) comes from the Crooks-Jarzynski non-equilibrium work relation by which we can evaluate an equilibrium ensemble average from a set of non-equilibrium trajectories. Note that, to avoid possible numerical errors in the exponential function, the following form of \(W_k\) is instead used in PLUMED,
\[ {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]} }} \]
where
\[ {w_{\min }}\left( t \right) = {\rm Min}_k\left[ {{w_k}\left( t \right)} \right] \]
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 \(N_r\to\infty\) ).
Note that LogPD calculations should always be initiated with an equilibrium \({\bf q}\)-configuration in each replica, because the Crooks-Jarzynski non-equilibrium work relation is invoked. Also note that LogPD is currently available only with Gromacs, while LogMFD can be performed with LAMMPS, Gromacs, and NAMD.
Introducing a thermostat on \({\bf X}\) is often recommended in LogMFD/PD to maintain the adiabatic decoupling between \({\bf q}\) and \({\bf X}\). In the framework of the LogMFD approach, the Nose-Hoover type thermostat and the Gaussian isokinetic (velocity scaling) thermostat can be used to control the kinetic energy of \({\bf X}\).
The equation of motion for \(X_i\) coupled to a Nose-Hoover thermostat variable \(\eta\) (single heat bath) is
\[ {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 \]
The equation of motion for \(\eta\) is
\[ Q\ddot \eta = \sum\limits_{i = 1}^N {\frac{{P_{{X_i}}^2}}{{{M_i}}} - N{k_B}T} \]
where \(N\) is the number of the CVs. Since the following pseudo-Hamiltonian is a constant of motion in Nose-Hoover LogMFD/PD,
\[ 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} \]
\(F({\bf X}(t))\) is obtained at each MFD step as
\[ 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] \]
The velocity scaling algorithm (which is equivalent to the Gaussian isokinetic dynamics in the limit \(\Delta t\to 0\)) can also be employed to control the velocity of \({\bf X}\), \(\bf{V}_x\).
The following algorithm is introduced to perform isokinetic LogMFD calculations [73];
\begin{eqnarray*} {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]\\ 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)} }}} \\ {V_{{X_i}}}^\prime \left( {{t_{n + 1}}} \right) &=& S\left( {{t_{n + 1}}} \right){V_{{X_i}}}\left( {{t_{n + 1}}} \right)\\ {X_i}\left( {{t_{n + 1}}} \right) &=& {X_i}\left( {{t_n}} \right) + \Delta t V_{X_i}^\prime \left( {{t_{n + 1}}} \right)\\ {\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)\\ F\left( {{t_{n + 1}}} \right) &=& \frac{1}{\alpha} \left[ \exp \left\{ \Psi_{\rm log} \left( t_{n+1} \right) / \gamma \right\} - 1 \right] \end{eqnarray*}
Note that \(V_{X_i}^\prime\left( {{t_0}} \right)\) is assumed to be initially given, which meets the following relation,
\[ \sum\limits_{i = 1}^N M_i V_{X_i}^{\prime 2} \left( t_0 \right) = N{k_B}{T} \]
It should be stressed that, in the same way as in the NVE and Nose-Hoover LogMFD/PD, \(F({\bf X}(t))\) can be evaluated at each MFD step (on-the-fly free energy reconstruction) in Velocity Scaling LogMFD/PD.
By default this Action calculates the following quantities. These quantities can be referenced elsewhere in the input by using this Action's label followed by a dot and the name of the quantity required from the list below.
| Quantity | Description |
| bias | the instantaneous value of the bias potential |
| _fict | For example, the fictitious collective variable for LogMFD is specified as ARG=dist12 and LABEL=logmfd in LOGMFD section in Plumed input file, the associated fictitious dynamical variable can be specified as PRINT ARG=dist12,logmfd.dist12_fict FILE=COLVAR |
| _vfict | For example, the fictitious collective variable for LogMFD is specified as ARG=dist12 and LABEL=logmfd in LOGMFD section in Plumed input file, the velocity of the associated fictitious dynamical variable can be specified as PRINT ARG=dist12,logmfd.dist12_vfict FILE=COLVAR |
| INTERVAL | Period of MD steps ( \(N_m\)) to update fictitious dynamical variables. |
| DELTA_T | Time step for the fictitious dynamical variables (MFD step). |
| THERMOSTAT | Type of thermostat for the fictitious dynamical variables. NVE, NVT, VS are available. |
| KAPPA | Spring constant of the harmonic restraining potential for the fictitious dynamical variables. |
| FICT_MAX | Maximum values reachable for the fictitious dynamical variables. The variables will elastically bounce back at the boundary (mirror boundary). |
| FICT_MIN | Minimum values reachable for the fictitious dynamical variables. The variables will elastically bounce back at the boundary (mirror boundary). |
| FLOG | The initial free energy value in the LogMFD/PD run.The origin of the free energy profile is adjusted by FLOG to realize \(F({\bf X}(t)) > 0\) at any \({\bf X}(t)\), resulting in enhanced barrier-crossing. (The value of \(H_{\rm log}\) is automatically set according to FLOG). |
| NUMERICAL_DERIVATIVES | ( default=off ) calculate the derivatives for these quantities numerically |
| ARG | the input for this action is the scalar output from one or more other actions. The particular scalars that you will use are referenced using the label of the action. If the label appears on its own then it is assumed that the Action calculates a single scalar value. The value of this scalar is thus used as the input to this new action. If * or *.* appears the scalars calculated by all the proceeding actions in the input file are taken. Some actions have multi-component outputs and each component of the output has a specific label. For example a DISTANCE action labelled dist may have three components x, y and z. To take just the x component you should use dist.x, if you wish to take all three components then use dist.*.More information on the referencing of Actions can be found in the section of the manual on the PLUMED Getting Started. Scalar values can also be referenced using POSIX regular expressions as detailed in the section on Regular Expressions. To use this feature you you must compile PLUMED with the appropriate flag. You can use multiple instances of this keyword i.e. ARG1, ARG2, ARG3... |
| TEMP | Temperature of the fictitious dynamical variables in LogMFD/PD thermostat. If not provided or provided as 0, it will be taken from the temperature of the MD system. |
| TAMD | When TAMD=1, TAMD/d-AFED calculations can be performed instead of LogMFD. Otherwise, the LogMFD protocol is switched on (default). |
| ALPHA | Alpha parameter for LogMFD. If not provided or provided as 0, it will be taken as 1/gamma. If gamma is also not provided, Alpha is set as 4, which is a sensible value when the unit of kcal/mol is used. |
| GAMMA | Gamma parameter for LogMFD. If not provided or provided as 0, it will be taken as 1/alpha. 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. |
| FICT | The initial values of the fictitious dynamical variables. If not provided, they are set equal to their corresponding CVs for the initial atomic configuration. |
| VFICT | The initial velocities of the fictitious dynamical variables. If not provided, they will be taken as 0. |
| MFICT | Masses of each fictitious dynamical variable. If not provided, they will be taken as 10000. |
| XETA | The initial eta variable of the Nose-Hoover thermostat for the fictitious dynamical variables. If not provided, it will be taken as 0. |
| VETA | The initial velocity of eta variable. If not provided, it will be taken as 0. |
| META | Mass of eta variable. If not provided, it will be taken as \(N*kb*T*100*100\). |
| WORK | The initial value of the work done by fictitious dynamical variables in each replica. If not provided, it will be taken as 0. |
The following input file tells plumed to restrain collective variables to the fictitious dynamical variables in LogMFD/PD.
plumed.dat
UNITS TIME=fs LENGTH=1.0 ENERGY=kcal/mol MASS=1.0 CHARGE=1.0 phi: TORSION ATOMS=5,7,9,15 psi: TORSION ATOMS=7,9,15,17 # LogMFD LOGMFD ... LABEL=logmfd ARG=phi,psi KAPPA=100.0,100.0 DELTA_T=0.5 INTERVAL=500 TEMP=300.0 FLOG=5.0 MFICT=5000000.0,5000000.0 # FICT=0.8,0.8 VFICT=3.5e-4,3.5e-4 ALPHA=4.0 THERMOSTAT=NVE VETA=0.0 META=20000.0 FICT_MAX=3.1,3.1 FICT_MIN=-3.1,-3.1 ... LOGMFD
To submit this simulation with Gromacs, use the following command line to execute a LogMFD run with Gromacs-MD. Here TOPO/topol0.tpr is an input file which contains atomic coordinates and Gromacs parameters.
gmx_mpi mdrun -s TOPO/topol0.tpr -plumed
This command will output files named logmfd.out and replica.out.
The output file logmfd.out records free energy and all fictitious dynamical variables at every MFD step.
logmfd.out
# LogMFD
# CVs : phi psi
# Mass for CV particles : 5000000.000000000 5000000.000000000
# Mass for thermostat : 20000.000000000
# 1:iter_md, 2:Flog, 3:2*Ekin/gkb[K], 4:eta, 5:Veta,
# 6:phi_fict(t), 7:phi_vfict(t), 8:phi_force(t),
# 9:psi_fict(t), 10:psi_vfict(t), 11:psi_force(t),
1 4.99918574 308.24149708 0.00000000 0.00000000 -2.85605938 0.00035002 5.19074544 2.79216364 0.00035000 -0.53762989
2 4.99836196 308.26124159 0.00000000 0.00000000 -2.85588436 0.00035005 4.71247605 2.79233863 0.00035000 -0.00532474
3 4.99743572 308.28344595 0.00000000 0.00000000 -2.85570932 0.00035007 5.34358230 2.79251363 0.00035000 -0.05119816
...
The output file replica.out records all collective variables at every MFD step.
replica.out
# Replica No. 0 of 1.
# 1:iter_md, 2:work, 3:weight,
# 4:phi(q)
# 5:psi(q)
1 -8.142952e-04 1.000000e+00 -2.80432694 2.78661234
2 -1.638105e-03 1.000000e+00 -2.80893462 2.79211039
3 -2.564398e-03 1.000000e+00 -2.80244854 2.79182665
...
Use the following command line to execute a LogPD run using two MD replicas (note that only Gromacs is currently available for LogPD). Here TOPO/topol0.tpr and TOPO/topol1.tpr are input files which contain atomic coordinates of each replica and Gromacs parameters.
mpirun -np 2 gmx_mpi mdrun -s TOPO/topol -plumed -multi 2
This command will output files named logmfd.out, replica.out.0 and replica.out.1.
The output file logmfd.out records free energy and all fictitious dynamical variables at every MFD step.
logmfd.out
# LogPD, replica parallel of LogMFD
# number of replica : 2
# CVs : phi psi
# Mass for CV particles : 5000000.000000000 5000000.000000000
# Mass for thermostat : 20000.000000000
# 1:iter_md, 2:Flog, 3:2*Ekin/gkb[K], 4:eta, 5:Veta,
# 6:phi_fict(t), 7:phi_vfict(t), 8:phi_force(t),
# 9:psi_fict(t), 10:psi_vfict(t), 11:psi_force(t),
1 5.00224715 308.16814691 0.00000000 0.00000000 -0.95937173 0.00034994 -12.91277494 0.78923967 0.00035000 0.07353010
2 5.00476934 308.10774854 0.00000000 0.00000000 -0.95919679 0.00034989 -11.20093553 0.78941467 0.00034999 -3.21098229
3 5.00702463 308.05376594 0.00000000 0.00000000 -0.95902187 0.00034983 -10.81712171 0.78958965 0.00034998 -2.07196718
...
The output file replica.out.0 records all collective variables of replica No.0 at every MFD step.
replica.out.0
# Replica No. 0 of 2.
# 1:iter_md, 2:work, 3:weight,
# 4:phi(q)
# 5:psi(q)
1 1.843110e-03 5.003389e-01 -1.10929125 0.83348865
2 3.466179e-03 5.010942e-01 -1.05020764 0.78731283
3 4.927870e-03 5.017619e-01 -1.04968867 0.79635198
...
The output file replica.out.1 records all collective variables of replica No.1 at every MFD step.
replica.out.1
# Replica No. 1 of 2.
# 1:iter_md, 2:work, 3:weight,
# 4:phi(q)
# 5:psi(q)
1 2.651173e-03 4.996611e-01 -1.06802968 0.74605205
2 6.075530e-03 4.989058e-01 -1.09264741 0.72681448
3 9.129358e-03 4.982381e-01 -1.08517238 0.74084241
...
1.8.14