Chapter 0. Introduction

0. Learning objectives

Upon completion of this chapter, the student will understand

The motivation of enhanced sampling methods.
Key terminologies such as phase/configuration space, ensemble, Boltzmann distribution, collective variable, free energy profile, and potential of mean force.
The derivation of the Boltzmann distribution, potential of mean force, and free energy profile.

1. Motivation of enhanced sampling methods

As supercomputing resources have become increasingly available in recent years¹, molecular dynamics (MD) has become a standard tool for studying atomistic-level behaviors of molecular systems in various scientific disciplines². However, an MD simulation can accurately describe the thermodynamic and structural properties of a system only if it achieves ergodicity during its simulation timescale. This means that the simulation must explore the phase space comprehensively, or equivalently, ensure that all possible configurations are sampled. This is usually a challenging task in real-world systems, as the metastable states of interest are frequently separated by substantial free energy barriers that are difficult to overcome by random fluctuations within a reasonable amount of time (see Figure 1 below, or the leftmost panel of Figure 1 in Chapter 2).

Figure 1. A demonstration of sampling a 2D free energy surface $f (x_{1}, x_{2}) = x_{1}^{4} - 6 x_{1}^{2} + x_{1} + x_{2}^{2}$ , with the dots representing the position of the system in the 2D space $(x_{1}, x_{2})$ at different simulation steps. The figure shows that the system was unable to cross the free energy barrier to reach the other metastable state if the standard MD approach is used.

To overcome this limitation, researchers have developed various enhanced sampling methods in recent decades. One subset of these methods involves sampling along a set of predefined collective variables (CVs). These CVs, as we denote as a $d$ -dimensional vector $ξ = (ξ_{1} (x), ξ_{2} (x), \dots, ξ_{d} (x))$ , can be any functions of the system coordinates ( $x$ ), such as torsions, coordination numbers, or radii of gyration. Ideally, CVs should represent a small number of degrees of freedom that capture the slowest degrees of freedom of the system and distinguish metastable states of interest. Some examples of CV-based methods include metadynamics, adaptive biasing force, umbrella sampling, and replica-exchange umbrella sampling, which will all be covered in this mini-course.

2. Terminologies

Before diving into the following chapters, we will review some key concepts in statistical mechanics by introducing a few terminologies one may come across in enhanced sampling methods. This section serves as a quick review of the essential concepts to better understand the tutorials in later chapters. Readers are encouraged to visit Section 2.2 (Glossary of essential notions) in the review by Hénin et al.³ for a more exhaustive list of terms.

Timescale

The real timescale of a molecular system can vary significantly depending on the size and complexity of the system, as well as the process of interest. As illustrated in Figure 2, the timescales of molecular motions can range from femtoseconds to seconds or even longer. For instance, bond vibrations, the fastest motion in a molecular system, typically occur on the order of femtoseconds. In contrast, biological processes such as allosteric transitions or protein-ligand binding have timescales considerably longer than bond vibrations, making the transitions between metastable states of interest rare events. This long timescale is due to the high free energy barriers separating metastable states, as discussed in the previous section.

Importantly, to ensure that the motion of atoms in a system is captured and no short-lived events are missed, such as bond vibrations and molecular rotations, the time step used in molecular dynamics simulations must be smaller than the shortest timescale in the system. This inevitably small time step makes simulations that study biological processes computationally expensive, easily requiring $10^{8}$ steps or more. To reduce the computational cost of simulating long-timescale processes, researchers frequently use enhanced sampling methods to bias the system towards the metastable states of interest. These methods make rare events more probable in a shorter timescale, thus reducing the number of required simulation steps.

Figure 2. The timescales of different molecular motions.

Phase space and configuration space

In statistical mechanics, the phase space of a system is a space where all the possible states of the system can be represented. In three-dimensional Cartesian coordinate systems, a state is characterized by the coordinates $x \in R^{3 N}$ and momenta $p \in R^{3 N}$ of the system, where $N$ is the number of particles of the system. Consequently, each state can be represented by a unique point in the phase space, which is the joint space $(x, p)$ .

In contrast, the configuration space is defined by all the possible configurations (defined solely by $x$ ) of the system, so its dimensionality ( $3 N$ ) is half of that of the phase space ( $6 N$ ). When simulating long-timescale events using molecular dynamics, it is crucial to have sufficient sampling in the configuration space to capture the distribution of all possible configurations of the system. This ensures that rare events and transitions between different configurations are adequately represented. It is this need for comprehensive sampling in the configuration space motivates the development of enhanced sampling methods, which are the focus of this mini-course.

Microstate and macrostate

