Time-independent free energies from metadynamics via Mean Force Integration

TL;DR

This paper introduces Mean Force Integration (MFI), a novel method to compute time-independent free energy profiles from biased simulations, enhancing accuracy and efficiency in free energy calculations.

Contribution

MFI provides a new approach that avoids ensemble bias averaging and can integrate data from multiple simulations, improving free energy estimation in metadynamics.

Findings

Validated against analytical potentials showing accurate free energy profiles.

Successfully applied to conformational free energy landscape of ibuprofen.

Demonstrated ability to combine independent simulations for convergence.

Abstract

Inspired by thermodynamic integration, we propose a method for the calculation of time-independent free energy profiles from history-dependent biased simulations via Mean Force Integration (MFI). MFI circumvents the need for computing the ensemble average of the bias acting on the system c(t) and can be applied to different variants of metadynamics. Moreover, MFI naturally extends to aggregate information obtained from independent metadynamics simulations, allowing to converge free energy surfaces without the need to sample recrossing events in a single continuous trajectory. We validate MFI against one- and two-dimensional analytical potentials and by computing the conformational free energy landscape of ibuprofen in the bulk of its most common crystal phase.

Tables

Table 1. Table 1: Simulation settings used for the calculation of the FES of ibuprofen conformational isomerism using MFI.

Set	Initial state	Final State	Average Length
			[ns]
1	c1	any other	1.7
2	c2	any other	0.05
3	c3	any other	0.013
4	c4	any other	0.015
5	c5	any other	0.05
6	c6	any other	0.0008
7	c1	c4,c5 or c6	15
8	c4	c1,c2 or c3	0.13
9	c5	c1,c2 or c3	0.05
10	c6	c1,c2 or c3	0.009

Equations

$$F(\mathbf{s})=-{\beta }^{-1}\ln p(\mathbf{s})$$

$$F(\mathbf{s})=-{\beta }^{-1}\ln p^b(\mathbf{s})-V(\mathbf{s})-⟨V(\mathbf{s}){⟩}_u$$

$$⟨V(\mathbf{s}){⟩}_u={\beta }^{-1}\ln \frac{{\int }_{\Omega }{e}^{-\beta F(\mathbf{s})+\beta V(\mathbf{s})}d\mathbf{s}}{{\int }_{\Omega }{e}^{-\beta F(\mathbf{s})}d\mathbf{s}}$$

$$\frac{d{F}_t(s)}{ds}=-\frac{d{\beta }^{-1}\ln p_t^b(s)}{ds}-\frac{d{V}_t(s)}{ds}$$

$$\bigg{\langle}{\frac{dF_{t}(s)}{d{s}}}\bigg{\rangle}_{t}=\frac{\sum_{t^{\prime}=1}^{t}p_{t^{\prime}}^{b}(s)\frac{{dF_{t^{\prime}}(s)}}{d{s}}}{\sum_{t^{\prime}=1}^{t}p_{t^{\prime}}^{b}(s)}$$

$$\frac{d{V}_t(s)}{ds}=\sum _{t^{\prime }=1}^t-\frac{w_t(s-s_{t^{\prime }})}{{{\sigma }_{M,t}}^2}\exp \left[-\frac{1}{2}\frac{{\left(s-s_{t^{\prime }}\right)}^2}{{{\sigma }_{M,t}}^2}\right]$$

$$p_t^b(s)=\frac{1}{n_{\tau }h\sqrt{2\pi }}\sum _{t^{\prime }=t}^{t+\tau }\exp \left[-\frac{{\left(s-s_{t^{\prime }}\right)}^2}{2h^2}\right]$$

$$\frac{d{\beta }^{-1}\ln p_t^b(s)}{ds}=\frac{{\sum }_{t^{\prime }=t}^{t+\tau }-\frac{s-s_{t^{\prime }}}{\beta h^2}\exp \left[-\frac{{\left(s-s_{t^{\prime }}\right)}^2}{2h^2}\right]}{{\sum }_{t^{\prime }=t}^{t+\tau }\exp \left[-\frac{{\left(s-s_{t^{\prime }}\right)}^2}{2h^2}\right]}$$

