Enhanced sampling of transition states

TL;DR

This paper introduces a novel enhanced sampling method based on variational techniques that dynamically targets transition states, improving sampling efficiency in free energy landscapes with high barriers.

Contribution

It proposes a dynamic target distribution using free energy derivatives to better sample transition regions in complex free energy landscapes.

Findings

Increased sampling of transition states in chemical reactions.

Effective enrichment of transition configurations in nucleation processes.

Demonstrated improvement over existing methods.

Abstract

The free energy landscapes of several fundamental processes are characterized by high barriers separating long-lived metastable states. In order to explore these type of landscapes enhanced sampling methods are used. While many such methods are able to obtain sufficient sampling in order to draw the free energy, the transition states are often sparsely sampled. We propose an approach based on the Variationally Enhanced Sampling Method to enhance sampling in the transition region. To this effect, we introduce a dynamic target distribution which uses the derivative of the instantaneous free energy surface to locate the transition regions on the fly and modulate the probability of sampling different regions. Finally, we exemplify the effectiveness of this approach in enriching the number of configurations in the transition state region in the cases of a chemical reaction and of a nucleation process.

Tables

Table 1. Table 1: Percentage of configurations observed in the TS region and the number of transitions between the two minimum energy states for 500 ps long simulations.

Simulation	Configurations in TS region (%)	Number of recrossings
Well tempered	19.02	229
Uniform	23.19	253
$p_L\left(s\right)$	28.38	233
$p\left(s\right)$	63.68	339

Equations

$$\Omega [V(\mathbf{s})]=\frac{1}{\beta }log\frac{\int d\mathbf{s}{e}^{-\beta (V(\mathbf{s})+F(\mathbf{s}))}}{\int d\mathbf{s}{e}^{-\beta F(\mathbf{s})}}+\int d\mathbf{s}p(\mathbf{s})V(\mathbf{s}),$$

$$V(\mathbf{s})=-F(\mathbf{s})-\frac{1}{\beta }\log p(\mathbf{s}).$$

$$P_V(\mathbf{s})=\frac{e^{-\beta (V(\mathbf{s})+F(\mathbf{s}))}}{\int d\mathbf{s}{e}^{-\beta (V(\mathbf{s})+F(\mathbf{s}))}},$$

$$\frac{\partial\Omega({\boldsymbol{\alpha}})}{\partial\alpha_{i}}=-\ \Big{\langle}{\frac{\partial V({\boldsymbol{s};\boldsymbol{\alpha}})}{\partial\alpha_{i}}}\Big{\rangle}_{V({\boldsymbol{\alpha}})}\ +\ \Big{\langle}{\frac{\partial V({\boldsymbol{s};\boldsymbol{\alpha}})}{\partial\alpha_{i}}}\Big{\rangle}_{p}\,,$$

$$V(\mathbf{s};\mathbf{\alpha })=\sum _k{\alpha }_k{f}_k(\mathbf{s}),$$

$$p(\mathbf{s})=\frac{e^{-\frac{\beta }{\gamma }F(\mathbf{s})}}{\int d\mathbf{s}{e}^{-\frac{\beta }{\gamma }F(\mathbf{s})}},$$

$$\begin{array}{rl}{F}^{(k+1)}(\mathbf{s})& =-V^{(k)}(\mathbf{s})-\frac{1}{\beta }\log p^{(k)}(\mathbf{s})\\ & =-V^{(k)}(\mathbf{s})-\frac{1}{\gamma }{F}^{(k)}(\mathbf{s}),\end{array}$$

$$p_L(s)=\frac{1}{Z_p}\frac{\zeta }{{\zeta }^2+{\left(\beta \frac{dF(s)}{ds}\right)}^2},$$

$$S_{a,\lambda }(s)=\frac{1-a}{1+e^{\lambda \beta \frac{d^{2}F(s)}{d{s}^2}}}+a,$$