In the case where a system is composed of a larger number of particles, such as atoms or molecules, a microstate is a particular spatial arrangement of the particles in the system, i.e., a configuration. On the other hand, a macrostate, which can be defined by macroscopic thermodynamics variables (e.g., temperature, pressure, or total energy), is a collection of many possible microstates. For example, in a system composed of particles A and B that can either have an energy of 0 or 1, there are only three possible macrostates, which correspond to total energies of 0, 1, and 2, respectively. In the macrostate having a total energy of 2, there are two possible microstates: particles A and B having energies of 0 and 1, and the other way around.

Boltzmann distribution

To accurately calculate thermodynamic properties from molecular dynamics, we need to sample the Boltzmann distribution so that the system is at equilibrium. The Boltzmann distribution, sometimes called a Gibbs distribution, is a probability distribution that characterizes the distribution of states of a system in thermal equilibrium at a given temperature. It describes the probability of finding a particle in a particular energy state given its energy and the temperature of the system. Specifically, in the case where the energy levels are discrete, the normalized probability $p_{i}$ of finding the system being in state $i$ can be expressed as follows: $\begin{matrix} (1) & p_{i} = \frac{e^{- β E_{i}}}{\sum_{i} e^{- β E_{i}}} \propto e^{- β E_{i}} \end{matrix}$ where $E_{i}$ is the total energy (the sum of potential and kinetic energy) of state $i$ , and $β = 1 / k_{B} T$ is the inverse temperature (where $k_{B}$ is the Boltzmann constant ( $1.38 \times 10^{- 23} J \cdot K^{- 1}$ ), and $T$ is the system temperature). Frequently, we denote the normalization constant as $Q = \sum_{i} e^{- β E_{i}}$ , which is the so-called partition function.

In a phase space, where a state is defined by its continuous variables, coordinates $x$ and momenta $p$ , we can rewrite the Boltzmann distribution as $\begin{matrix} (2) & μ (x, p) = \frac{e^{- β (U (x) + K (p))}}{\int \int e^{- β (U (x) + K (p))} d x d p} \propto e^{- β (U (x) + K (p))} \end{matrix}$ where $μ (x, p)$ , analogous to $p_{i}$ in the discrete space, is the probability of finding the system to have coordinates $x$ and momenta $p$ . In this case, we have the partition function of the phase space $Q = \int \int e^{- β (U (x) + K (p))} d x d p = (\int e^{- β U (x)} d x) (\int e^{- β K (p)} d p)$ . Note that in Equation $2$ , we express the total energy as the sum of the potential energy $U (x)$ and the kinetic energy $K (p)$ .

From Equation $2$ , we can derive the configurational distribution (the Boltzmann distribution in the configuration space) by integrating $μ (x, p)$ with respect to the momenta $p$ : $\begin{matrix} (3) & ν (x) = \int μ (x, p) d p = \frac{1}{Q} \int e^{- β (U (x) + K (p))} d p = \frac{e^{- β U (x)} (\int e^{- β K (p)} d p)}{(\int e^{- β U (x)} d x) (\int e^{- β K (p)} d p)} \end{matrix}$ Finally, we have $\begin{matrix} (4) & ν (x) = \frac{e^{- β U (x)}}{\int e^{- β U (x)} d x} \propto e^{- β U (x)} \end{matrix}$ Frequently, Equation $4$ is written as $ν (x) = \frac{1}{Z} e^{- β U (x)}$ , where $Z = \int e^{- β U (x)} d x$ is the configurational partition function.

Supplementary note: Maxwell-Boltzmann distribution

Different from the Boltzmann distribution, the Maxwell-Boltzmann distribution is a special case of the Boltzmann distribution. It is a probability distribution that describes the speed ( $v$ ) of non-interacting gas particles in thermal equilibrium at a given temperature. In a three-dimensional velocity space, the Maxwell-Boltzmann distribution is given by $\begin{matrix} (S1) & f (v) = {(\frac{m}{2 π k_{B} T})}^{3 / 2} e^{- \frac{m v^{2}}{2 k_{B} T}} \propto e^{- \frac{m v^{2}}{2 k_{B} T}} \end{matrix}$

Supplementary note: Derivation of the Boltzmann distribution

