OPES (on-the-fly probability enhanced sampling)

💭Catalogs

• Fundamental knowledge
• Reconstruct the unbiased probability
• Update bias based on probability
• On-the-fly probability?
• OPES-explore
• Differences between OPES and OPES-explore

Fundamental knowledge

Considering a system with an interaction potential $U(\mathbf{R})$, where $\mathbf{R}$ denotes the atomic coordinates. Sampling is accelerated by adding a bias potential $V(\mathbf{s})$ that depends on $\mathbf{R}$ via a set of collective variables (CVs), $\mathbf{s}=\mathbf{s}(\mathbf{R})$. The CVs are chosen so as to describe the modes of the system that are more difficult to sample. The choice of a proper set of CVs is critical, as it determines the efficiency of the method. The properties of the unbiased system are then calculated by using a reweighting procedure.

In fact the unbiased probability density $$\begin{array}{c}P(\mathbf{s})=⟨\delta [\mathbf{s}-\mathbf{s}(\mathbf{R})]⟩\propto \int d\mathbf{R}{e}^{-\beta U(\mathbf{R})}\delta [\mathbf{s}-\mathbf{s}(\mathbf{R})]\end{array}$$ can be written as an average over the biased ensemble: $$\begin{array}{c}P(\mathbf{s})=\frac{⟨\delta [\mathbf{s}-\mathbf{s}(\mathbf{R})]e^{\beta V(\mathbf{s})}{⟩}_V}{⟨e^{\beta V(\mathbf{s})}{⟩}_V}\end{array}$$ where $\beta$ is the inverse temperature.

In this way it is also possible to reconstruct the free energy surface (FES), defined as $$\begin{array}{c}F(\mathbf{s})=-\frac{1}{\beta }\log P(\mathbf{s})\end{array}$$