$$\bigg{\langle}{\frac{dF_{t}(s)}{d{s}}}\bigg{\rangle}_{t}=\frac{1}{{\sum_{t^{\prime}=1}^{t}p_{t^{\prime}}^{b}(s)}}\left\{{\sum_{t^{\prime}=1}^{t}\sum_{t^{\prime\prime}=t^{\prime}}^{t^{\prime}+\tau}\frac{s-s_{t^{\prime\prime}}}{\beta\,n_{\tau}{h^{3}}\sqrt{2\pi}}\exp\left[{-\frac{\left({s-s_{t^{\prime\prime}}}\right)^{2}}{2h^{2}}}\right]}+{\sum_{t^{\prime}=1}^{t}p_{t^{\prime}}^{b}(s)\frac{{dV_{t^{\prime}}(s)}}{d{s}}}\right\}$$

$$⟨(\mathbf{R}){⟩}_u=⟨O(\mathbf{R})\exp \left[\beta V_t(\mathbf{s})-\beta ⟨V_t(\mathbf{s}){⟩}_u\right]{⟩}_t$$

$$\nabla F_t(\mathbf{s})=-{\beta }^{-1}\nabla \ln p_t^b(\mathbf{s})-\nabla V_t(\mathbf{s})$$

$$\begin{gathered} p_t^b(\mathbf{s})=\frac{1}{n_{\tau }{\left(2\pi \right)}^{d/2}|\mathbf{h}{|}^{1/d}}\times \\ \times \sum _{t^{\prime }=t}^{t+\tau }\exp \left[-\frac{1}{2}{\left(\mathbf{s}-{\mathbf{s}}_{t^{\prime }}\right)}^{⊺}{\mathbf{h}}^{-2}\left(\mathbf{s}-{\mathbf{s}}_{t^{\prime }}\right)\right] \end{gathered}$$

Full text

I Introduction

Several enhanced sampling molecular simulation methods are aimed at computing free energy surfaces (FES) as a function of low-dimensional sets of collective variables (CVs). A strategy common to several of these methods is the introduction of an artificial bias potential that perturbs the Hamiltonian of the system, enhances the sampling of rare transitions, and facilitates the exploration of high energy regions of phase space.chipot2007free ; torrie1977nonphysical ; marsili2006self ; grubmuller1995predicting ; voter1997hyperdynamics ; darve2001calculating ; Laio2002 ; Barducci2008 . Metadynamics Barducci2011 is a method that implements this concept by introducing a history-dependent bias potential, iteratively updated as a sum of Gaussian contributions defined in the space of CVs. As discussed in the original publications Barducci2008 ; Laio2002 and in several reviews on the topicBarducci2011 ; valsson2016enhancing ; laio2008metadynamics , the FES recovered from the metadynamics bias potential cannot be considered inherently time-independent. The time-dependence of a FES computed from metadynamics is captured by the work performed by the metadynamics algorithmTiwary2015A , usually indicated as the time-dependent constant $c(t)=\langle{V(\mathbf{s},t)}\rangle$Barducci2008 ; bonomi2009reconstructing ; Tiwary2015A . A popular approach at the calculation of time-independent free energy landscapes from time-dependent metadynamics simulations has been introduced by Tiwary and Parrinello Tiwary2015A , who proposed an explicit expression for $c(t)$, particularly suited to the analysis of well-tempered metadynamicsBarducci2008 calculations.

Here we propose a different approach to tackle this problem, based on the observation that while promoting the exploration of phase space, metadynamics probes the gradient of the free energy hypersurface projected in CV space. By analysing metadynamics simulations from this perspective we propose a method for the calculation of free energy surfaces through mean force integration (MFI). MFI applies to different variants of metadynamics, and provides a framework to consistently patch the sampling obtained from independent simulations in a single FES, enabling an efficient use of independent, asynchronous replicas. In this paper we discuss the MFI method, we validate it against model potentials and by computing the conformational free energy landscape of ibuprofen in the bulk of its crystal phase.