$$\begin{array}{rl}p(s)& =\frac{1}{Z_p}\frac{\zeta }{{\zeta }^2+{\left(\beta \frac{dF(s)}{ds}\right)}^2}\left(\frac{1-a}{1+e^{\lambda \beta \frac{d^{2}F(s)}{d{s}^2}}}+a\right)\\ & =p_L(s)S_{a,\lambda }(s).\end{array}$$

$$\epsilon =\sqrt{\frac{\int ds[{\bar{F}}_r(s)-\bar{F}(s)]}{\int ds}}$$

$${\bar{F}}_r(s)=F_r(s)-\frac{\int ds{F}_r(s)}{\int ds}$$

$$\bar{F}(s)=F(s)-\frac{\int dsF(s)}{\int ds}$$

Full text

]Department of Chemistry and Applied Biosciences, ETH Zurich, c/o USI Campus, Via Giuseppe Buffi 13, CH-6900, Lugano, Switzerland \alsoaffiliationFacoltà di Informatica, Istituto di Scienze Computazionali, Università della Svizzera Italiana (USI), Via Giuseppe Buffi 13, CH-6900, Lugano, Switzerland

]Department of Physics, ETH Zurich, c/o USI Campus, Via Giuseppe Buffi 13, CH-6900, Lugano, Switzerland \alsoaffiliationFacoltà di Informatica, Istituto di Scienze Computazionali, Università della Svizzera Italiana (USI), Via Giuseppe Buffi 13, CH-6900, Lugano, Switzerland

\alsoaffiliationFacoltà di Informatica, Istituto di Scienze Computazionali, Università della Svizzera Italiana (USI), Via Giuseppe Buffi 13, CH-6900, Lugano, Switzerland \alsoaffiliationIstituto Italiano di Tecnologia, Via Morego 30, 16163 Genova, Italy

1 Introduction

One of the crucial problems in the atomistic simulation of systems of chemical-physical relevance is the study of rare events. Typical rare events are chemical reactions, phase transitions, and biomolecule conformational changes. These events take place on a time scale that cannot be reached by straightforward molecular dynamics (MD) and require enhanced sampling methods to be studied. A vast class of such methods 1, 2, 3, 4, 5, 6, 7 relies on the identification of an appropriate set of collective variables (CVs) whose fluctuations are then suitably enhanced. The CVs, here indicated as ${\boldsymbol{s}}={\boldsymbol{s}}({\boldsymbol{R}})$ , are functions of the atomic coordinates ${\boldsymbol{R}}$, and are meant to encode the slow degrees of freedom of the system. In our group we have focused on two such approaches, namely Metadynamics 2 and Variationally Enhanced Sampling (VES) 6. One of the outputs of these calculations is the free energy $F({\boldsymbol{s}})$ as a function of the chosen CVs. With adroitly chosen CVs, the $F({\boldsymbol{s}})$ separates the metastable states and the apparent TS does coincide with the real one. Luckily much progress has been made in finding 8 and ameliorating CVs 9, 10, 11 and we shall assume here that we are in a situation such that the CVs pass close to the real transition state (TS).

However, even with good CVs, in a standard Metadynamics or VES calculation the transition state is not as well sampled as one would like it to be. In fact in these methods based on a dynamical approach in which the $F({\boldsymbol{s}})$ is progressively filled by a bias, most of the sampling effort is spent in exploring the $F({\boldsymbol{s}})$ minima and the transition states (TS) are sampled only when the system moves from one basin to another. This is an event that, while accelerated by the bias, is comparatively rare. Given the relevance that the transition state plays in understanding the process under study and eventually guiding it towards a desired end, it would be of great help to have a method that partitions sampling more evenly between the different stationary points of $F({\boldsymbol{s}})$ and harvests a much larger set of transition state configurations.

In order to achieve this result we shall use one of the most attractive features of VES, namely the possibility of searching for a bias such that the distribution in the biased ensemble is equal to a pre-assigned one. This leads to an enhanced sampling of the transition state region whose properties can then be studied in great detail. We apply the method thus developed to the study of a bimolecular nucleophilic substitution (S${}_{\text{N}}$2) reaction and to a liquid-vapor transition, which is a prototypical nucleation process.

