Molecular Simulation/Replica Exchange Molecular Dynamics

Replica Exchange molecular dynamics involves performing simultaneous simulations on many ($N$) copies of a system with a different value for a certain control variable. Parallel tempering is a replica exchange method in which temperature is used as the control parameter. A set of temperatures are assigned to the $N$ copies of a system satisfying $T_{N}>T_{N-1}>\dots >T_{1}$. The probability distribution for an independent system depends on the configuration of all the atoms $(\mathbf {R} )$ within a system and $\mathbf {R} _{1},\mathbf {R} _{2},\dots ,\mathbf {R} _{N}$ represents the complete set of configurations for $N$ copies of the system. Replica exchange simulations are performed on canonical ensemble (NVT) systems. An exchange of coordinates between a neighbouring pair of replica is attempted, and the attempted move is accepted or rejected based on Metropolis Monte Carlo criterion.

A replica exchange molecular dynamics simulation of the folding of a peptide. The lowest temperature replica is labeled in bold].

A replica exchange molecular dynamics simulation where there are 4 replicas at different temperatures. m steps of molecular dynamics are completed between exchange attempts.

A schematic of a replica exchange simulation carried out on 4 replicas at different temperatures is shown. The m steps of MD are done on each individual replica between exchange attempts. The first attempted exchange between replicas 3 and 4 is made, $(\mathbf {R} _{3},\mathbf {R} _{4})\rightarrow ({\tilde {\mathbf {R} }}_{3},{\tilde {\mathbf {R} }}_{4})$, and the attempted exchange is accepted, ${\tilde {\mathbf {R} }}_{3}=\mathbf {R} _{4}$ and ${\tilde {\mathbf {R} }}_{4}=\mathbf {R} _{3}$.The attempted exchange between replicas 1 and 2 is rejected, therefore, $(\mathbf {R} _{1},\mathbf {R} _{2})\rightarrow ({\tilde {\mathbf {R} }}_{1},{\tilde {\mathbf {R} }}_{2})$ where ${\tilde {\mathbf {R} }}_{1}=\mathbf {R} _{1}$ and ${\tilde {\mathbf {R} }}_{2}=\mathbf {R} _{2}$.

Advantages

The aim of replica exchange molecular dynamics is to ensure that high temperature configurations can cross barriers easily on the potential energy surface. The low temperature replicas sample the lowest energy minima, while the high temperature replicas are capable of sampling the entire surface, making lower temperature configurations accessible at higher temperatures. Furthermore, conventional molecular dynamics or Monte Carlo methods show extremely slow convergence when calculating equilibrium properties, therefore, it takes a long time for a configuration to escape a local minimum because the probability of crossing an energy barrier of height $q^{\ddagger }$ is proportional to $\exp \left(-{\frac {q^{\ddagger }}{k_{B}T}}\right)$. Replica exchange methods improve sampling efficiency by making it easier to cross energy barriers.

The Acceptance Probability

Simulation of two NVT systems where the second simulation is thermostated at a higher temperature, $T_{2}>T_{1}$.

The probability distribution of System 1:

${\begin{aligned}P_{1}(\mathbf {R} _{1})={\frac {\exp \left(-{\mathcal {V}}(\mathbf {R} _{1})/k_{B}T_{1}\right)}{\int \exp \left(-{\mathcal {V}}(\mathbf {r} )/k_{B}T_{1}\right)d\mathbf {r} }}\\\end{aligned}}$

$\mathbf {R} _{1}$ represents the coordinates of all atoms in system 1.

The probability distribution of System 2:

${\begin{aligned}P_{2}(\mathbf {R} _{2})={\frac {\exp \left(-{\mathcal {V}}(\mathbf {R} _{2})/k_{B}T_{2}\right)}{\int \exp \left(-{\mathcal {V}}(\mathbf {r} )/k_{B}T_{2}\right)d\mathbf {r} }}\\\end{aligned}}$

$\mathbf {R} _{2}$ represents the coordinates of all atoms in system 2.

The probability of System 1 having coordinates $\mathbf {R} _{1}$ and System 2 having coordinates $\mathbf {R} _{2}$:

${\begin{aligned}P_{1}(\mathbf {R} _{1})P_{2}(\mathbf {R} _{2})={\frac {\exp \left(-{\mathcal {V}}(\mathbf {R} _{1})/k_{B}T_{1}\right)}{\int \exp \left(-{\mathcal {V}}(\mathbf {r} )/k_{B}T_{1}\right)d\mathbf {r} }}{\frac {\exp \left(-{\mathcal {V}}(\mathbf {R} _{2})/k_{B}T_{2}\right)}{\int \exp \left(-{\mathcal {V}}(\mathbf {r} )/k_{B}T_{2}\right)d\mathbf {r} }}\\\end{aligned}}$

The probability of System 1 having coordinates $\mathbf {R} _{2}$ and System 2 having coordinates $\mathbf {R} _{1}$:

${\begin{aligned}P_{1}(\mathbf {R} _{2})P_{2}(\mathbf {R} _{1})&={\frac {\exp \left(-{\mathcal {V}}(\mathbf {R} _{2})/k_{B}T_{1}\right)}{\int \exp \left(-{\mathcal {V}}(\mathbf {r} )/k_{B}T_{1}\right)d\mathbf {r} }}{\frac {\exp \left(-{\mathcal {V}}(\mathbf {R} _{1})/k_{B}T_{2}\right)}{\int \exp \left(-{\mathcal {V}}(\mathbf {r} )/k_{B}T_{2}\right)d\mathbf {r} }}\\\end{aligned}}$

The probability of accepting the attempted exchange of coordinates between System 1 and 2 becomes:

${\begin{aligned}{\frac {P_{\textit {swapped}}}{P_{\textit {old}}}}&={\frac {{\frac {\exp \left(-{\mathcal {V}}(\mathbf {R} _{2})/k_{B}T_{1}\right)}{\int \exp \left(-{\mathcal {V}}(\mathbf {r} )/k_{B}T_{1}\right)d\mathbf {r} }}{\frac {\exp \left(-{\mathcal {V}}(\mathbf {R} _{1})/k_{B}T_{2}\right)}{\int \exp \left(-{\mathcal {V}}(\mathbf {r} )/k_{B}T_{2}\right)d\mathbf {r} }}}{{\frac {\exp \left(-{\mathcal {V}}(\mathbf {R} _{1})/k_{B}T_{1}\right)}{\int \exp \left(-{\mathcal {V}}(\mathbf {r} )/k_{B}T_{1}\right)d\mathbf {r} }}{\frac {\exp \left(-{\mathcal {V}}(\mathbf {R} _{2})/k_{B}T_{2}\right)}{\int \exp \left(-{\mathcal {V}}(\mathbf {r} )/k_{B}T_{2}\right)d\mathbf {r} }}}}\\&={\frac {\exp \left(-{\mathcal {V}}\left(\mathbf {R} _{2}\right)/k_{B}T_{1}\right)}{\exp \left(-{\mathcal {V}}\left(\mathbf {R} _{1}\right)/k_{B}T_{1}\right)}}{\frac {\exp \left(-{\mathcal {V}}\left(\mathbf {R} _{1}\right)/k_{B}T_{2}\right)}{\exp \left(-{\mathcal {V}}\left(\mathbf {R} _{2}\right)/k_{B}T_{2}\right)}}\\&=\exp \left({\frac {-{\mathcal {V}}\left(\mathbf {R} _{2}\right)}{k_{B}T_{1}}}\right)\exp \left({\frac {+{\mathcal {V}}\left(\mathbf {R} _{1}\right)}{k_{B}T_{1}}}\right)\exp \left({\frac {-{\mathcal {V}}\left(\mathbf {R} _{1}\right)}{k_{B}T_{2}}}\right)\exp \left({\frac {+{\mathcal {V}}\left(\mathbf {R} _{2}\right)}{k_{B}T_{2}}}\right)\\&=\exp \left({\frac {[-{\mathcal {V}}(\mathbf {R} _{2})+{\mathcal {V}}(\mathbf {R} _{1})]}{k_{B}T_{1}}}\right)\exp \left({\frac {[-{\mathcal {V}}(\mathbf {R} _{1})+{\mathcal {V}}(\mathbf {R} _{2})]}{k_{B}T_{2}}}\right)\end{aligned}}$

Letting $\Delta {\mathcal {V}}={\mathcal {V}}(\mathbf {R} _{2})-{\mathcal {V}}(\mathbf {R} _{1})$

${\begin{aligned}{\frac {P_{\textit {swapped}}}{P_{\textit {old}}}}&=\exp \left({\frac {-\Delta {\mathcal {V}}}{k_{B}T_{1}}}\right)\exp \left({\frac {+\Delta {\mathcal {V}}}{k_{B}T_{2}}}\right)\\&=\exp \left({\frac {-\Delta {\mathcal {V}}}{k_{B}T_{1}}}+{\frac {\Delta {\mathcal {V}}}{k_{B}T_{2}}}\right)\\&=\exp \left(-\Delta {\mathcal {V}}\beta _{1}+\Delta {\mathcal {V}}\beta _{2}\right)\\&=\exp \left(-\Delta {\mathcal {V}}\left[\beta _{1}-\beta _{2}\right]\right)\\\end{aligned}}$

Where $\beta _{1}={\frac {1}{k_{B}T_{1}}}$ and $\beta _{2}={\frac {1}{k_{B}T_{2}}}$.

${\begin{aligned}P_{\textit {acceptance}}={\text{min}}\left[1,\exp \left(-\Delta {\mathcal {V}}\left({\frac {1}{k_{B}T_{1}}}-{\frac {1}{k_{B}T_{2}}}\right)\right)\right]\\\end{aligned}}$