LOGMFD
 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.

# LogMFD

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], [74] (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 [80];

${\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).

# LogPD

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) [81].

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.

# Using LogMFD/PD with a thermostat

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}$$.

## Nose-Hoover LogMFD/PD

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]$

## Velocity scaling LogMFD/PD

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 [82];

\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.

Examples

# Examples

## Example of LogMFD

The following input file tells plumed to restrain collective variables to the fictitious dynamical variables in LogMFD/PD.

plumed.dat

Click on the labels of the actions for more information on what each action computes
UNITS TIMEthe units of time. =fs LENGTHthe units of lengths. =1.0 ENERGYthe units of energy. =kcal/mol MASSthe units of masses. =1.0 CHARGEthe units of charges. =1.0
phi: TORSION ATOMSthe four atoms involved in the torsional angle =5,7,9,15
psi: TORSION ATOMSthe four atoms involved in the torsional angle =7,9,15,17
# LogMFD
logmfd: LOGMFD ...

ARGthe input for this action is the scalar output from one or more other actions. =phi,psi
KAPPAcompulsory keyword
Spring constant of the harmonic restraining potential for the fictitious dynamical
variables. =100.0,100.0
DELTA_Tcompulsory keyword
Time step for the fictitious dynamical variables (MFD step). =0.5
INTERVALcompulsory keyword
Period of MD steps (\f$N_m\f$) to update fictitious dynamical variables. =500
TEMPTemperature of the fictitious dynamical variables in LogMFD/PD thermostat. =300.0
FLOGcompulsory keyword
The initial free energy value in the LogMFD/PD run.The =5.0
MFICTMasses of each fictitious dynamical variable. =5000000.0,5000000.0
VFICTThe initial velocities of the fictitious dynamical variables. =3.5e-4,3.5e-4
ALPHAAlpha parameter for LogMFD. =4.0
THERMOSTATcompulsory keyword
Type of thermostat for the fictitious dynamical variables. =NVE
VETAThe initial velocity of eta variable. =0.0
METAMass of eta variable. =20000.0
FICT_MAXcompulsory keyword
Maximum values reachable for the fictitious dynamical variables. =3.1,3.1
FICT_MINcompulsory keyword
Minimum values reachable for the fictitious dynamical variables. =-3.1,-3.1
...


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
...


## Example of LogPD

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
...

Glossary of keywords and components
Description of components

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
Compulsory keywords
 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).
Options
 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.