2 Method

2.1 Variationally Enhanced Sampling

In VES, the bias potential $V({\boldsymbol{s}})$ is obtained by optimizing the functional

[TABLE]

where $T$ is the temperature, $k_{B}$ is the Boltzmann constant, $\beta=1/k_{B}T$ and $p({\boldsymbol{s}})$ is a target distribution which is to be normalized. It has been shown 6 that $\Omega[V({\boldsymbol{s}})]$ is a convex functional and at its minimum the value for the bias $V({\boldsymbol{s}})$ is, within an irrelevant constant,

[TABLE]

From equation (2) it follows that at the minimum the CV distribution in the biased ensemble,

[TABLE]

is equal to the pre-assigned target distribution $p({\boldsymbol{s}})$.

Optimization

In order to use the variational principle and minimize the functional, the bias $V({\boldsymbol{s}})$ is expressed in terms of a set of parameters ${\boldsymbol{\alpha}}=\{\alpha_{1},\alpha_{2},\alpha_{3},...\}$. This turns the functional $\Omega[V({\boldsymbol{\alpha}})]$ into a function $\Omega({\boldsymbol{\alpha}})$. This function $\Omega({\boldsymbol{\alpha}})$ is then optimized using the gradient

[TABLE]

where the terms $\langle\cdots\rangle_{V({\boldsymbol{\alpha}})}$ and $\langle\cdots\rangle_{p}$ are expectation values calculated in the ensemble biased by ${V({\boldsymbol{\alpha}})}$ and over the distribution $p({\boldsymbol{s}})$ respectively. However, the gradient evaluated from such a statistical calculation is inherently noisy and thus, a stochastic optimization algorithm 12 is necessary to obtain a faster and smoother convergence.

Here, we take the bias potential $V({\boldsymbol{s}};{\boldsymbol{\alpha}})$ to be a linear expansion in a finite number of basis functions,

[TABLE]

and use the Bach and Moulines averaged stochastic gradient descent algorithm 13 to update iteratively the parameters of the bias expansion.

Target distribution

A straightforward choice of $p({\boldsymbol{s}})$ is the uniform distribution, that at convergence gives $V({\boldsymbol{s}})=-F({\boldsymbol{s}})$. However, this distribution forces the system to visit with equal probability all points in the CV space, including uninteresting high energy regions, often resulting in poor convergence. Although it is possible to construct a $p({\boldsymbol{s}})$ that focuses on important regions of the CV space, this requires apriori information on $F({\boldsymbol{s}})$, a quantity that is in general unknown at the beginning of the calculation. However, Valsson and Parrinello 7 introduced a self-consistent approach to include $F({\boldsymbol{s}})$ in the definition of $p({\boldsymbol{s}})$. They suggested a $p({\boldsymbol{s}})$ that at convergence results in,

[TABLE]

where $\gamma\geq 1$ is a parameter called bias factor. Unlike VES with a fixed $p({\boldsymbol{s}})$, this approach requires two coupled iterative schemes. In addition to updating the parameters for $V({\boldsymbol{s}};{\boldsymbol{\alpha}})$ every $\Delta t_{\text{V}}$ MD steps, the $p({\boldsymbol{s}})$ is updated every $\Delta t_{\text{TD}}$ number of iterations of the optimization of $\Omega({\boldsymbol{\alpha}})$. At first, $p^{(0)}({\boldsymbol{s}})$ is typically taken to be a uniform distribution and the first guess for $F({\boldsymbol{s}})$ is zero. The $F({\boldsymbol{s}})$ at any given step $(k+1)$ is then obtained as,

[TABLE]

In the asymptotic limit, one has $V({\boldsymbol{s}})=-(1-\frac{1}{\gamma})\ F({\boldsymbol{s}})$ modulo a constant as in the well-tempered ensemble.