II Theory

The free energy profile along a suitably defined set of collective variables $\mathbf{s(R)}$, function of the atomic coordinates $\mathbf{R}$, can be expressed as:

[TABLE]

where $\beta=(k_{B}T)^{-1}$, in which $k_{B}$ is the Boltzmann constant and $T$ is the temperature; and $p({\mathbf{s}})$ is the equilibrium probability density projected on $\mathbf{s}$. Under the effect of a perturbation of the system’s Hamiltonian introduced by the bias potential $V(\mathbf{s})$, the unperturbed free energy profile $F(\mathbf{s})$ is Zwanzig1954 :

[TABLE]

Where $p^{b}(\mathbf{s})$ represents the equilibrium probability density under the effect of the bias potential $V(\mathbf{s})$, and $\langle{V(\mathbf{s})}\rangle_{u}$ is the ensemble average of the bias in the unperturbed ensemble:

[TABLE]

It should be noted that Eq. 2 provides an implicit expression for $F(\mathbf{s})$, which appears on the right hand side within $\langle{V(\mathbf{s})}\rangle_{u}$. In Eq. 2 the term $\langle{V(\mathbf{s})}\rangle_{u}$ is non-local, i.e. it contributes to the absolute value of $F(\mathbf{s})$, but the ensemble average operation makes it independent with respect to $\mathbf{s}$. The calculation of this term is essential in Umbrella Sampling (US) haydock1990tryptophan , for estimating free energy profiles from multiple biased simulations that sample different regions of $\mathbf{s}$. The estimate of this term in US simulations is commonly carried out iteratively via the Weighted Histogram Analysis Method (WHAM) algorithmroux1995calculation ; kumar1992weighted . In the context of adaptive enhanced sampling methods such as metadynamics, in which the bias potential changes in time according to the sampling history of the system, i.e. $V(\mathbf{s},t)$, the term $\langle{V(\mathbf{s})}\rangle_{u}$ is a function of time and corresponds to the constant $c(t)$ for which an explicit formulation has been proposed by Tiwary and Parrinello(Tiwary2015A, ). In the following we illustrate how MFI does not require an estimate of $c(t)$ to obtain a time-independent estimate of $F(\mathbf{s})$.

We shall begin by noting that in metadynamics the bias potential $V(\mathbf{s},t)$ is evolved discretely in time, through updates performed at regular time intervals of length $\tau$.

Between two consecutive updates of the bias potential, performed at times $t$ and $t+\tau$, the system evolves under the effect of the stationary bias $V_{t}(\mathbf{s})$ and samples the biased probability density $p_{t}^{b}(\mathbf{s})$. The estimate of $p_{t}^{b}(\mathbf{s})$ obtained during the sampling time $\tau$ is typically localised in a small subregion of $\mathbf{s}$, moreover in different iterations of the bias update algorithm, the perturbation of the Hamiltonian introduced by the biasing potential $V_{t}(\mathbf{s})$ is different. Hence, in order to reconstruct a global free energy surface $F(\mathbf{s})$ from Eq. 2, the term $\langle{V_{t}(\mathbf{s})}\rangle_{u}$ is necessary and has to be evaluated at every update of the bias potential.

In MFI we approach this problem by taking inspiration from the Umbrella Integration (UI) methodKastner2005 . In UI, instead of straightforwardly applying Eq. 2, the estimate of the non-local term $\langle{V(\mathbf{s})}\rangle_{u}$ is circumvented by computing the mean force in CV space, $\nabla{F(\mathbf{s})}$. The free energy surface $F(\mathbf{s})$ is then obtained by numerical integration of the mean force.

In order to gradually introduce complexity, in the following section we shall outline the details of the method for a mono-dimensional CV space. We then discuss the generalisation of MFI to CV spaces of higher dimensionality, and finally we outline how MFI provides the means to consistently merge the sampling obtained from independent simulations into a single estimate of the free energy surface.

Mean Force Integration in 1D CV spaces