The derivation of the Boltzmann distribution (Equation $2$ ) is based on the second law of thermodynamics, which states that any spontaneous process will tend to maximize the total entropy of the system. By definition, the entropy of a system can be expressed as $\begin{matrix} (S2) & S = k_{B} \ln W \end{matrix}$ where $W$ is the multiplicity, i.e., the number of microstates corresponding to a particular macrostate of a thermodynamic system. Now, we consider a system of $N$ particles, where $n_{i}$ particles have an energy of $E_{i}$ , i.e., $\sum_{i} n_{i} = N$ . Then, the multiplicity $W$ can be expressed as $W = \frac{N!}{n_{1}! \cdot n_{2}! \cdot n_{3}! \cdot \dots}$ . Hence, $\begin{matrix} (S3) & \ln W = \ln N! - \sum_{i} \ln n_{i}! \end{matrix}$ By Stirling’s approximation ( $\ln n! = n \ln n - n$ ), we have $\begin{matrix} (S4) & \ln W = (N \ln N - N) - (\sum_{i} n_{i} \ln n_{i} - \sum_{i} n_{i}) = N \ln N - \sum_{i} n_{i} \ln n_{i} \end{matrix}$ Let the probability of finding a particle having an energy $E_{i}$ be $p_{i} = n_{i} / N$ , we then have $\ln W = N \ln N - \sum_{i} (N p_{i}) \ln (N p_{i}) = N \ln N - N \sum_{i} p_{i} (\ln N + \ln p_{i})$ $\begin{matrix} (S5) & = N \ln N - N \ln N \sum_{i} p_{i} - N \sum_{i} p_{i} \ln p_{i} = - N \sum_{i} p_{i} \ln p_{i} \end{matrix}$ Given the entropy $S = k_{B} \ln W$ , we have $\begin{matrix} (S6) & S / k_{B} N = - \sum_{i} p_{i} \ln p_{i} \end{matrix}$ Now, the goal is to find $p_{i}$ that maximizes the entropy term $S / k_{B} N$ under the following two constraints:

The sum of the probabilities should be 1, i.e., $\sum_{i} p_{i} = 1$
The average energy of the particles should be $\bar{E}$ , i.e., $\sum_{i} p_{i} E_{i} = \bar{E}$

To solve for $p_{i}$ that maximize $f$ , we consider the following function: $\begin{matrix} (S7) & f (p_{i}) = - \sum_{i} p_{i} \ln p_{i} - α (\sum_{i} p_{i}) - β (\sum_{i} p_{i} E_{i}) \end{matrix}$ where $α$ , $β$ are Language multipliers. (Note that $E_{i}$ ’s are a bunch of constants, not variables as $p_{i}$ .) Then, we equate the partial derivatives of $f$ with respect to $p_{i}$ to 0, namely, $\frac{\partial f}{\partial p_{i}} = 0 \Rightarrow - \sum_{i} \ln p_{i} - \sum_{i} p_{i} \cdot \frac{1}{p_{i}} - α - β (\sum_{i} E_{i}) = 0$ $\begin{matrix} (S8) & \Rightarrow \sum_{i} (- \ln p_{i} - 1 - α - β E_{i}) = 0 \Rightarrow p_{i} = e^{- (α + 1)} e^{- β E_{i}} \propto e^{- β E_{i}} \end{matrix}$ Applying the normalization constraint $\sum_{i} p_{i} = 1$ , we have $\begin{matrix} (S9) & e^{- (α + 1)} \sum_{i} e^{- β E_{i}} = 1 \Rightarrow e^{- (α + 1)} = \frac{1}{e^{- β E_{i}}} \end{matrix}$ which means that $e^{- (α + 1)}$ is the normalization constant. To get the expression of $β$ , we consider the Maxwell relation ${(\frac{\partial S}{\partial U})}_{V, N} = \frac{1}{T}$ .

Specifically, we have $\begin{matrix} (S10) & d U = d (N \bar{E}) = N d (\sum_{i} p_{i} E_{i}) = N \sum_{i} E_{i} d p_{i} \end{matrix}$ $d S = - k_{B} N d (\sum_{i} p_{i} \ln p_{i}) = - k_{B} N \sum_{i} (\ln p_{i} + 1) d p_{i}$ $\begin{matrix} (S11) & \Rightarrow d S = - k_{B} N \sum_{i} (- α - β E_{i}) d p_{i} \end{matrix}$ Given that $\sum_{i} p_{i} = 1 \Rightarrow \sum_{i} d p_{i} = 0$ , $d S$ can be simplied as follows: $\begin{matrix} (S12) & d S = α k_{B} N {\sum_{i} d p_{i}}^{0} + β k_{B} (N \sum_{i} E_{i} d p_{i}) = β k_{B} d U \end{matrix}$ Hence, $\begin{matrix} (S13) & {(\frac{\partial S}{\partial U})}_{V, N} = β k_{B} = \frac{1}{T} \Rightarrow β = \frac{1}{k_{B} T} \end{matrix}$ Altogether, Equations $S8$ and $S13$ derive the expression of the Boltzmann distribution, i.e., Equation $1$ .

Ensemble

In statistical mechanics, an ensemble is a collection of hypothetical copies of a system that all share the same macroscopic properties (such as temperature, pressure, volume, and chemical potential) but may differ in their microscopic details (such as the positions and velocities of individual particles). Some of the most common thermodynamic ensembles include

