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=1/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], [71] (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 [77];

\[ {\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 of \(N_m\) steps 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) [78].

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, free energy profiles can be reconstructed more efficiently (requiring less elapsed computing time) in LogPD than with a single MD system in LogMFD. In the case that there exist more than one stable state separated by high energy barriers in the hidden subspace orthogonal to the CV-subspace, LogPD is particularly of use to incorporate all the contributions from such hidden states with appropriate weights (in the limit \(N_r\to\infty\) ).

Note that LogPD calculations should always be initiated with an equilibrium \({\bf q}\)-configuration in each replica, because the Crooks-Jarzynski non-equilibrium work relation is invoked. Also note that LogPD is currently available only with Gromacs, while LogMFD can be performed with LAMMPS, Gromacs, Quantum Espresso, NAMD, and so on.

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

\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
tested on master
UNITS 
TIME
the units of time.
=fs
LENGTH
the units of lengths.
=1.0
ENERGY
the units of energy.
=kcal/mol
MASS
the units of masses.
=1.0
CHARGE
the units of charges.
=1.0 phi: TORSION
ATOMS
the four atoms involved in the torsional angle
=5,7,9,15 psi: TORSION
ATOMS
the four atoms involved in the torsional angle
=7,9,15,17 # LogMFD logmfd: LOGMFD ...
ARG
compulsory keyword the labels of the scalars on which the bias will act
=phi,psi
KAPPA
compulsory keyword Spring constant of the harmonic restraining potential.
=1000.0,1000.0
DELTA_T
compulsory keyword Time step for the fictitious dynamical variables (DELTA_T=1 often works).
=1.0
INTERVAL
compulsory keyword Period of MD steps (N_m) to update fictitious dynamical variables.
=200
TEMP
Target temperature for the fictitious dynamical variables in LogMFD/PD.
=300.0
FLOG
compulsory keyword The initial free energy value in the LogMFD/PD run.The
=2.0
MFICT
Masses of each fictitious dynamical variable.
=5000000,5000000
VFICT
The initial velocities of the fictitious dynamical variables.
=3.5e-4,3.5e-4
ALPHA
Alpha parameter for LogMFD.
=4.0
THERMOSTAT
compulsory keyword Type of thermostat for the fictitious dynamical variables.
=NVE
FICT_MAX
compulsory keyword Maximum values reachable for the fictitious dynamical variables.
=3.15,3.15
FICT_MIN
compulsory keyword Minimum values reachable for the fictitious dynamical variables.
=-3.15,-3.15 ...

To submit this simulation with Gromacs, use the following command line to execute a LogMFD run with Gromacs-MD. Here topol.tpr is the input file which contains atomic coordinates and Gromacs parameters.

gmx_mpi mdrun -s topol.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.000000 5000000.000000
# Mass for thermostat   :   11923.224809
# 1:iter_mfd, 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),
       0       2.000000     308.221983       0.000000       0.000000      -2.442363       0.000350       5.522717       2.426650       0.000350       7.443177
       1       1.995466     308.475775       0.000000       0.000000      -2.442013       0.000350      -4.406246       2.427000       0.000350      11.531345
       2       1.992970     308.615664       0.000000       0.000000      -2.441663       0.000350      -3.346513       2.427350       0.000350      15.763196
       3       1.988619     308.859946       0.000000       0.000000      -2.441313       0.000350       6.463092       2.427701       0.000351       6.975422
...

The output file replica.out records all collective variables at every MFD step.

replica.out

 Replica No. 0 of 1.
# 1:iter_mfd, 2:work, 3:weight,
# 4:phi(q)
# 5:psi(q)
       0    0.000000e+00     1.000000e+00       -2.436841       2.434093
       1   -4.539972e-03     1.000000e+00       -2.446420       2.438531
       2   -7.038516e-03     1.000000e+00       -2.445010       2.443114
       3   -1.139672e-02     1.000000e+00       -2.434850       2.434677
...

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 0/topol.tpr and 1/topol.tpr are the input files, each containing the atomic coordinates for the corresponding replica and Gromacs parameters. All the directories, 0/ and 1/, should contain the same plumed.dat.

mpirun -np 2 gmx_mpi mdrun -s topol -plumed -multidir 0 1

This command will output files named logmfd.out in 0/, and replica.out.0 and replica.out.1 in 0/ and 1/, respectively.

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.000000 5000000.000000
# Mass for thermostat   :   11923.224809
# 1:iter_mfd, 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),
       0       2.000000     308.221983       0.000000       0.000000      -2.417903       0.000350       4.930899       2.451475       0.000350      -3.122024
       1       1.999367     308.257404       0.000000       0.000000      -2.417552       0.000350       0.851133       2.451825       0.000350      -1.552718
       2       1.999612     308.243683       0.000000       0.000000      -2.417202       0.000350      -1.588175       2.452175       0.000350       1.601274
       3       1.999608     308.243922       0.000000       0.000000      -2.416852       0.000350       4.267745       2.452525       0.000350      -4.565621
...

The output file replica.out.0 records all collective variables and the Jarzynski weight of replica No.0 at every MFD step.

replica.out.0

# Replica No. 0 of 2.
# 1:iter_mfd, 2:work, 3:weight,
# 4:phi(q)
# 5:psi(q)
       0    0.000000e+00     5.000000e-01       -2.412607       2.446191
       1   -4.649101e-06     4.994723e-01       -2.421403       2.451318
       2    1.520985e-03     4.983996e-01       -2.420269       2.455049
       3    1.588855e-03     4.983392e-01       -2.411081       2.447394
...

The output file replica.out.1 records all collective variables and the Jarzynski weight of replica No.1 at every MFD step.

replica.out.1

# Replica No. 1 of 2.
# 1:iter_mfd, 2:work, 3:weight,
# 4:phi(q)
# 5:psi(q)
       0    0.000000e+00     5.000000e-01       -2.413336       2.450516
       1   -1.263077e-03     5.005277e-01       -2.412009       2.449229
       2   -2.295444e-03     5.016004e-01       -2.417322       2.452512
       3   -2.371507e-03     5.016608e-01       -2.414078       2.448521
...
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
ARG the labels of the scalars on which the bias will act
INTERVAL Period of MD steps (N_m) to update fictitious dynamical variables.
DELTA_T Time step for the fictitious dynamical variables (DELTA_T=1 often works).
THERMOSTAT Type of thermostat for the fictitious dynamical variables. NVE, NVT, VS are available.
KAPPA Spring constant of the harmonic restraining potential.
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({X}(t)) > 0 at any X(t), resulting in enhanced barrier-crossing. (The value of H_log is automatically set according to FLOG). In fact, F({X}(t)) is correctly estimated in PLUMED even when F({X}(t)) < 0 in the LogMFD/PD run.
Options
NUMERICAL_DERIVATIVES

( default=off ) calculate the derivatives for these quantities numerically

TEMP Target temperature for the fictitious dynamical variables in LogMFD/PD. It is recommended to set TEMP to be the same as the temperature of the MD system in any thermostated LogMFD/PD run. If not provided, it will be taken from the temperature of the MD system (Gromacs only).
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. If THERMOSTAT=VS, they are instead automatically adjusted according to TEMP.
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.
TEMPPD Temperature of the Boltzmann factor in the Jarzynski weight in LogPD (Gromacs only). TEMPPD should be the same as the temperature of the MD system, while TEMP may be (in principle) different from it. If not provided, TEMPPD is set to be the same value as TEMP.