For the sake of clarity, let us begin by considering a simple case in which $\mathbf{s}$ is a mono-dimensional CV space, thus indicated as the scalar $s$. The derivative of the free energy profile with respect to $s$ is:

[TABLE]

where $\frac{dF_{t}(s)}{d{s}}$ is the mean force in CV space obtained from the sampling performed in the time interval $[t;t+\tau]$, the term $\frac{d{V_{t}(s)}}{d{s}}$ is the derivative of the bias potential updated at time $t$, which is stationary during the time interval $[t;t+\tau]$. Finally, the term $\frac{d\beta^{-1}\ln{p_{t}^{b}(s)}}{d{s}}$ corresponds to the mean force in $s$ under the effect of the perturbation due to the bias potential $V_{t}(s)$, sampled during the time interval $[t;t+\tau]$. It should be noted that the term $\frac{d{V_{t}(s)}}{d{s}}$ is accumulated from all the updates of the bias potential performed up to time $t$. On the contrary, the term associated with $p_{t}^{b}(s)$ is estimated anew after every iterative update of the bias potential.

During a metadynamics simulation the bias is updated frequently, usually in thousands of iterations. Each update of the bias potential will yield a mean force estimate $\frac{dF_{t}(s)}{d{s}}$. Following the approach proposed in Umbrella Integration(Kastner2005, ), the average mean force realization at time $t$ is estimated as:

[TABLE]

From \big{\langle}{\frac{dF_{t}(s)}{d{s}}}\big{\rangle}_{t} a time-independent estimate of $F(s)$ is obtained through numerical integration. It should be noted also that we indicate with $\langle{...}\rangle_{t}$ the estimate at time $t$ of the mean force in $s$, however $\langle{\frac{dF_{t}(s)}{d{s}}}\rangle_{t}$ is an inherently time-independent quantity.

In order to apply Eq. 5, we derive an analytical expression for the terms $\frac{d{V_{t}(s)}}{d{s}}$ and $\frac{d\beta^{-1}\ln{p_{t}^{b}(s)}}{d{s}}$.

The former can be straightforwardly computed as the derivative of the sum of Gaussians accumulated up to time $t$:

[TABLE]

where $w_{t}$, and $\sigma_{M,t}$ are the values of the Gaussian height and width at time $t$, and $s_{t}$ is the position in CV space that corresponds to the mean value of the Gaussian deposited at time $t$.

As demonstrated in the results section, this expression holds regardless of the protocol followed to update $V_{t}(\mathbf{s})$, and is therefore applicable to any metadynamics variant including standard MetaD ($w_{t}$, and $\sigma_{M,t}$ constant), WTmetaD and TTmetaD Dama2014 ($w_{t}$, updated at every iteration, $\sigma_{M,t}$ constant), and adaptive Gaussians metaDBranduardi2012 ($\sigma_{M,t}$ updated at every iteration).

In order to express the term $\frac{d\beta^{-1}\ln{p_{t}^{b}(\mathbf{s})}}{d{s}}$ in a general form, we apply a kernel density estimation of $p_{t}^{b}(\mathbf{s})$, the biased probability density sampled in the time interval $[t;t+\tau]$. Using Gaussian kernels $p_{t}^{b}(\mathbf{s})$ takes the form:

[TABLE]

where $n_{\tau}$ is the number of frames sampled in the time interval $[t;t+\tau]$, $h$ is the kernel bandwidth, $s_{t}$ is the instantaneous value of $\mathbf{s}$. Thus, the mean force contribution associated with the biased probability density term is:

[TABLE]

We note that in order to apply MFI to metadynamics we need to significantly depart from the hypothesis of Kastner et al. in Umbrella Integration, i.e. of $p_{t}^{b}(\mathbf{s})$ being a mono-modal probability density, normally distributed around the average value of $\mathbf{s}$. This hypothesis holds for the Umbrella Sampling protocol, where the bias potential confines sampling in a specific region of $\mathbf{s}$, localised around a certain target position. However it breaks down in the case of metadynamics, where each repulsive Gaussian contribution tends to push the system away from its center. This yields biased distributions that are far from being mono-modal, even on the short timescale of $\tau$, which can nevertheless be faithfully captured by Eq. 7.

