Metadynamics

Independent metadynamics simulations (replicas) can be coupled together to improve usability and parallel performance. There are several such methods proposed: the multiple walker MTD,[20] the parallel tempering MTD,[21] the bias-exchange MTD,[22] and the collective-variable tempering MTD.[23] 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[24] and Suwa-Todo[25] algorithms give better replica exchange rates.[26]

Typical (single-replica) MTD simulations can include up to 3 CVs, even using the multi-replica approach, it is hard to exceed 8 CVs in practice. This limitation comes from the bias potential, constructed by adding Gaussian functions (kernels). It is a special case of the kernel density estimator (KDE). The number of required kernels, for a constant KDE accuracy, increases exponentially with the number of dimensions. So MTD simulation length has to increase exponentially with the number of CVs to maintain the same accuracy of the bias potential. Also, the bias potential, for fast evaluation, is typically approximated with a regular grid.[27] The required memory to store the grid increases exponentially with the number of dimensions (CVs) too.

A high-dimensional generalization of metadynamics is NN2B.[28] It is based on two machine learning algorithms: the nearest-neighbor density estimator (NNDE) and the artificial neural network (ANN). NNDE replaces KDE to estimate the updates of bias potential from short biased simulations, while ANN is used to approximate the resulting bias potential. ANN is a memory-efficient representation of high-dimensional functions, where derivatives (biasing forces) are effectively computed with the backpropagation algorithm.[28][29]

An alternative method, exploiting ANN for the adaptive bias potential, uses mean potential forces for the estimation.[30] This method is also a high-dimensional generalization of the Adaptive Biasing Force (ABF) method.[31] Additionally, the training of ANN is improved using Bayesian regularization,[32] and the error of approximation can be inferred by training an ensemble of ANNs.[30]

In 2015, White, Dama, and Voth introduced experiment-directed metadynamics, a method that allows for shaping molecular dynamics simulations to match a desired free energy surface. This technique guides the simulation towards conformations that align with experimental data, enhancing our understanding of complex molecular systems and their behavior.[33]

In 2020, an evolution of metadynamics was proposed, the on-the-fly probability enhanced sampling method (OPES),[34][35][36] which is now the method of choice of Michele Parrinello's research group.[37] The OPES method has only a few robust parameters, converges faster than metadynamics, and has a straightforward reweighting scheme.[38] OPES has been implemented in the PLUMED library since version 2.7.[39]

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 ${\displaystyle V_{\text{bias}}}$:

${\displaystyle 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 ${\displaystyle ({\vec {s}}\equiv {\vec {s}}(\{{\vec {r}}_{i}\}))}$. The bias potential is continuously updated by adding bias at rate ${\displaystyle \omega }$, where ${\displaystyle {\vec {s}}_{t}}$ is an instantaneous collective variable value at time ${\displaystyle t}$:

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

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

${\displaystyle 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 ${\displaystyle \tau }$ time intervals (${\displaystyle \lfloor \;\rfloor }$ denotes the floor function) and ${\displaystyle \delta }$-function is replaced with a localized positive kernel function ${\displaystyle K}$. The bias potential becomes a sum of the kernel functions centred at the instantaneous collective variable values ${\displaystyle {\vec {s}}_{j}}$ at time ${\displaystyle \tau j}$:

${\displaystyle 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, whose covariance matrix has diagonal non-zero elements only:

${\displaystyle 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 ${\displaystyle \tau }$, ${\displaystyle \omega }$, and ${\displaystyle {\vec {\sigma }}}$ are determined a priori and kept constant during the simulation.

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

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

every MD step:
    compute CV values:
        ${\displaystyle {\vec {s}}_{t}:={\vec {s}}(\{{\vec {r}}\})}$
    
    every ${\displaystyle n}$ MD steps:
        update bias potential:
            ${\displaystyle 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:
        ${\displaystyle {\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 ${\displaystyle \{{\vec {r}}\}}$ and ${\displaystyle \{{\vec {v}}\}}$ by ${\displaystyle \Delta t}$

The finite size of the kernel makes the bias potential to fluctuate around a mean value. A converged free energy can be obtained by averaging the bias potential. The averaging is started from ${\displaystyle t_{\text{diff}}}$, when the motion along the collective variable becomes diffusive:

${\displaystyle {\bar {F}}({\vec {s}})=-{\frac {1}{t_{\text{sim}}-t_{\text{diff}}}}\int _{t_{\text{diff}}}^{t_{\text{sim}}}\!\!\!\!\!V_{\text{bias}}({\vec {s}},t)\,dt+C}$

PLUMED[46] is an open-source library implementing many MTD algorithms and collective variables. It has a flexible object-oriented design[47][48] and can be interfaced with several MD programs (AMBER, GROMACS, LAMMPS, NAMD, Quantum ESPRESSO, DL_POLY_4, CP2K, and OpenMM).[49][50]

Other MTD implementations exist in the Collective Variables Module [51] (for LAMMPS, NAMD, and GROMACS), ORAC, CP2K,[52] EDM,[53] and Desmond.