Metadynamics

Metadynamics (MTD; also abbreviated as METAD or MetaD) is a computer simulation method in computational physics, chemistry and biology. It is used to compute free energy and other state functions of a system, where ergodicity is hindered by the form of the system's energy landscape. It was first suggested by Alessandro Laio and Michele Parrinello in 2002 and is usually applied within molecular dynamics simulations. MTD closely resembles a number of recent methods such as adaptively biased molecular dynamics, adaptive reaction coordinate forces and local elevation umbrella sampling. More recently, both the original and well-tempered metadynamics were derived in the context of importance sampling and shown to be a special case of the adaptive biasing potential setting. MTD is related to the Wang-Landau sampling.

Introduction

The technique builds on a large number of related methods including (in a chronological order) the deflation, tunneling, tabu search, local elevation, conformational flooding, Engkvist-Karlström and adaptive biasing force methods.

Metadynamics has been informally described as "filling the free energy wells with computational sand". The algorithm assumes that the system can be described by a few collective variables. During the simulation, the location of the system in the space determined by the collective variables is calculated and a positive Gaussian potential is added to the real energy landscape of the system. In this way the system is discouraged to come back to the previous point. During the evolution of the simulation, more and more Gaussians sum up, thus discouraging more and more the system to go back to its previous steps, until the system explores the full energy landscape -at this point the modified free energy becomes a constant as a function of the collective variables which is the reason for the collective variables to start fluctuating heavily. At this point the energy landscape can be recovered as the opposite of the sum of all Gaussians.

The time interval between the addition of two Gaussian functions, as well as the Gaussian height and Gaussian width, are tuned to optimize the ratio between accuracy and computational cost. By simply changing the size of the Gaussian, metadynamics can be fitted to yield very quickly a rough map of the energy landscape by using large Gaussians, or can be used for a finer grained description by using smaller Gaussians. Usually, the well-temperated metadynamics is used to change the Gaussian size adaptively. Also, the Gaussian width can be adapted with the adaptive Gaussian metadynamics.

Metadynamics has the advantage, upon methods like adaptive umbrella sampling, of not requiring an initial estimate of the energy landscape to explore. However, it is not trivial to choose proper collective variables for a complex simulation. Typically, it requires several trials to find a good set of collective variables, but there are several automatic procedure proposed: essential coordinates, Sketch-Map, and non-linear data-driven collective variables.

Multi-replica approach

Independent metadynamics simulations (replicas) can be coupled together to improve usability and parallel performance. There are several such methods proposed: the multiple walker MTD, the parallel tempering MTD, the bias-exchange MTD, and the collective-variable tempering MTD. The last three are similar to the parallel tempering method and use replica exchanges to improve sampling. Typically, the Metropolis–Hastings algorithm is used for replica exchanges, but the infinite swapping and Suwa-Todo algorithms give better replica exchange rates.

Algorithm

Assume, we have a classical ${\textstyle N}$-particle system with positions at ${\textstyle \{ \vec r_i \}}$ ${\textstyle (i \in 1...N)}$ in the Cartesian coordinates ${\textstyle (\vec r_i \in \mathbb{R}^3)}$. The particle interaction are described with a potential function ${\textstyle V \equiv V(\{ \vec r_i \})}$. The potential function form (e.g. two local minima separated by a high-energy barrier) prevents an ergodic sampling with molecular dynamics or Monte Carlo methods.

Original metadynamics

A general idea of MTD is to enhance the system sampling by discouraging revisiting of sampled states. It is achieved by augmenting the system Hamiltonian ${\textstyle H}$ with a bias potential $V_\text{bias}$:

$H = T + V + V_\text{bias}$.

The bias potential is a function of collective variables ${\textstyle (V_\text{bias} \equiv V_\text{bias}(\vec s\,))}$. A collective variable is a function of the particle positions $(\vec s \equiv \vec s(\{ \vec r_i\}))$. The bias potential is continuously updated by adding bias at rate $\omega$, where $\vec s_t$ is an instantaneous collective variable value at time $t$:

$\frac{\partial V_\text{bias}(\vec s\,)}{\partial t} = \omega\, \delta(|\vec s - \vec s_t|)$.

At infinitely long simulation time $t_\text{sim}$, the accumulated bias potential converges to free energy with opposite sing (and irrelevant constant $C$):

$V_\text{bias}(\vec s\,) = \!\!\int_0^{t_\text{sim}} \!\!\!\omega\, \delta(|\vec s - \vec s_t|)\; dt \quad\Rightarrow\quad F(\vec s\,) = -\!\!\!\!\lim_{t_\text{sim} \to \infty}\!\! V_\text{bias}(\vec s\,) + C$

For a computationally efficient implementation, the update process is discretised into $\tau$ time intervals ($\lfloor\;\rfloor$ denotes the floor function) and $\delta$-function is replaced with a localized positive kernel function $K$. The bias potential becomes a sum of the kernel functions centred at the instantaneous collective variable values $\vec s_j$ at time $\tau j$:

$V_\text{bias}(\vec s\,) \approx \tau \!\!\!\sum_{j=0}^{\left\lfloor \frac{t_\text{sim}}{\tau} \right\rfloor}\!\! \omega\, K(|\vec s - \vec s_j|)$.

Typically, the kernel is a multi-dimensional Gaussian function, which covariance matrix has diagonal non-zero elements only:

$V_\text{bias}(\vec s\,) \approx \tau \!\!\!\sum_{j=0}^{\left\lfloor \frac{t_\text{sim}}{\tau} \right\rfloor}\!\! \omega \exp\!\! \left(\! -\frac{1}{2} \left| \frac{\vec s - \vec s_j}{\vec \sigma} \right|^2 \right)$.

The parameter $\tau$, $\omega$, and $\vec \sigma$ are determined a priori and kept constant during the simulation.

Implementation

Below there is a pseudocode of MTD base on molecular dynamics (MD), where $\{\vec r\}$ and $\{\vec v\}$ are the $N$-particle system positions and velocities, respectively. The bias $V_\text{bias}$ is updated every $n = \tau/\Delta t$ MD steps, and its contribution to the system forces $\{\vec F\,\}$ is $\{\vec F_\text{bias}\}$.

set initial $\{\vec r\}$ and $\{\vec v\}$ 
set $V_\text{bias}(\vec s\,) := 0$

every MD step:
    compute CV values:
        $\vec s_t := \vec s(\{\vec r\})$
    
    every $n$ MD steps:
        update bias potential:
            $V_\text{bias}(\vec s\,) := V_\text{bias}(\vec s\,)
             + \tau \omega \exp\!\! \left(\! -\frac{1}{2} \left| \frac{\vec s - \vec s_t}{\vec \sigma} \right|^2 \right)$
    
    compute atomic forces:
        $\vec F_i := -\frac{\partial V(\{\vec r\,\})}{\partial \vec r_i}
         \overbrace{\left. -\frac{\partial V_\text{bias}(\vec s\,)}{\partial\vec s} \right|_{\vec s_t}\!\!\!
                     \frac{\partial \vec s(\{\vec r\,\})}{\partial \vec r_i}}^{\vec F_{\text{bias},i}}$
    
    propagate $\{\vec r\}$ and $\{\vec v\}$ by $\Delta t$

Applications

Metadynamics has been used to study:

• protein folding
• chemical reactions
• molecular docking
• phase transitions.

Implementations

PLUMED

PLUMED is an open-source library implementing many MTD algorithms and collective variables. It has a flexible object-oriented design and can interfaced with several MD programs (AMBER, GROMACS, LAMMPS, NAMD, Quantum ESPRESSO, and CP2K).

Other

Other MTD implementations exist in LAMMPS, NAMD, ORAC, CP2K, and Desmond.

External links

• Introduction to metadynamics
• PLUMED
• Colvars module website (NAMD and LAMMPS)
• Visual movie of metadynamics

See also

• Local Elevation
• Parallel tempering
• Umbrella sampling