Bias $V(\mathbf{s})$ can be builded by adding at fixed intervals repulsive Gaussians centered at the instantaneous point sampled. At the $n$-th iteration the bias is given by: $$\begin{array}{c}{V}_n(\mathbf{s})=\sum _k^n{e}^{-\beta V_{k-1}({\mathbf{s}}_k)/(\gamma -1)}G(\mathbf{s},{\mathbf{s}}_k)\end{array}$$ where the parameter $\gamma >1$ is called the bias factor, and the Gaussian function is defined as $$\begin{array}{c}G(\mathbf{s},{\mathbf{s}}^{\prime })=h\exp \left[-\frac{1}{2}(\mathbf{s}-{\mathbf{s}}^{\prime }{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{s}-{\mathbf{s}}^{\prime })\right]\end{array}$$ with height $h$ and variance $\mathbf{\Sigma }$ set by the user. Typically only diagonal variances ${\Sigma }_{ij}={\sigma }_i^2{\delta }_{ij}$ are employed. At convergence there is a simple relation between the bias and the free energy $$\begin{array}{c}V(\mathbf{s})=-(1-1/\gamma )F(\mathbf{s})\end{array}$$ and the sampled ensemble is a smoothed version of the unbiased one, with FES barriers lowered by a factor $\gamma$.

Reconstruct the unbiased probability

Enhanced sampling based on the probability reconstruction is not a new idea. It was first proposed by the adaptive umbrella sampling method. Typically, in such methods the bias at $n$-th iteration is defined as: $$\begin{array}{c}{V}_n(\mathbf{s})=\frac{1}{\beta }\log {\hat{P}}_n(\mathbf{s}),\end{array}$$ where ${\hat{P}}_n(\mathbf{s})$ is an estimate of the probability obtained via a weighted histogram or some more elaborate method.

In building our method we will introduce few key differences that come from the long experience with MetaD.

First we would like to introduce explicitly a target distribution $p^{tg}(\mathbf{s})$, that will be sampled once the method reaches convergence. This can be obtained with the following bias: $$\begin{array}{c}V(\mathbf{s})=-\frac{1}{\beta }\log \frac{p^{tg}(\mathbf{s})}{P(\mathbf{s})}\end{array}$$ In adaptive umbrella sampling the target distribution is uniform $$\begin{array}{c}{p}^{tg}(\mathbf{s})\propto 1\end{array}$$ while in MetaD it is the well-tempered distribution $$\begin{array}{c}{p}^{tg}(\mathbf{s})=p^{\text{WT}}(\mathbf{s})\propto [P(\mathbf{s}){]}^{1/\gamma }\end{array}$$

Here we express the target distribution as a function of the unbiased one, $p^{tg}(\mathbf{s})\propto [P(\mathbf{s}){]}^{1/\gamma }$, we only need to estimate $P(\mathbf{s})$ via reweighting in order to calculate the bias.

We build our probability distribution estimate on the fly by periodically depositing Gaussians, known as kernel density estimation (KDE). Each new Gaussian is weighted according to the previously deposited bias potential: $$\begin{array}{c}{\tilde{P}}_n(\mathbf{s})=\frac{{\sum }_k^n{w}_{k}G(\mathbf{s},{\mathbf{s}}_k)}{{\sum }_k^n{w}_k}\end{array}$$ where the weights $w_k$ are given by $$\begin{array}{c}{w}_k=e^{\beta V_{k-1}({\mathbf{s}}_k)}\end{array}$$

${\tilde{P}}_n(\mathbf{s})$ is not properly normalized, and we will take care of the normalization separately.

The $G(\mathbf{s},{\mathbf{s}}_k)$ are Gaussians as those defined previously for MetaD, with diagonal variance $$\begin{array}{c}{\Sigma }_{ij}={\sigma }_i^2{\delta }_{ij}\end{array}$$ and fixed height $$\begin{array}{c}h=\prod _i{\left({\sigma }_i\sqrt{2\pi }\right)}^{-1}\end{array}$$ In KDE the most relevant parameter is the bandwidth, i.e., the width of the Gaussians.

A good choice of the bandwidth should depend on the amount of available data: the larger the sampling the smaller the bandwidth. Thus we choose to shrink the bandwidth as the simulation proceeds according to the popular Silverman's rule of thumb.

At $n$-th iteration: $$\begin{array}{c}{\sigma }_i^{(n)}={\sigma }_i^{(0)}[N_{\text{eff}}^{(n)}(d+2)/4{]}^{-1/(d+4)}\end{array}$$ where ${\sigma }_i^{(0)}$ is the initial standard deviation estimated from a short unbiased simulation, $d$ is the dimensionality of the CV space, and $$\begin{array}{c}{N}_{\text{eff}}^{(n)}=\frac{{\left({\sum }_k^n{w}_k\right)}^2}{{\sum }_k^n{w}_k^2}\end{array}$$ is the effective sample size.

Update bias based on probability

We can now discuss the normalization problem.

${\tilde{P}}_n(\mathbf{s})$ should be normalized not with respect to the full CV space, but only over the CV space actually explored up to step $n$, that we call $|{\Omega }_n|$. Thus we introduce the normalization factor $$\begin{array}{c}{Z}_n=\frac{1}{|{\Omega }_n|}{\int }_{{\Omega }_n}{\tilde{P}}_n(\mathbf{s})d\mathbf{s},\end{array}$$ that will change over time, as the system explores new regions of the CV space, and it will have an impact in the biasing scheme. This impact becomes particularly relevant in CV spaces of dimension $d\gg 1$, since the volume explored ${\Omega }_n$ grows with a power of $d$.

To estimate $Z_n$ we take advantage of our compressed kernels representation, and consider the centers of the kernels as points for a Monte Carlo integration $$\begin{array}{c}{Z}_n=\frac{1}{N}\sum _k^N\tilde{P}({\mathbf{s}}_{\mathbf{k}})=\frac{1}{NS}\sum _{k,k^{\prime }}^{N}G({\mathbf{s}}_{\mathbf{k}},{\mathbf{s}}_{{\mathbf{k}}^{\prime }})\end{array}$$ where $G({\mathbf{s}}_{\mathbf{k}},{\mathbf{s}}_{{\mathbf{k}}^{\prime }})$ are the compressed Gaussians, $N$ is their total number, and $$\begin{array}{c}S=\sum _k^n{w}_k\end{array}$$ is the global normalization of the KDE.

Finally we can explicitly write the bias at the $n$-th step as: $$\begin{array}{c}{V}_n(\mathbf{s})=(1-1/\gamma )\frac{1}{\beta }\log \left(\frac{{\tilde{P}}_n(\mathbf{s})}{Z_n}+ϵ\right)\end{array}$$ where $ϵ\ll 1$ can be seen as a regularization term that ensures the argument of the logarithm is always greater than zero. We notice that the addition of this term is not merely a technicality to solve a numerical issue, but rather it allows one to set a limit to the bias, thus providing a better control over the desired exploration. It can be chosen to be $$\begin{array}{c}ϵ=e^{-\beta \Delta E/(1-1/\gamma )}\end{array}$$ where $\Delta E$ is the height of the free energy barriers one wishes to overcome.

On-the-fly probability?

If you are having difficulties following the above instructions, allow me to summarize six steps to aid in understanding the essential aspects of the OPES method.

1. The only input parameter considered is $\Delta E$ for the rare event.

2. The initial bias $V$, i.e., minimum bias, is obtained as $-\Delta E$.

3. The Gaussian kernel is constructed using this above bias.

4. The unbiased probability $P$ is then constructed based on a weighted KDE approach (Reweighting!).

5. Update the bias based on the above $P$ and guidelines (Biasing!).

Repeat steps 2 to 5 until the regime is quasi-static, meaning that the the unbiased probability $P$ in a on-the-fly way.

6. Finally, the free energy profile can be calculated through unbiased probability $P$.

OPES-explore

OPES method focuses on fast convergence (e.g. using DPMD to calculate the free energy profiles), but there are cases where fast exploration is preferred instead (e.g. using DP-GEN iterations to build training datesets).

For this reason, we introduce a new variant of the OPES method that focuses on quickly escaping metastable states at the expense of convergence speed.

In formulating OPES-explore, we restrict ourselves to the case of using as target the well-tempered distribution $$\begin{array}{c}{p}^{tg}(\mathbf{s})=p^{\text{WT}}(\mathbf{s})\propto [P(\mathbf{s}){]}^{1/\gamma }\end{array}$$

In OPES-explore, one builds the bias starting from the on-the-fly estimate of the distribution that is being sampled in the biased simulation: $$\begin{array}{c}{p}_n^{\text{WT}}(\mathbf{s})=\frac{1}{n}\sum _k^{n}G(\mathbf{s},{\mathbf{s}}_k),\end{array}$$ where ${\mathbf{s}}_k$ is the CVs value sampled at step $k$.

As the simulation converges, $p_n^{\text{WT}}(\mathbf{s})$ approaches the target well-tempered distribution $p^{\text{WT}}(\mathbf{s})$. Thus, we use the approximation $P(\mathbf{s})\propto [p_n^{\text{WT}}(\mathbf{s}){]}^{\gamma }$ and write the bias: $$\begin{array}{c}{V}_n(\mathbf{s})=(\gamma -1)\frac{1}{\beta }\log \left(\frac{p_n^{\text{WT}}(\mathbf{s})}{Z_n}+ϵ\right),\end{array}$$

Both OPES variants are applications of the general bias potential $$\begin{array}{c}V(\mathbf{s})=-\frac{1}{\beta }\log \frac{p^{tg}(\mathbf{s})}{P(\mathbf{s})}\end{array}$$ but OPES estimates on-the-fly $P(\mathbf{s})$ and uses it to calculate the bias, while OPES-explore does the same but with $p^{\text{WT}}(\mathbf{s})\propto [P(\mathbf{s}){]}^{1/\gamma }$.

The free energy surface as a function of the CVs can be estimated in two distinct ways, either directly from the probability estimate $$\begin{array}{c}{F}_n(\mathbf{s})=-\gamma \frac{1}{\beta }\log p_n^{\text{WT}}(\mathbf{s})\end{array}$$ or via importance sampling reweighting, e.g., using a weighted KDE, $$\begin{array}{c}{F}_n(\mathbf{s})=-\frac{1}{\beta }\log \sum _k^n{e}^{\beta V_{k-1}({\mathbf{s}}_k)}G(\mathbf{s},{\mathbf{s}}_k).\end{array}$$

Differences between OPES and OPES-explore

The idea of defining the bias potential is similar, i.e., using as target the well-tempered distribution $p^{tg}(\mathbf{s})=p^{\text{WT}}(\mathbf{s})\propto [P(\mathbf{s}){]}^{1/\gamma }$.

The way of estimating the probability distribution are different. Specifically, unbiased probability $P(\mathbf{s})$ is estimated on-the-fly using weighted KDE in OPES (NOTE: the bias is well-tempered in OPES), while the well-tempered (biased) probability $p^{\mathrm{WT}}(\mathbf{s})(\propto [P(\mathbf{s}){]}^{\frac{1}{\gamma }}$) is estimated based on averaged KDE in OPES-explore.

Therefore, OPES method focuses on fast convergence, while OPES-explore focuses on fast exploration.

Summary

📖 This Notebook summarizes the OPES methodology based on the author's reading of the original literature and practical experience.