Combining Eq. 4-8 we obtain an analytic expression for the mean force in $\mathbf{s}$:

[TABLE]

A graphical scheme representing the calculation procedure for the update of the mean force through Eq. 9 is reported in Fig. 1. It should be noted that, while the second term of Eq. 9 depends on the specific bias protocol, the first term is generally valid for any history-dependent biasing protocol based on discrete iterative updates of the bias potential. By numerically integrating \big{\langle}{\frac{dF_{t}(\mathbf{s})}{d{\mathbf{s}}}}\big{\rangle}_{t} one can obtain a time-independent estimate of the free energy surface $F(\mathbf{s})$. In the long time limit $p_{t}^{b}(\mathbf{s})$ approaches the limit distribution associated with the chosen sampling method. In the case of standard metadynamics, in the long time limit the biased distribution becomes flat, and the term $\frac{d\beta^{-1}\ln{p_{t}^{b}(\mathbf{s})}}{d{\mathbf{s}}}\rightarrow{0}$, thus recovering the standard estimator of the free energy $F(\mathbf{s})=-V(\mathbf{s})+C$. If needed, by using the integrated profile $F(\mathbf{s})$ in Eq. 3 one can compute the non-local time dependent constant $c(t)=\langle{V(\mathbf{s},t)}\rangle_{u}$ and perform on-the-fly reweighting for additional variables function of the system’s coordinates $O(\mathbf{R})$ as:

[TABLE]

where, following the notation of Ref.Tiwary2015A , with angular brackets on the right-hand side we indicate an average over the biased simulation. In contrast to the approach of Tiwary and Parrinello we do not invoke any assumption on the bias evolution, and Eq. 10 is valid for any biasing protocol, as long as the bias is updated at discrete time intervals of length $\tau$ that enable the local estimate of $p_{t}^{b}(\mathbf{s})$ under the effect of a stationary bias $V_{t}(\mathbf{s})$.

Generalization to a $d$-dimensional free energy surface

The result obtained for the monodimensional case can be straightforwardly generalised to the calculation of an arbitrary-dimensional free energy hypersurface $\mathbf{s}$, where Eq. 4 becomes Kastner2009 :

[TABLE]

In this case, the $d$-dimensional biased probability distribution $p_{t}^{b}(\mathbf{s})$ sampled in the time interval $\tau$ can be obtained with a multivariate kernel density estimation. Using multivariate Gaussian kernels $p_{t}^{b}(\mathbf{s})$ takes the form:

[TABLE]

where $d$ is the dimensionality of the CV space $\mathbf{s}$, $\mathbf{h}$ is the variance/covariance matrix of the multivariate Gaussian kernel, and $|{\mathbf{h}}|$ is the determinant of the variance/covariance matrix. For the typical case of $d$=2, and diagonal $\mathbf{h}$ we have that $|{\mathbf{h}}|^{1/d}=\sqrt{h_{1}h_{2}}$, and $\left(\mathbf{s}-\mathbf{s}_{t^{\prime}}\right)^{\intercal}\mathbf{h}^{-2}\left(\mathbf{s}-\mathbf{s}_{t^{\prime}}\right)=\sum_{i=1}^{2}\left(\frac{s_{i}-s_{i,{t^{\prime}}}}{h_{i}}\right)^{2}$ where $h_{1}$ and $h_{2}$ are the bivariate Gaussian kernel bandwidths in dimensions 1 and 2 respectively.

In a $d$-dimensional metadynamics simulation the bias $V_{t}(\mathbf{s})$ is represented as the sum of multivariate Gaussian contributions in $d$ dimensions. It should be noted that while in standard and well-tempered metadynamics the variance/covariance matrix is kept constant, in the case of adaptive Gaussians metaD, the variance/covariance matrix of the multidimensional Gaussian kernels is not diagonal and its terms are adaptively estimated based on sampling.