NVT ensemble: The NVT ensemble is also called the canonical ensemble, where the number of particles ( $N$ ), volume ( $V$ ), and temperature ( $T$ ) are fixed. The probability distribution and the corresponding partition function for the NVT ensemble are given by the Boltzmann distribution, namely, $\begin{matrix} (5) & p_{i} = \frac{e^{- β E_{i}}}{\sum_{i} e^{- β E_{i}}}, Z_{N V T} = \sum_{i} e^{- β E_{i}} \end{matrix}$
NPT ensemble: The NPT ensemble is also called the isothermal-isobaric ensemble, or the Gibbs ensemble, where the number of particles ( $N$ ), pressure( $P$ ), and temperature ( $T$ ) are fixed. The probability distribution and partition function of the NPT ensemble is given below: $\begin{matrix} (6) & p_{i} = \frac{e^{- β E_{i}} e^{- β P V_{i}}}{\sum_{i} e^{- β E_{i}} e^{- β P V_{i}}}, Z_{N P T} = \sum_{i} e^{- β E_{i}} e^{- β P V_{i}} \end{matrix}$
NVE ensemble: The NVE ensemble is also called the microcanonical ensemble. The NVE ensemble is defined by assigning an equal probability to every microstate whose energy falls within a range centered at $E$ , and all the other microstates are given a probability of 0. That is, the probability in the NVE ensemble can be written as $p_{i} = 1 / W$ , where $W$ is the number of microstates within the range of energy.
$μ$ VT ensemble: The $μ$ VT ensemble is also called the grand canonical ensemble, where the chemical potential ( $μ$ ), volume ( $V$ ) and temperature ( $T$ ) are fixed. The probability distribution and partition function of the $μ V T$ ensemble is given below: $\begin{matrix} (7) & p_{i} = \frac{e^{- β E_{i}} e^{β μ N_{i}}}{\sum_{i} e^{- β E_{i}} e^{β μ N_{i}}}, Z_{μ V T} = \sum_{i} e^{- β E_{i}} e^{β μ N_{i}} \end{matrix}$

Notably, given the copies of the system in an ensemble, we can calculate the ensemble average of an observable $O (x, p)$ of interest, which is simply taking the average of $O (x, p)$ considering all points $(x, p)$ in the phase space: $\begin{matrix} (8) & ⟨ O ⟩ = \int O (x, p) μ (x, p) d x d p \end{matrix}$

Ergodicity

The dynamics of a system is said to be ergodic if the ensemble average of the observable $O (x, p)$ of interest is equal to the time average of the observable over an infinitely long trajectory, namely, $\begin{matrix} (9) & ⟨ O ⟩ = \int O (x, p) μ (x, p) d x d p = lim_{M \to \infty} \sum_{t = 1}^{M} O (x_{t}, p_{t}) \end{matrix}$ In simple terms, if one monitors an obervable of interest for a single particle over a sufficiently long period of time and take the average of the observable at all time frames, the outcome would be equivalent to the average of the same observable computed over all possible configurations of the system at a single instant. In practice, given that it is impractical to identify all possible configurations at any given time frame, reaching ergodicity in molecular dynamics is essential as it allows one to estimate the ensemble average by computing the time average over a sufficiently long trajectory.

Metastable state

Metastable states are the states having a higher potential energy than the ground state but are stable enough for the system to reside for an extended (but finite) amount of time. For example, there are two metastable states in Figure 1. One has values of $x_{2}$ around -1.8, while the other has values of $x_{2}$ around 1.8.

Collective variable (CV)

Here, we refer to collective variables (CVs) as a $d$ -dimensional vector $ξ = (ξ_{1} (x), ξ_{2} (x), \dots, ξ_{d} (x))$ , which can be seen as a lower-dimensional representation reduced/projected from the full $3 N$ -dimensional configurations $x$ via the mapping function $Φ (x) = ξ$ , with $d ≪ 3 N$ . Traditionally, CVs can be defined as any functions of the coordinates of the system, but could additionally include the alchemical variable, as shown in a recent paper proposing alchemical metadynamics⁴, the method we will explore in Chapter 7.

Ideally, an optimal set of CVs should be able to capture the slowest dynamics of the system and distinguish the metastable states of interest. For example, in Figure 1, $x_{1}$ is a reliable CV for distinguishing the two metastable states of interest. Conversely, $x_{2}$ is orthogonal to the slowest degree of freedom of the system and, therefore, is not suitable for describing transitions between the two metastable states. Identifying the optimal set of CVs is often challenging in systems where the slowest-relaxing coordinates are not intuitive. In such cases, machine learning-based algorithms like RAVE⁵ or SGOOP⁶ can be helpful.

Free energy

