# Methods to handle non-ergodic systems with Molecular Dynamics

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?

• Can you give an example of what you mean? Because I wouldn't say that the MD simulation taking too long in practice is a problem of non-ergodicity. Like if some protein takes one second to fold, you probably can't use MD but this is not because the system is non-ergodic. Usually to circumvent this type of issue, people use something like replica-exchange MD, some type of MC, umbrella sampling, or metadynamics. All of these are methods that destroy the meaning of the sequence of frames in exchange for enhanced sampling of the system. – jheindel Oct 26 '20 at 17:41
• I am thinking of situations where there are energetic barriers that prevent the exploration of the entire relevant phase space on the timescale of the MD simulations, though where the physical system does explore these regions on experimental timescales (as an example; other similar situations are interesting too). While using the phrasing "non-ergodic" to describe such a situation may be using a slightly loose definition of ergodicty, I believe I have heard it used by some people in the field. – 2ndQuantized Oct 26 '20 at 21:21
• @jheindel You could write up an answer including the different techniques you mentioned (replica-exchange, umbrella sampling, etc.) as those are the kinds of techniques I am thinking of. If you or others know of more approaches, they could be included in answers too. – 2ndQuantized Oct 29 '20 at 21:35
• @jheindel it would be great if you could write up an answer! – Nike Dattani Nov 20 '20 at 22:46
• Problems with the ForceField accuracy are more likely to occur than ergodicity, particularly in liquid simulations. Ergodicity is more of a problem in glasses. Ultimately I think any method to overcome a true ergodicity issue, will incorporate Monte Carlo in some fashion. – Wesley Jul 23 at 16:43

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)

## 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.

• +10 for a very thorough addition of content here! Please remember to label your equations in case someone else has to refer to them, for example when they use our "Cite" button to generate a BibTeX formatted citation which is picked up by Google Scholar and other citation aggregators. – Nike Dattani Jul 22 at 21:53
• @NikeDattani I had no idea these equations could be cited and picked up by Google Scholar. Thanks for the info. – WaterMolecule Jul 22 at 23:40
• WaterMolecule: I don't know of any examples at MMSE yet, but you can see that MathOverflow posts have been cited at least 204 times in academic papers. 33 of those were since May 2021, so it's been becoming more and more popular to site SE posts in the academic literature. – Nike Dattani Jul 23 at 1:27
• Should that be $T_i$ in the definition of $\beta_i$? – 2ndQuantized Jul 23 at 17:17
• @2ndQuantized Yes, thanks for the correction. – WaterMolecule Jul 23 at 17:32