III Results and Discussion

In this section we assess the accuracy of MFI by computing the free energy surface in the case of one and two dimensional model potentials using metadynamics. After that we demonstrate, for the case of ibuprofen conformational isomerism in its crystal bulk, how MFI enables the efficient and accurate calculation of free energy surfaces by merging the information obtained from independent simulations that do not include recrossing events.

All metadynamics simulations analysed in this section, both on analytical model potentials and on ibuprofen, were performed using PLUMED 2.4 tribello2014plumed . In the case of ibuprofen GROMACS 5.1.4 VanDerSpoel2005 was used as a molecular dynamics engine with the Generalised Amber Force Field (GAFF) used to define the potential energy of ibuprofen Wang2004 ; Marinova2018 . A detailed description of the simulations setup for ibuprofen conformational isomerism in the crystal bulk can be found in Supplementary Materials, as well as in Ref. Marinova2018 . All the data and PLUMED input files required to reproduce the results reported in this paper are available on PLUMED-NEST (www.plumed-nest.org), the public repository of the PLUMED consortium NEST , as plumID:19.071.

Model Potentials

1D double well.

In order to quantify the accuracy and the convergence rate of the free energy surface obtained by MFI we start from a simple 1D double well model potential. We perform Langevin dynamics, biased with different metadynamics variants, i.e. standard, well-tempered, transition-tempered and adaptive Gaussians.

In Fig.2(a-c) we report, as a function of the number of Gaussians added to define the bias potential, $n_{G}$, the quantity $\overline{\epsilon{(s)}}\times{\sqrt{n_{G}}}$. Where $\overline{\epsilon{(s)}}$ is the average absolute error in the estimate of the free energy profile $F(s)$. The exact free energy profile is defined by the expression $F_{exact}(\mathbf{s})=-5s^{2}+s^{4}$. In the absence of systematic errors $\overline{\epsilon{(s)}}\times{\sqrt{n_{G}}}$ converges to a flat plateauBarducci2008 . As can be seen in panels $a$-$c$ such condition, is met in plain metadynamics ($a$), Well Tempered Adaptive Bias metadynamics ($b$), and Well Tempered metadynamics ($c$). The relative position of the horizontal plateau for the quantity $\overline{\epsilon{(s)}}\times{\sqrt{n_{G}}}$ allows to compare the relative rate of convergence of different methods. For instance, the lower the plateau, the faster the convergence. From panels $a$ and $b$ one can see that MFI (red) provides faster convergence compared to the standard estimator based on the total bias deposited (gray) for both standard metadynamics and adaptive bias metadynamics. Similar results are reported in the Supplementary material for Well Tempered and Transition Tempered Metadynamics.

Focussing on Fig.2(a,b), and in particular on the free energy profiles obtained at $t_{1}$ and $t_{2}$ (indicated with dashed lines in Fig. 2a) and b), one can see that the estimate of the double-well FES obtained from MFI (solid red curve) provides a more accurate estimate of the analytical FES (dashed blue line) than that obtained from the bias potential (solid gray curve) after the same number of bias updates. In shaded red (MFI) or gray (bias) is represented the position-dependent absolute error, demonstrating that the error associated with MFI is significantly smaller than the associated with FES estimates based only on the Bias potential.

In Fig. 2c we report for Well Tempered metadynamics, a comparison between the convergence rate of MFI and three typical reweighting methods. We consider the reweighting approaches proposed by Tiwary and ParrinelloTiwary2015A , Bonomi et al. bonomi2009reconstructing , and a simplistic reweighting strategy based on considering all the sampling during the WTmetaD as if it was performed under the effect of the final bias.Gimondi2018 It can be seen that, like MFI, both the reweighting methods of Bonomi et al. bonomi2009reconstructing , and Tiwary et al. Tiwary2015A converge faster than estimates of the FES based solely on the bias potential. Among all methods compared in Fig. 2c, MFI shows the fastest convergence rate.