In the canonical (NVT) ensemble, Helmholtz free energy is defined as $\begin{matrix} (10) & F = - \frac{1}{β} \ln Z_{N V T} = - \frac{1}{β} \int e^{- β U (x)} d x \end{matrix}$ while in the Gibbs (NPT) ensemble, the Gibbs free energy can be expessed as $\begin{matrix} (11) & G = F + P V = - \frac{1}{β} [\ln Z_{N V T} - (\frac{\partial \ln Z_{N V T}}{\partial \ln V})] \end{matrix}$ Practically, in a soluted system where the solvent is nearly imcompressible, i.e., $\frac{\partial V}{\partial P} \approx 0$ , the difference in Gibbs free energy is almost equal to the difference in Helmholtz free energy.

Free energy profile

A free energy profile is meaningful only if there is a collective variable $ξ$ of interest. A free energy profile, sometimes called a free energy landscape or a free energy surface, can be computed as follows: $\begin{matrix} (12) & F (ξ) = - \frac{1}{β} \ln P (ξ) = - \frac{1}{β} \ln \int δ (ξ - Φ (x)) ν (x) d x = - \frac{1}{β} ⟨ δ (ξ - ξ (x)) ⟩ \end{matrix}$ where $ν (x) = e^{- β U (x)} / \int e^{- β U (x)} d x$ is the Boltzmann distribution and $Φ$ is a function mapping the $3 N$ -dimensional vector $x$ to the $d$ -dimensional CV vector $ξ$ . The function $δ$ is the multivariate Dirac delta function, i.e., $\begin{matrix} (13) & δ (ξ - Φ (x)) = \prod_{i}^{d} δ^{'} (ξ_{i} - Φ_{i} (x)), δ^{'} (ξ_{i} - Φ_{i} (x)) = {\begin{cases} + \infty, if ξ_{i} = Φ_{i} (x) \\ 0, otherwise \end{cases} \end{matrix}$ where $δ^{'}$ is a single-variable Dirac delta function and $Φ_{i}$ is a function mapping $x$ to the $i$ -th CV, $ξ_{i}$ . While Equations $12$ and $13$ may seem daunting, they essentially mean that the free energy profile along the CVs can be computed from the marginal probability density $P (x)$ along the CVs, which is obtained by integrating the Boltzmann distribution $ν (x)$ over all variables except for the CVs $Φ (x)$ . This interpretation may be understood more easily with the example in the Supplementary note below.

Supplementary note: An example of analytical free energy profile

In this example, given the potential energy $U (x_{1}, x_{2}) = x_{1}^{4} - 6 x_{1}^{2} + x_{1} + x_{2}^{2}$ (the same functional form as the one in the free energy surface in Figure 1), we want to determine its free energy profile along the collective variable $x_{1}$ . First, we can calculate the Boltzmann distribution $\begin{matrix} (S14) & ν (x) = \frac{1}{Z} \exp (- β (x_{1}^{4} - 6 x_{1}^{2} + x_{1})) \exp (- β x_{2}^{2}) \end{matrix}$ where $Z = \int \int \exp (- β (x_{1}^{4} - 6 x_{1}^{2} + x_{1})) \exp (- β x_{2}^{2}) d x_{1} d x_{2}$ is a spearable integral. To save space, we let $f (x_{1}) = x_{1}^{4} - 6 x_{1}^{2} + x_{1}$ , so $ν (x) = \frac{1}{Z} e^{- β f (x_{1})} e^{- β x_{2}^{2}}$ and $Z = (\int e^{- β f (x_{1})} d x_{1}) (\int e^{- β x_{2}^{2}} d x_{2})$ . Then, given that $x_{1}$ is the only CV, we can compute the free energy profile $F (x_{1})$ as follows: $F (x_{1}) = - \frac{1}{β} \ln (\int \int δ^{'} (x_{1} - Φ_{1} (x)) ν (x_{1}, x_{2}) d x_{1} d x_{2})$ $\begin{matrix} (S15) & = - \frac{1}{β} \ln [(\int_{- \infty}^{+ \infty} δ^{'} (x_{1} - Φ_{1} (x)) e^{- β f (x_{1})} d x_{1}) (\int_{- \infty}^{+ \infty} e^{- β x_{2}^{2}} d x_{2})] \end{matrix}$ where $Φ_{1} (x)$ takes in $x$ (and maps it to a specific value of $x_{1}$ ), but is a function of $x_{1}$ and therefore, can be grouped into the first integral in Eequation $S15$ . Importantly, the value $x_{1}$ runs from $- \infty$ to $+ \infty$ and whenever it is not equal to $Φ_{1} (x)$ , the detal function $δ^{'} (x_{1} - Φ_{1} (x))$ returns 0. In fact, given $\int f (x) δ^{'} (x - y) d x = f (y)$ (and $(\int_{- \infty}^{+ \infty} e^{- β x_{2}^{2}} d x_{2}) = \sqrt{π / β}$ ), we have $\begin{matrix} (S16) & F (x_{1}) = - \frac{1}{β} \ln (e^{- β f (x_{1})} \cdot \sqrt{\frac{π}{β}}) \Rightarrow F (x_{1}) = x_{1}^{4} - 6 x_{1}^{2} + x_{1} - \frac{1}{β} \ln (\sqrt{\frac{π}{β}}) \end{matrix}$ The same result can be reached by integrating the Boltzmann distribution overall all variables but not $x_{1}$ , i.e., $\begin{matrix} (S17) & F (x_{1}) = - \frac{1}{β} \ln (\int e^{- β U (x_{1}, x_{2})} d x_{2}) = x_{1}^{4} - 6 x_{1}^{2} + x_{1} - \frac{1}{β} \ln (\sqrt{\frac{π}{β}}) \end{matrix}$

Free energy estimator

A free energy estimator is an expression or algorithm that is used to estimate free energy differences or free energy surfaces from simulation data, such as configurations, energies, and forces. Other than estimating free energy by directly measuring the probability ratios or from transition matrices, common free energy estimators include weighted histogram analysis method (WHAM)⁷, thermodynamic integration (TI)⁸, Bennett acceptance ratio (BAR)⁹, and multistate Bennett acceptance ratio (MBAR)¹⁰. These free energy estimators will be briefly introduced in the following chapters.

Potential of mean force (PMF)

While the term potential of mean force (PMF) is frequently used interchangeably with free energy surface, the historical definition of PMF is defined between two particles 1 and 2 and represents the effective potential energy that particle 1 experiences due to the presence of particle 2, when the positions of all other particles in the system are averaged over. It is defined based on the radial distribution function $g (r)$ : $\begin{matrix} (14) & W (r) = - \frac{1}{β} \ln g (r) = F (r) + \frac{2}{β} \ln r + C \end{matrix}$ where $F (r)$ is the free energy surface along the interparticle distance and $C$ is an arbitrary constant. That is, the historic definition of PMF is related to the free energy surface but is still a distinct quantity. Check the Supplementary note below to see how PMF relates to the radial distribution function. In this mini-course, we treat the term “potential of mean force” equivalently as the term “free energy profile”.

Supplementary note: Relationship between the potential of mean force and distribution/correlation functions

1. General distribution functions
Here we consider a canonical (NVT) ensemble characterized by

N

particles in a volume

V

(hence the average number density

ρ = N / V

) at temperature

T

. (Note that the theory also works for other thermodynamic ensembles.) Given

U_{N} = U_{N} (r_{1}, . . ., r_{N})

, the potential energy due to the interaction between particles, the configurational partition function is given by

\begin{matrix} (S18) & Z_{N} = \int . . . \int e^{- β U_{N} (r_{1}, . . ., r_{N})} d r_{1} . . . d r_{N} \end{matrix}

where

r_{1}, . . ., r_{N}

are the position vectors of the particles. With this, the probability of finding particle 1 in the volume element

d r_{1}

at point

r_{1}

, etc., can be expressed as

\begin{matrix} (S19) & P^{(N)} (r_{1}, . . ., r_{N}) d r_{1} . . . d r_{N} = \frac{e^{- β U_{N}}}{Z_{N}} d r_{1} . . . d r_{N} \end{matrix}

To determine the probability for the case where the positions of

n < N

particles are fixed in the volume elements

d r_{1}, . . ., d r_{n}

at the point

r_{1}, . . ., r_{n}

, with the remaining

N - n

particles unconstrained, we simply integrate Equation

S19

over the particles with indices

n + 1, . . ., N

, namely,

$\begin{matrix} (S20) & P^{(n)} (r_{1}, . . ., r_{n}) = \int . . . \int \frac{e^{- β U_{N}}}{Z_{N}} d r_{n + 1} . . . d r_{N} \end{matrix}$

Assuming that the particles are indistinguishable, here we define the n-particle density function:

ρ^{(n)} (d r_{1} . . . d r_{N}) = \frac{N!}{(N - n)!} P^{(N)} (r_{1}, . . ., r_{n})

\begin{matrix} (S21) & \Rightarrow ρ^{(n)} (d r_{1} . . . d r_{N}) = \frac{N!}{(N - n)!} \int . . . \int \frac{e^{- β U_{N}}}{Z_{N}} d r_{n + 1} . . . d r_{N} \end{matrix}

which represents the probability density of finding any particle at

r_{1}

, any particle at

r_{2}

, ..., and any particle at

r_{n}

. (Note that for

n = 1

(one-particle density), we have

ρ^{(1)} (r) = ρ

, the bulk number density. See Wikipedia or LibreTexts for the proof.)

2. General correlation functions
A general correlation function can be defined in terms of the distribution function $ρ^{(n)} (r_{1}, . . ., r_{n})$ as $g^{(n)} (r_{1}, . . ., r_{n}) = ρ^{(n)} (r_{1}, . . ., r_{n}) / ρ^{n}$ $\begin{matrix} (S22) & \Rightarrow g^{(n)} (r_{1}, . . ., r_{n}) = {(\frac{V}{N})}^{n} \frac{N!}{(N - n)!} \int . . . \int \frac{e^{- β U_{N}}}{Z_{N}} d r_{n + 1} . . . d r_{N} \end{matrix}$ Specifically, for $n = 2$ , we have the so-called pair correlation functions. Note that $g^{(2)} (r_{1}, r_{2})$ can be rewritten as $g (r) = g (| r_{1} - r_{2} |)$ , so that a pair of particles (with one of them being at the origin) separated by a distance $r$ is considered. $\begin{matrix} (S23) & g^{(2)} (r_{1}, r_{2}) = g (r) = \frac{1}{ρ^{2}} N (N - 1) \int . . . \int \frac{e^{- β U_{N}}}{Z_{N}} d r_{3} . . . d r_{N} \end{matrix}$ 3. Potential of mean force (PMF)
Finally, here we relate the pair correlation function/radial distribution function $g (r)$ to the mean force $⟨ F ⟩$ , which is the negative derivative of the potential of mean force (PMF) $W$ . Here, we consider the mean force that particle 1 experiences due to the presence of particle 2: $⟨ F_{1} ⟩ = - ⟨ \frac{d}{d r_{1}} U_{N} (r_{1}, . . ., r_{N}) ⟩_{r_{1}, r_{2}} = \frac{- \int . . . \int (\frac{d U_{N}}{d r_{1}}) e^{- β U_{N}} d r_{3} . . . d r_{N}}{\int . . . \int e^{- β U} d r_{3} . . . d r_{N}}$ $\begin{matrix} (S24) & = \frac{1}{β} \frac{\frac{d}{d r_{1}} (\int . . . \int e^{- β U_{N}} d r_{3} . . . d r_{N})}{\int . . . \int e^{- β U_{N}} d r_{3} . . . d r_{N}} = \frac{1}{β} \frac{d}{d r_{1}} \ln (\int . . . \int e^{- β U_{N}} d r_{3} . . . d r_{N}) \end{matrix}$ From Equation $S23$ , we have $\begin{matrix} (S25) & \int . . . \int e^{- β U_{N}} d r_{3} . . . d r_{N} = \frac{ρ^{2} Z_{N}}{N (N - 1)} g^{(2)} (r_{1}, r_{2}) \end{matrix}$ Hence $\begin{matrix} (S26) & ⟨ F_{1} ⟩ = \frac{1}{β} \frac{d}{d r_{1}} [\ln (\frac{ρ^{2} Z_{N}}{N (N - 1)}) + \ln (g^{(2)} (r_{1}, r_{2}))] = \frac{1}{β} \frac{d}{r_{1}} \ln g (r) \end{matrix}$ Finally, we have $\begin{matrix} (S27) & W (r) = - \int ⟨ F_{1} ⟩ d r_{1} = \int (\frac{1}{β} \frac{d}{d r_{1}} \ln g (r)) d r_{1} \Rightarrow W (r) = - \frac{1}{β} \ln g (r) \end{matrix}$ (Check this handout to see how the radial distribution function $g (r)$ can be used to calculate other thermodynamic quantities such as the energy or pressure.)

3. Exercises

👨‍💻 Exercise 1

Given the functional form of the free energy surface in Figure 1 ( $f (x_{1}, x_{2}) = x_{1}^{4} - 6 x_{1}^{2} + x_{1} + x_{2}^{2}$ ), estimate the timescale of transitioning from the left to right metastable states at 300 K. (Hint: One can use the transition state theory (TST) to estimate the rate constant of the process.)

Click here for the solution.

First, we need to figure out the height of the free energy barrier going from left to right metastable states. This requires determining the local minima and maxima of the free energy surface by equating the partial derivative of $f (x_{1}, x_{2})$ with respect to $x_{1}$ : $\frac{\partial f}{\partial x_{1}} = 0 \Rightarrow 4 x_{1}^{3} - 12 x_{1} + 1 = 0$ As a result, at $x_{1} \approx - 1.7723$ , $f (x_{1})$ has the global minimum of $- 10.7524 k_{B} T$ , while a local maximum of $0.0417 k_{B} T$ occurs at $x_{1} \approx 0.0835$ , so the barrier height is $10.7941 k_{B} T$ . Given the free energy barrier height ( $Δ G^{‡} = 10.7941 k_{B} T$ ), now we can use transition state theory (TST) to determine the rate constant as follows: $\begin{matrix} (E1) & k = \frac{k_{B} T}{h} \exp (\frac{- Δ G^{‡}}{k_{B} T}) \end{matrix}$ Accordingly, we have $\begin{matrix} (E2) & k = \frac{1.38 \times 10^{- 23} J / K \times 300 K}{6.626 \times 10^{- 34} J \cdot s} \exp (- 10.7524) \approx 1.34 \times 10^{8} s^{- 1} \end{matrix}$ Therefore, the timescale of the process is $t = 1 / k \approx 7.5 \times 10^{- 9} s = 7.5 ns$ . That is, with a time step of 2 femtoseconds, it would require at least 3.7 million steps to simulate this process in a standard molecular dynamics simulation. Given currently available high-end computing resources, such a computational cost is not too demanding. However, for an even larger free energy barrier, say, $20 k_{B} T$ , the timescale quickly goes up to 77.6 microseconds. For a barrier of $30 k_{B} T$ , which is not uncommon,

👨‍💻 Exercise 2

A one-dimensional simple harmonic potential with a force constant $K$ can be expressed as $U (x) = \frac{1}{2} K (x - x_{0})^{2}$ , where $x_{0}$ is the equilibrium position on the $x$ -axis. Given this, derive the probability density distribution and Helmholtz free energy expressions for the system.

Click here for the solution.

The one-dimensional harmonic potential is given by $\begin{matrix} (E3) & U (x) = \frac{1}{2} K (x - x_{0})^{2} \end{matrix}$ Substituting Equation $E3$ into Equation $4$ , we have $\begin{matrix} (E4) & ν (x) = \frac{\exp (- \frac{1}{2} β K (x - x_{0})^{2})}{\int_{- \infty}^{\infty} \exp (- \frac{1}{2} β K (x - x_{0})^{2}) d x} \end{matrix}$ Let $a = \sqrt{\frac{β K}{2}} (x - x_{0}) \Rightarrow d a = \sqrt{\frac{β K}{2}} d x$ , then the denominator $Z$ of Equation $E4$ can be calculated as below: $\begin{matrix} (E5) & Z = \int_{- \infty}^{\infty} \exp (- a^{2}) d x = \sqrt{\frac{2}{β K}} \int_{- \infty}^{\infty} \exp (- a^{2}) d a = \sqrt{\frac{2 π}{β K}} \end{matrix}$ Therefore, the probability distribution of a one-dimensional harmonic oscillator is $\begin{matrix} (E6) & ν (x) = \sqrt{\frac{β K}{2 π}} \exp (- \frac{β K}{2} (x - x_{0})^{2}) \end{matrix}$ Importantly, this is a Gaussian distribution centered at $x_{0}$ , with the standard deviation being $σ = 1 / \sqrt{β K}$ . Then, using $F = - \frac{1}{β} \ln Z$ , we can compute the free energy as $\begin{matrix} (E7) & F = - \frac{1}{β} \ln (\sqrt{\frac{2 π}{β K}}) = - \frac{1}{2 β} \ln (2 π σ) \end{matrix}$

Padua, D. Encyclopedia of Parallel Computing; Springer Science & Business Media, 2011. ↩︎
Gelpi, J.; Hospital, A.; Goñi, R.; Orozco, M. Molecular Dynamics Simulations: Advances and Applications. Adv. Appl. Bioinform. Chem. 2015, 8, 37–47. ↩︎
HéninJ.; LelièvreT.; Shirts, M. R.; Valsson, O.; Delemotte, L. Enhanced Sampling Methods for Molecular Dynamics Simulations. Living J. Comput. Mol. Sci. 2022, 4 (1). ↩︎
Hsu, W.-T.; Piomponi, V.; Merz, P. T.; Bussi, G.; Shirts, M. R. Alchemical Metadynamics: Adding Alchemical Variables to Metadynamics to Enhance Sampling in Free Energy Calculations. Journal of Chemical Theory and Computation 2023. ↩︎
Ribeiro, J. M. L.; Bravo, P.; Wang, Y.; Tiwary, P. Reweighted Autoencoded Variational Bayes for Enhanced Sampling (RAVE). The Journal of Chemical Physics 2018, 149 (7), 072301. ↩︎
Tiwary, P.; Berne, B. J. Spectral Gap Optimization of Order Parameters for Sampling Complex Molecular Systems. Proceedings of the National Academy of Sciences 2016, 113 (11), 2839–2844. ↩︎
Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. The Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992, 13 (8), 1011–1021. ↩︎
Kirkwood, J. G. Statistical Mechanics of Fluid Mixtures. J. Chem. Phys. 1935, 3 (5), 300–313. ↩︎
Bennett, C. H. Efficient Estimation of Free Energy Differences from Monte Carlo Data. J. Comput. Phys. 1976, 22 (2), 245–268. ↩︎
Shirts, M. R.; Chodera, J. D. Statistically Optimal Analysis of Samples from Multiple Equilibrium States. J. Chem. Phys. 2008, 129 (12), 124105. ↩︎

Last updated on Sep 20, 2024