It has been shown that, even if the instantaneous estimate of $F({\boldsymbol{s}})$ is initially inaccurate, VES quickly converges to the correct $F({\boldsymbol{s}})$ 6, 7. Additionally, in many practical cases, the well-tempered variant was found to converge much faster than a VES simulation that uses a uniform $p({\boldsymbol{s}})$. In figure 1, we illustrate with the case of a one-dimensional double-well potential of barrier height $40\ k_{B}T$, the probability density of the CV obtained with uniform and well-tempered $p({\boldsymbol{s}})$ as compared to the unbiased Boltzmann probability distribution. It can be seen that in a well-tempered simulation, the system visits the TS less frequently than in a uniform target simulation. However, as discussed above one has to pay the price of a slower convergence.

2.2 Targeting the transition state

For simplicity sake, we deal here only with cases in which the CV is unidimensional. However, an extension to a multidimensional set of CVs is possible. For a well chosen CV, the TS is located at the apparent $F(s)$ TS. Hence, sampling of the TS region can be enhanced by designing a $p({{s}})$ that gives a higher probability of sampling in a region around the apparent TS.

To do so, we utilize the instantaneous $F({{s}})$ similarly to what is done in well-tempered VES. We take into account the fact that the apparent TS is a $F({{s}})$ stationary point, where $\frac{dF({s})}{ds}=0$. A $p(s)$ constructed to enhance the probability of sampling the stationary points should, in principle, be able to boost TS sampling. The most obvious choice of such a distribution would be a Gaussian-like distribution similar to the well-tempered $p(s)$ that has as argument $\frac{dF({s})}{ds}$. However, we found it difficult to describe the rapidly varying tail of such a function with a finite basis expansion unless $\gamma$ was taken so large that the resulting distribution was practically uniform. Hence, we ameliorated the behaviour at the tail by using a Lorentzian-like function that decays polynomially

[TABLE]

where $Z_{p}$ is a normalization factor, and $\zeta$ is a parameter of units $s^{-1}$ that regulates the broadening of the distribution. This $p_{L}(s)$ will weight uniformly all the stationary states of $F(s)$: minima, TS, and if present, inflection points. However for the purpose of our calculation, it is desirable to give more weight to the TS. In order to do this we take advantage of the fact that at a local maximum not only $\frac{dF({s})}{ds}$ vanishes but also the second derivative satisfies the condition $\frac{d^{2}F({s})}{ds^{2}}<0$. For this reason, we multiply $p_{L}({s})$ in equation (8) by the function,

[TABLE]

where $0\leq a\leq 1$. This function goes from the value $S_{a,\lambda}(s)=1$, when $\frac{d^{2}F({s})}{ds^{2}}$ is large and negative to $S_{a,\lambda}(s)=a$, when $\frac{d^{2}F({s})}{ds^{2}}$ is large and positive. The parameter $\lambda$ determines the sharpness of this transition. In practice for $\lambda\to\infty$ and $a=0$, $S_{a,\lambda}(s)$ behaves like a rectangular pulse and has a non-zero value in the region in which $\frac{d^{2}F({s})}{ds^{2}}\leq 0$. This region from now on is taken to define the TS region. The case $a=1$ is included for generality and corresponds to using simply $p_{L}(s)$. This can come useful when an equilibrated sampling between maxima and minima is required. Thus, the resulting $p(s)$ is,

[TABLE]

In figure 2, we show the effect of changing the parameters $a$ and $\lambda$ in the symmetric double well case.

In order to obtain these results and the others reported below we had to solve the problem of computing the first and second derivatives of $F(s)$. If not done with care this can introduce unwanted noise. To this effect we expand $F({s})$ in the same basis set $\{f_{k}({s})\}$ used to express $V(s)$ (see equation (5)). In this way, the $s$ dependence appears only in the functions $f_{k}(s)$ whose derivatives can be computed analytically.

We now illustrate how this approach works in practice on the paradigmatic cases of a bimolecular nucleophilic substitution reaction and the condensation of a Lennard-Jones gas.