The analysis of the results obtained for this initial test case shows that MFI provides an accurate estimates of the model FES, which rapidly converges to the analytical results once a crossing event is sampled (i.e. around $n_{G}$=200 in Fig. 2a-c). We note finally that the rate of convergence and the quality of the FES obtained with MFI are very moderately affected by the choice of particular metadynamics variant, and that in the case of WT metadynamics, where several reweighting strategies can be comparedTiwary2015A ; bonomi2009reconstructing ; Gimondi2018 , MFI display the fastest asymptotic convergence rate.

2D double well.

In order to further demonstrate the applicability of MFI to 2D surfaces and assess its sensitivity to key parameters appearing in the definition of the thermodynamic force expression reported in Eq. 11, we perform Langevin dynamics simulations on a 2D double-well model potential ($F_{exact}(\mathbf{s})=-3s_{1}^{2}+s_{1}^{4}-3s_{1}s_{2}+s_{2}^{4}$).

In Fig. 3a) and c) we report the FES obtained with MFI and compare it to the FES obtained as $F(\mathbf{s})=-V(\mathbf{s})$ at the same simulation time. It can be seen that, the FES obtained through MFI provides a better representation of the exact analytical potential given the same sampling. This is particularly evident in the transition region between the two local minima. In Fig. 3b) we show the time evolution of the quantity $\overline{\epsilon{(s)}}\times{\sqrt{n_{G}}}$ demonstrating that MFI does not display any systematic error, and that it converges faster than the negative of the bias potential to $F_{exact}$. The reason for the faster convergence is that, similarly to what happens to the method proposed by Tiwary et al., MFI provides an expression for the mean force that holds for any time, and not only in the long time limit. Hence realisations at short times, when the bias potential deposited in the CV space between minima is scarce, are already representative of the exact FES.

In Fig.3d-f), we conduct a systematic investigation of the error associated to the parameters that can be freely selected to inform Eq. 9. Such parameters are: i) the bandwidth of the Gaussian kernels used to construct the biased probability density $p_{b}(\mathbf{s})$ (Fig.3d), ii) the stride used to extract data points from the system’s evolution to build a kernel density estimator of $p_{b}(\mathbf{s})$ (Fig.3e), and iii) the number of grid points used to numerically integrate the mean force in CV space, and obtain an estimate of the FES (Fig.3f).

In all cases the dependence of the error on parameters is weak, with a mean absolute error of the order of $k_{B}T$ over the entire of the parameter space investigated. Nevertheless, the dependence of the error on each of the three parameters is different. For instance, the error dependence on the number of grid-points and on M the stride display trends typical of numerical convergence, in which the error decreases monotonically with a smaller stride (i.e. more data points) and a finer grid used for numerical integration. The error dependence on the bandwidth used to estimate $p_{b}(\mathbf{s})$ instead displays a non-monotonic behaviour, which for the 2D Langevin simulation has a minimum for a bandwidth comprised between 0.1 and 0.2. This range corresponds to half of the mean fluctuation of the CV in the time interval $\tau$, that separates successive updates of the metadynamics potential. This observation confirms the validity of the heuristics typically implemented in the selection of the Gaussian width in setting up a metadynamics simulation also in defining a sensible bandwidth for the calculation of the FES through MFI.

Patching independent metadynamics simulations: the case of ibuprofen.

To demonstrate the effectiveness of MFI in making the most of the sampling carried out by independent simulations, we compute the FES associated with ibuprofen conformational rearrangements in the crystal bulk. Ibuprofen is a small organic molecule consisting of a phenyl ring with two para-substituents and a chiral centre. In a recent workMarinova2018 , we have demonstrated how, unlike in solution, when ibuprofen is embedded in the crystal bulk its conformational rearrangement is restricted and the escape from its crystallographic conformational state (named c1 in Fig. 4a) is a rare event with an associated timescale of around 100 ns.

In this work, we compute the free energy landscape associated to ibuprofen in the crystal bulk by performing a series of independent metadynamics simulations initialized in each of its conformational states. Simulations are stopped once a prescribed crossing event is observed. The sets of simulations, together with their initial and final configuration, and average duration are reported in Tab. 1. Additional details on the configurational landscape and the simulation setup are provided as Supplementary Material.

