Methods to handle non-ergodic systems with Molecular Dynamics

Question

Given the large difference between simulation timescales and the timescales on which we normally interact with ensembles of molecules, the time average of a molecular dynamics simulation does not always match what we would expect to observe experimentally. Effectively, there is an ergodicity problem as the molecular dynamics simulation is not exploring the entirety of the relevant phase space.

What techniques exist to improve the handling of non-ergodic systems for molecular dynamics simulations?

Answer 1

Metadynamics

One way to make your dynamics explore more of phase space is to force it out of regions it has already explored. Metadynamics does this by intermittently adding bias potentials (typically Gaussians) to the system Hamiltonian. These bias potentials are functions of user defined collective variables, which should be chosen such that states of interest can be clearly distinguished from each other. The figure below gives a good example of how metadynamics proceeds; a common analogy is that you are exploring a potential energy surface that is slowly being filled with sand/water around the current state. In the long time limit, the bias potentials should sum up to the free energy, but with the opposite sign.

The main difficultly of this method is determining what collective variables need to be considered in order to fully describe your system. Its an area of active of research to develop fairly general collective variables (or at least guidelines for their selection) for a variety of systems. As an example of what these collective variables can look like, consider [1], a paper by Michele Parrinello, who originally proposed the concept of metadynamics. In this paper, they define an enthalpic and an entropic collective variables to model crystal formation of aluminum and sodium.

1. Pablo M. Piaggi, Omar Valsson, and Michele Parrinello PRL 119, 015701 (2017)

Answer 2

Replica Exchange

Methods like temperature replica exchange and Hamiltonian replica exchange can help to sample phase space more uniformly. For temperature replica exchange, we have the base temperature that we are interested in (often room temperature or body temperature) and a number of other replicas with higher temperatures. For instance, we could have 10 replicas with temperatures of 300, 350, 400, 450, 500, 550, 600, 650, 700, and 750 K. We run simulations at these temperatures simultaneously. At the higher temperatures, enthalpic barriers easier to cross (they are smaller relative to thermal energy, $k_\mathrm{B}T$) and the kinetics are faster. For example, while the cis and trans amide conformations of proline have similar energy, there is a large barrier between (>10 kcal/mol) them. At 300 K, the transitions between the cis and trans will never occur on the simulation timescale (microseconds). You could probably run continuously for a year and never see such a transition. However, transitions can be observed within a few nanoseconds at 750 K. So, our high temperature simulations will explore more of conformational space. We can take advantage of this by applying a Monte Carlo procedure that improves sampling at the base temperature (300 K), while maintaining correct $NVT$ ensemble at this temperature.

For example, we can let the simulations run for each replica and every 1000 steps attempt exchanges of atomic coordinates between the replicas. We calculate the probability of the exchange by:

$$ p_\mathrm{exchange} = \min \left\{ 1,\frac{\exp(- \beta_i V_j - \beta_j V_i)}{\exp(-\beta_i V_i - \beta_j V_j)} \right\},\tag{1}$$

where $i$ and $j$ are the indices of two replicas, $V_i$ is the potential energy of replica $i$ before the exchange attempt, and $\beta_i = 1/(k_\mathrm{B} T_i)$ is the inverse temperature of the ensemble for replica $i$. We can see $\frac{\exp(- \beta_i V_j - \beta_j V_i)}{\exp(-\beta_i V_i - \beta_j V_j)}$, as the ratio of the probability of observing the states of the two replicas after exchange to the probability of observing the states of the two replicas before exchange. If this ratio is larger than one, then we accept the exchange immediately. If this ratio is less than one, we generate a random number on [0,1] and accept the exchange if that number is less than the probability ratio. If the atomic configurations are exchanged, we reset the velocities from a Maxwell-Boltzmann distribution at the destination temperature (alternatively, we can include the kinetic energy in the probability calculation and exchange velocities too).

Hence, barriers can be crossed at high temperature and the resulting atomic configurations can be passed down to lower temperatures. The reason that we often need many temperatures is that, if the temperatures are too far apart, the energies distributions will not overlap, $p$ will be very small, and the exchanges will rarely, if ever, be accepted. We need exchanges to be accepted to get the enhanced sampling we desire. So the temperature range and number of replicas must be chosen carefully to both obtain better exploration of phase space (usually promoted by a higher maximum temperature) and a reasonably high acceptance rate (promoted by having a small difference in temperature between adjacent replicas).

A major problem with the temperature replica exchange scheme is that it does not scale well as the number of degrees of freedom in the system is increased. The energy fluctuation (standard deviation of the energy distribution) of a system in the $NVT$ ensemble is $\sigma_E = \sqrt{k_\mathrm{B} T^2 C_V}$, where $C_V$ is the heat capacity at constant volume. The heat capacity is typically proportional to the number of particles $N$, so the width of the energy distribution grows like $\sqrt{N}$. However, the mean energy typically grows like $N$, so if we double the size of the system, the energy distributions between simulations at 300 and 350 K will overlap less.

Hamiltonian replica exchange can be more targeted and is thus better suited for systems with large numbers of degrees of freedom. In this case, we have a base Hamiltonian, which is our unmodified and presumably realistic force field. The other replicas have an increasingly more biased potential energy. The biasing potential is chosen to help us to sample the transitions we want (for example, conformational transitions in a peptide), while not changing interactions for other degrees of freedom (for example, water–water interactions). Again, we expect faster exploration of phase space for the most biased replicas and the resulting conformations can be exchanged down to the base (unbiased) replica. We attempt exchanges every so many steps and accept them with the probability

$$ p_\mathrm{exchange} = \min \left\{ 1,\frac{\exp(-\beta V_i(\mathbf{X}_j) - \beta V_j(\mathbf{X}_i)}{\exp(-\beta V_i(\mathbf{X}_i) - \beta V_j(\mathbf{X}_j)} \right\},\tag{2}$$

where $\beta$ is the constant temperature for all replicas, $V_i$ is the potential energy function for replica $i$, and $\mathbf{X}_i$ is the set of all atomic coordinates for replica $i$ before the exchange attempt.