3 Bimolecular nucleophilic substitution reaction

Nucleophilic substitution reactions constitute a popular class of chemical reactions. The asymmetric substitution reaction of fluoromethane and chloromethane CH${{}_{\text{3}}}$F + Cl{{}^{\text{-}}}\ {\rightleftharpoons}\ CH${{}_{\text{3}}}$Cl + F${{}^{\text{-}}}$ is a simple example of such reactions and has been well studied over the years. Following a previous enhanced sampling study of the same reaction 14, we chose the difference between the C-Cl ($d_{Cl}$) and C-F ($d_{F}$) atomic distances as CV.

Here, we contrast four different simulations that use VES with different target distributions. The first is a well-tempered one with $\gamma=50$, the second uses a $p(s)$ that is uniform in the interval $-2.0\leq s\leq 2.9$. The third simulation uses $p_{L}(s)$ as defined in equation (8) with $\zeta=20$ Å*-1*. Finally, in the last one we use $p_{L}(s)$ given by $\zeta=40$ Å*-1* and multiplied by a $S_{a,\lambda}(s)$ that has $a=0.01$ and $\lambda=0.1$. In the last two cases the values of the parameters have been chosen such that the resulting $p(s)$ does not present very sharp features, thus limiting the number of basis functions that need to be included in the $V(s)$ expansion. Our choice of basis functions was to use orthonormal Legendre polynomials. A detailed description of the VES parameters and the simulation details is provided in the supplementary information (SI).

We first converge the bias coefficients (${\boldsymbol{\alpha}}$) and use the resulting static potential to perform the simulation. We then obtain the $F(s)$ using an umbrella-like reweighting. As expected all the four calculations give the same $F(s)$ within statistical errors. In table 1, we quantify the fraction of configurations in the TS region and the number of times the system makes a transition between the two minima. One can observe that the VES calculation with the $p_{L}(s)\ S_{a,\lambda}(s)$ target distribution has indeed resulted in a larger number of configurations sampled in the TS region. Furthermore, it can be seen from the table that the number of transitions is greatly enhanced.

4 Liquid-vapour phase transition

Nucleation is a phenomena of interest in a broad range of fields. A commonly used model in homogeneous nucleation studies is that of the Classical Nucleation Theory (CNT). It states that the energetics of growth of a stable nucleus (or cluster) in a metastable bulk is shaped by an interplay of two factors. While the creation of the interface between this cluster and the bulk phase hinders the cluster growth, the formation of these more stable regions favours it. These two terms are proportional to the surface area and the volume of the cluster respectively. As a result of these competing contributions, cluster growth is hindered until a critical size is reached. Beyond this critical size, the stability gained from the volume term surpasses the cost of interface formation and the cluster grows spontaneously. Due to the existence of this barrier, it is often difficult to determine the critical cluster size 1516 and obtain sufficient sampling in the region even using enhanced sampling techniques.

Here, we compare the results of our TS-target VES simulation with uniform VES simulation in the case of condensation of Argon. The details of the simulation can be found in the SI. Following Piaggi et al.17, we use a CV defined as $s=(n_{l})^{1/3}$ where $n_{l}$ is the number of liquid-like atoms. In the case in which there is only one cluster present $s$ is roughly proportional to the radius of the cluster. We further restrict the growth of clusters beyond $120$ atoms using a harmonic restraining potential having $\kappa=0.1$. The parameters for the uniform VES simulation are similar to those of Piaggi et al. 17 while that of the TS-target are $a=0.02$, $\lambda=0.1$ and $\zeta=30$. Since we want to focus predominantly in the critical region we use a very small $a$. Here also, we first converge the bias coefficients $({\boldsymbol{\alpha}})$ for the two cases using $10$ walker simulations and then perform a single walker simulation with a static bias potential. The $F(s)$ is obtained by an umbrella-like reweighting and is shown in figure 4.