In order to compute the sampling error associated to MFI we divide the simulations in six groups, each containing five randomly selected simulations from every set reported in Tab. 1. In Fig. 4a) the FES (blue to red colormap) obtained by averaging the results of each group is reported. The FES is represented in the space of two torsional angles. A global one, $s_{1}$, which describes the rearrangement of the para-substituents of the phenyl ring, and a local one, $s_{2}$, capturing the rotation of the methyl groups within the isobutanyl substituentMarinova2018 . Starting from the crystallographic conformer, c1, the rotation of the local torsional angle generates isomers c2 and c3. Rotation along the global torsional angle in each of the conformers c1, c2 and c3 results in respectively conformers c4, c5 and c6 as shown in Fig. 3 in the SI. In Fig. 4a) we also report, as a term of comparison, the isocontours of a reference FES obtained with standard post-processing of a 120 ns long WTmetaD simulation performed as reported in the SI. The position-dependent standard error in the FES computed with MFI is generally rather small, as shown in Fig. 4b), with a maximum of 0.5 kJ/mol and an average of $\sim 0.1$ kJ/mol in the region of interest ($\Delta{F}<45$ kJ/mol).

We note that the transition pathways between states need to be sufficiently sampled to obtain an accurate estimate of the free energy difference between them. In Fig. 4c) we show the sampling achieved when using simulation sets 1 to 6, which correspond to simulations that are stopped as soon as the starting configuration transforms in any of the other stable conformers. On the plot each set is represented in a different colour, according to its starting configuration. These sets of simulations allow an accurate reconstruction of the free energy profile for the CV space occupied by conformers c1, c2 and c3 as the sampling in the channels between them is sufficient to connect the corresponding regions of the free energy, but are insufficient to generate the full FES accurately.

By including sets 7 to 10 the sampling in the transition channels c2 $\Longleftrightarrow$ c4 and c5 $\Longleftrightarrow{}$ c1 is improved by forcing a cross over along the $s_{1}$ direction at shown in In Fig. 4d) yielding a fully converged FES.

By employing MFI we have successfully reconstructed an accurate FES associated with the conformational rearrangement of ibuprofen in the crystal bulk from independent simulations without recrossings. MFI proves to be a powerful tool in obtaining free energy profiles without the need of sampling recrossing transitions along the same, continuous trajectory.

IV Conclusions

In this work, we have introduced MFI as a method for the calculation of time-independent free energy surfaces from history-dependent metadynamics simulations. Inspired by Umbrella integration, MFI is based on the analytic evaluation of the mean force in CV space and does not require the explicit calculation of the ensemble average of the deposited bias $\langle{V(\mathbf{s})}\rangle$. MFI applies to any history-dependent biasing schedule, provided that the bias is updated in discrete time steps, separated by a time interval $\tau$. We have shown that MFI provides accurate and rapidly converging estimates of analytical free energy profiles in one and two dimensions. Furthermore, we have demonstrated the applicability of MFI to the calculation of free energy surfaces from ensembles of independent metadynamics trajectories without the requirement of sampling recrossing events within the same continuous trajectory. We anticipate that MFI will be useful to obtain and systematically refine FES when metadynamics trajectories independently sampling realizations of a rare event are available. An example of this application would be the recovery of FES from ensembles of biased trajectories generated with an infrequent metadynamics protocol.salvalaglio2014assessing ; tiwary2013metadynamics

Supplementary Material

Convergence for 1D model potential in the case of Well Tempered and Transition Tempered Metadynamics, additional information on the Ibuprofen simulation setup, and on the associated free energy landscape, including convergence and error analysis.

Acknowledgements

This work was financially supported by Pfizer, and by the Engineering and Physical Sciences Research Council (EPSRC) grant EP/R018820/1. We acknowledge the Legion High Performance Computing Facility for access to Legion@UCL and associated support services in the completion of this work.

V References