Further, we analyze the shapes of the largest clusters obtained during the simulations. To do so, first we identify the largest clusters of liquid-like atoms obtained from both the simulations using the depth first search clustering algorithm described in Ref. 18. Clusters spanning more than half the box size and clusters with less than 6 atoms are rejected. We calculate the moment of inertia tensor matrix, ${\boldsymbol{T}}$ for all the clusters using $T_{\alpha\beta}=\sum^{\prime}m(x_{\alpha\beta}\ -\ \hat{x}_{\alpha\beta})$, where the primed sum is over the atoms contained in the largest cluster, $m$ and $x$ are the mass and position of these atoms while $\hat{x}$ is the position of the center of mass of the cluster. We diagonalize the matrix to obtain the three eigenvalues $(\lambda_{1}\geq\lambda_{2}\geq\lambda_{3})$ and subsequently calculate the shape anisotropy value defined as $k=\frac{3}{2}\frac{\lambda_{1}^{4}\ +\ \lambda_{2}^{4}\ +\ \lambda_{3}^{4}}{(\lambda_{1}^{2}\ +\ \lambda_{2}^{2}\ +\ \lambda_{3}^{2})^{2}}-\frac{1}{2}$. The value $k=0$ indicates a perfectly spherical cluster whereas the value $k=1$ corresponds to a cluster where all the particles are linearly arranged. Figure 5 shows how even clusters with the same number of particles could have very different shapes and hence, a range of anisotropy values. In figure 6, the probability density of the CV space is shown. The probability density has been estimated using multivariate gaussian kernels.

It can be concluded from figure 6 that the proposed target distribution has been able to sample more effectively the critical cluster size. With the uniform target distribution, the system tends to sample more the uninteresting region with very small clusters due to the low energy of the region while with the TS-target distribution, the high-energy region around the critical cluster has been sampled more in the same computational time. Also, it is evident from both the figures that as the cluster size increases the clusters tend to be more spherical. Hence, while the results obtained conform to previously observed trends, we are now able to quantify more confidently about the critical cluster properties.

Last but not least in this case an improvement in convergence rate is observed. In order to quantify the convergence of $F(s)$, we use an error metric similar to that used by Valsson and Parrinello 7.

[TABLE]

where

[TABLE]

[TABLE]

Multiple $F(s)$ are obtained by reweighting the Uniform and TS-target 10 walkers VES simulations at regular intervals using the reweighting scheme proposed in Ref. 19. The reference $F(s)$ is obtained by an umbrella-like reweighting of the simulation that used the converged static bias obtained from the 10 walker Uniform VES simulation. It can be seen from figure 7 that the TS-target does converges faster.

5 Conclusion

We have proposed a very general approach based on the Variationally Enhanced Sampling method to enhance sampling in the transition region, and exemplified with the cases of a S${}_{\text{N}}$2 reaction and condensation of Lennard-Jones gas how our method can be successfully applied to very diverse scenarios. Furthermore, our approach is capable of locating the transition states on the fly and acquiring an increased sampling of the transition region despite the process being a rare-event. This of course under the condition of dealing with a good quality CV. Taking into account the ability of our method to increase the frequency with which the transition region is visited, we believe that this could be especially useful in studying reactions that take place in condensed phases where the effect of the environment on the transition state is relevant and yet difficult to describe with a single conformation. Nucleation is one example of such class of reactions. We believe our approach in these cases may obtain a much faster $F(s)$ exploration. {acknowledgement}

The authors would like to thank Tarak Karmakar for carefully reading the manuscript and for useful discussions. The authors thank the European Union Grant No. ERC-2014-ADG-670227 and the National Center for Computational Design and Discovery of Novel Materials MARVEL for funding. The calculations were performed using the Euler HPC Cluster at ETH Zurich and the Mönch cluster of the Swiss National Computing Center (CSCS).

{suppinfo}

The computational details for the one-dimensional potential, the S${}_{\text{N}}$2 reaction and the condensation of Argon are included as Supporting Information.