Guidelines for the analysis of free energy calculations

1.1 Free energy calculations assist drug discovery

Complex chemical and biological systems pose a key challenge for modern molecular and computational science. We seek computational models which can provide quantitative predictions, not just qualitative insight. Researchers seek to answer questions such as “how much?”, “how big?”, “how tight?” and so on, and increasingly apply physically-detailed computation to help answer these questions. Models seek to mimic or simulate the processes in question, helping reveal and provide new understanding of mechanisms and phenomena which might be challenging or impossible to probe experimentally [20, 23, 38, 43].

Free energy calculations [11, 7, 25, 9, 18] provide a good example of a computational technique which provides a quantitative answer to a specific question – in this case, “what is the free energy difference between the two thermodynamic end states of the system?” This question arises, for example, in drug discovery [16] where drugs need to be ranked by, among other criteria, their binding affinity [9] to a target protein. To use free energy calculations to answer this question, one must first build a model of the system, and then identify the end states between which the free energy difference is to be computed.

1.2 Free energy calculations begin with a definition of the end states

The thermodynamic end states (for example, Fig. 1 states 1 and 2) are the key starting point in free energy calculations. In principle, free energy differences between the end states can be computed simply from simulations conducted in one or both states [11]. But in practice, this is typically not possible for biomolecular systems on reasonable timescales. To compute accurate free energy differences between states, their phase space integrals must have sufficient overlap which in practice is attainable only when both states are extremely similar. When this is not the case, it can be impossible to directly compute the free energy difference between end states 1 and 2. In such cases, we can instead compute free energy differences between a series of intermediate states which do have sufficient overlap, leading from state 1 to state 2. These intermediate states are typically artificial, unphysical states constructed to link the physical states of interest, and form part of a thermodynamic cycle (see Fig. 1) linking the two end states of interest. For most free energy calculations relevant to binding, solvation, and solubility, this alternate, unphysical pathway involves effectively deleting and/or inserting some atoms, while possibly also making parameter changes to those and other atoms.

1.3 The thermodynamic cycle depicts alternate paths between the end states

The thermodynamic cycle for standard hydration free energy calculations is comprised of the four legs joining the four states of interest (Fig. 1): (2) a molecule interacting with a box of water, (1) the same molecule present alone in the gas phase, (4) the molecule in water but not interacting with the surrounding water, and (3) the non-interacting molecule again alone in the gas phase. We thus compute the solvation free energy by modifying the solute molecule in each of its environments. To do so, we compute the free energy of turning off the solute’s non-bonded interactions with its environment (called decoupling) or turning off both internal non-bonded interactions and interactions with the environment (called annihilation) (Fig. 1)¹. Specifically, we compute the free energies associated with Figure 1, 1 → 3 and 2 → 4, turning off the molecule’s interactions in gas phase and in water, respectively. Most commonly, both of these transformations are carried out by first scaling (usually linearly) the solute charges to zero and then turning off the solute’s Lennard-Jones (LJ) interactions (usually via the “soft-core” scheme [2]). These transformations are done in a series of steps by introducing a parameter λ which modulates the potential energy of the system, so that as λ goes from 0 to 1 the potential energy transitions between that of the initial state and that of the target final state. Simulations are then run at a set of different λ values connecting the two states. In other words, each of the two transformation pathways is subdivided into a variety of individual steps, where each step involves a transition between two λ values. The number and spacing of λ values is chosen to assure adequate overlap between the conformational spaces of the two states being considered. These intermediate states, and their corresponding λ values, are normally said to be alchemical, since they correspond to unphysical states, often involving a change in the chemical identity of the species considered.

Following alchemical transformation of the molecule in gas and solution, it remains to connect the two end states ((3) and (4) in Fig. 1). However, the free energy of the non-interacting molecule does not depend on the nature of its environment, and so the transfer free energy of the non-interacting molecule between environments is 0. Thus, there is no associated free energy change going from states (3) to (4) in Fig. 1.

1.4 The interactions can be decoupled

Changes in electrostatic interactions are often separated from changes in LJ interactions to avoid inaccuracy in the free energy estimate and sampling challenges. Specifically, if electrostatic interactions are retained while an atom’s LJ interactions are being removed, the associated charge becomes more and more exposed, and can create huge electrostatic forces leading to large and expensive-to-converge free energy differences. In extreme cases this can result in the lack of separation between positive and negative charges, which is especially problematic, potentially leading to numerical instabilities and simulation crashes [3, 31]. An additional benefit of separating these transformations is it provides a mechanism to maintain optimally efficient linear scaling of the charge interactions [27] with λ while using alternative scaling schemes for LJ interactions. Specifically, to avoid situations when the derivative ∂U/∂λ (needed for the TI analysis) would have been discontinuous, the potential for LJ transformations is typically treated via the soft-core potential [2, 42, 36, 39]. Although we focus here on the decoupled scheme, our general analysis would apply to the combined case (the one that would require no electrostatic decoupling [5, 12]) as well.

1.5 The intermediate lambda states can be controlled by lambda vectors

The MD packages like GROMACS [33] and DESMOND [37] can handle free energy calculations via multiple λ values controlling progress of different interaction types so that, for example, Coulomb, LJ, and restraining transformations can be controlled separately. Each step along the transformation path is associated with a unique set of λ values that is often referred to as the λ vector. For the thermodynamic cycle in Fig. 1 the λ vector has two components that control the Coulomb and LJ interactions. Each of the transformation paths 1 → 3 and 2 → 4 can then be presented as a train of intermediate coupled states (λ_coul, λ_LJ), the initial and final states being (0, 0) and (1,1) as depicted in Fig 2. If λ controls the strength of Coulomb and LJ interactions with the solute’s environment, then as λ progresses, solute-solvent interactions gradually decrease until, at the end state, the system consists of pure solvent, overlaid by a parallel system consisting of the non-interacting, isolated solute with full internal interactions. In calculations where internal solute non-bonded interactions are removed as well, the solute end state is slightly different, consisting of an assembly of atoms which interact only via their bonded interactions.

Once λ states are selected, equilibrium simulations are carried out, storing the necessary information for analysis (typically ∂U/∂λ, the derivative of the potential with respect to λ, and ΔU_i,j the potential energy differences between states at the different λ values evaluated from individual trajectories).

1.6 The automated analysis should be an essential part of the free energy calculations

Free energy calculations have typically been an experts-only endeavor, and one reason for this is that both their setup and analysis require, as a rule, substantial manual intervention. Analysis often involves in-house scripts and as more researchers get involved with the free energy calculations, standard analysis tools become increasingly important both to help ensure best practices are followed, and to avoid duplication of effort.

Our focus here is on the analysis of free energy calculations, which typically consists of a series of sequential steps. These free energy calculations themselves can be conducted with a variety of different sampling techniques, and our focus here is primarily on the analysis stage, regardless of sampling technique.

Here we present what we believe are current best practices for analysis of alchemical free energy calculations. Conceptually, we break analysis into four main stages:

1. subsampling the data to retain uncorrelated samples

2. calculating free energy differences along with the corresponding statistical errors via a variety of TI-and FEP-based methods

3. producing textual and graphical outputs of the computed data

4. inspecting for

   • –

     for convergence and identifying the equilibrated portion of the simulation

   • –

     good phase space overlap for all pairs of adjacent lambda states

2.1 Obtaining input data

As noted, the key input information needed for full, general analysis of free energy calculations includes potential energy differences between (at least) adjacent lambda values, as well as ∂U/∂λ values at all lambda values. To be specific, we will give an example of the calculation with GROMACS-formatted input files. In particular, GROMACS currently (v3.3 through v5.0) store all energies to binary energy files, but also write out all the potential energies and differences thereof needed for analysis to human-readable text files with the .xvg file format. In GROMACS, these are formatted as shown in Fig. 3. For standard simulations (in contrast to expanded ensemble simulations discussed below), there are several such files-one for each λ value. In GROMACS, the precise number of ∂U/∂λ fields varies with the number of different types of λ value which are utilized, which corresponds to the number of dimensions in the λ vector.²

Expanded ensemble simulations [22, 21, 24, 29, 26] are an approach which allows for simultaneous exploration of both λ and coordinate space in a single simulation, potentially allowing for faster sampling across alchemical states provided that the kinetic barriers that divide conformations important are lower at other λ states. In expanded ensemble simulations, a single simulation samples all states, and thus produces a single energy file. At each time step, the simulation is in one specific λ state, which is stored in the second field of the output energy file only in the expanded ensemble case. This allows determination of which λ state stored ∂U/∂λ and ΔU_i,j values belong to.

2.2 Alchemical analysis techniques can be divided into families

As noted, a variety of methods can take output from alchemical free energy calculations and yield free energy differences. Conceptually, these methods can be divided into two categories based on the quantity used to compute ΔG: thermodynamic integration [17] (TI) methods and free energy perturbation [44] (FEP) methods.

In TI, the free energy change along the path composed of K states is computed as a weighted sum of the ensemble averages of the derivative of potential energy function with respect to the coupling parameter λ:

where W_i are the weighting factors that depend on the numerical integration scheme used [28].

Several different schemes are available for numerical integration in TI. In our provided tool, alchemical-analysis.py, we implement TI-1 and TI-3 [28] which differ in how they interpolate between data points for integration. TI-1 uses the trapezoidal rule (a first-order polynomial), while TI-3 uses a (natural) cubic spline. The relative performance of these different TI methods will depend on the nature of the underlying data and the shape of the $\frac{\partial U}{\partial \mathrm{\lambda }}$ curve being integrated – and thus, it depends on the alchemical path chosen.

Perturbation-based methods include a broad range of techniques loosely related to FEP. In our prototype tool, these include Deletion Exponential Averaging (DEXP), Insertion Exponential Averaging (IEXP), Gaussian Deletion (GDEL), Gaussian Insertion (GINS), Bennett Acceptance Ratio (BAR), Unoptimized Bennett Acceptance Ratio (UBAR), Range-based Bennett Acceptance Ratio (RBAR), and Multistate Bennett Acceptance Ratio (MBAR). Some of them, like BAR [1] and MBAR [35], are in common use and are deeply entrenched in the parlance of the field, while others either do not have customary, generally accepted names, or remain little-known. For these methods we will use naming conventions suggested by Paliwal and Shirts [28]. The first two, DEXP and IEXP, are based on the exponential averaging scheme of the potential energies between two adjacent states—the so-called Zwanzig relationship [44]:

Depending on the direction of the transformation the process can be interpreted as either “deletion” or “insertion”, hence the first letters in the acronyms. Typically, one of these processes, DEXP, proceeds in the direction of increasing entropy, while the other, IEXP, proceeds in the direction of decreasing entropy.

If potential energy differences are distributed in a Gaussian manner (GINS and GDEL [28]), then the Zwanzig relationship reduces to [14, 13]:

where ${\sigma }_{\Delta U_{ij}}^2=\langle \Delta U_{ij}^2\rangle -{\langle \Delta U_{ij}\rangle }^2$ Again, here, the estimator is referred to as either GDEL when ΔU_ij is used in the direction of increasing entropy, or GINS when ΔU_ij is used in the direction of decreasing entropy.

Methods based directly on the Zwanzig relationship, such as those just discussed, yield alternate estimations for the free energy difference depending on the direction of the transformation. These discrepancies originate from undersampling in the tail regions of the ΔU_ij distributions [32], which yields biased free energy estimates. BAR [1] eliminates the bias in ΔG estimation by including both forward, ΔU_ij, and reverse, ΔU_ji, potential energy differences in the analysis. In BAR, the free energy change between intermediate states i and j (comprised of N_i and N_j microstates, respectively) is found by solving numerically the implicit function of ΔU_ij:

where $C=\Delta G_{ij}+\frac{1}{\beta }\ln \left(\frac{N_j}{N_i}\right)$.

In addition to this full version of BAR, there are two BAR-related methods that are advantageous in that there is no need to retain all potential energy differences for postprocessing. These methods focus on accumulating the averages in eq 4 as the simulation progresses. This is achieved by either setting the constant C = β⁻¹ ln(N_j/N_i) and thus avoiding the self-consistency procedure entirely (UBAR [28]), or picking a range of starting values of C and obtaining a range of ΔG_ij estimates from which the one having minimum variance is chosen as an input value for the constant C thus making the self-consistent solution essentially precalculated (RBAR [28]). UBAR can have issues when the free energy is significantly different from zero, while RBAR is essentially as accurate as BAR as long as the true value of C is within ≈ 1 – 2k_BT of the one of the range of trial C values.

MBAR [35] constitutes a further development of the BAR method. In BAR, the free energy change between the two adjacent states is computed to yield the minimum variance given data collected at that single pair of states alone, while MBAR finds the best estimate of free energy changes between all states simultaneously by optimizing the matrix of the ΔG variances, thus making use of all available data. MBAR can also be considered [35] as a limiting case of the WHAM [19] method in which the histogram width is set to zero.

2.3 Cross-comparison of different analysis techniques can highlight problems

TI and perturbation-based analysis techniques have different limitations. Specifically, the accuracy of TI is not a direct function of overlap in energy distributions but instead is a function of the average curvature $\langle \frac{{\partial }^{2}U}{\partial {\mathrm{\lambda }}^2}\rangle$ [30]. On the other hand, perturbation-based techniques do not depend on smoothness of the integrand, but rather on overlap in the sampled energy distributions. Given these differences in input information and limitations, consistency checks across these method families can be a valuable tool for identifying analysis or sampling problems.

In our experience, comparing results from different methods can serve as a warning sign of either insufficient sampling or a λ spacing which is too wide, so it is useful to have a family of analysis methods. Fig. 4 shows hydration free energies for 3-methylindole computed via a variety of methods. As seen from Fig 4(b), the discrepancies between results from the alternate interpolation schemes are most prominent in the vicinities of the rapid change in the $\langle \frac{\partial U}{\partial \mathrm{\lambda }}\rangle$ derivative (the van der Waals lambda states 3–4, 4–5, 5–6, and 7–8). A good practice, thus, is to ensure dense lambda spacing for these regions.

In general, approaches based on Zwanzig’s relation are expected to break down earlier (in terms of phase space overlap) than other approaches, so IEXP and DEXP will tend to become inconsistent even when data may still be sufficient to obtain accurate free energy estimates from many of the other methods (see Fig. 3). In our work, we primarily look for disagreement between TI-based methods and BAR-based methods.

2.4 The free energy change is often broken down into components

Alchemical transformations are usually comprised of several conceptual steps which modify different terms in the potential. For example, solvation or binding free energy calculations are often separated at least into electrostatic and Lennard-Jones components. Thus, the free energy associated with modifying each of these terms in the potential can be computed separately, which can provide some qualitative insight. However, it is important to note that the free energies of each component are path-dependent observables. ΔG_Coul will change in value depending on in what order the electrostatic transformation is carried out, so it cannot be directly considered the electrostatic component of the free energy. In GROMACS, since λ is a vector, multiple steps can be handled within a single set of simulations, and our alchemical-analysis.py tool automatically handles separation of different free energy components. It automatically identifies the number of charging states and prints out along with the total free energy change its breakdown into free energy components, ΔG_Coul and ΔG_vdW. The van der Waals contribution is computed as ΔG_Total −ΔG_Coul, where ΔG_Coul is the free energy change between states differing only in their charge state (in GROMACS, controlled by coul-lambda or fep-lambda). Whenever there are other types of transformations involved (for example, λ controlling restraints holding a ligand in a binding site in a binding free energy calculation [4]) the difference ΔG_Total −ΔG_Coul will not equal ΔG_vdW.³

2.5 Analysis should be carried out on uncorrelated samples

In general, the free energy expressions above are intended to be applied to a set of uncorrelated samples of the relevant observables, but simulation data will be stored more frequently. Thus, we typically need to re-sample the relevant energy derivatives ∂U/∂λ and/or energy differences ΔU_i,j to obtain uncorrelated samples. Particularly, samples used for computing the averages in eqs (1)–(4) are statistically independent samples. There are several techniques to identify and retain independent samples [11, 40, 10] and we find analysis of autocorrelation times to be particularly useful in this regard.

Autocorrelation analysis begins with calculation of the autocorrelation function, which is fairly standard procedure in time series analysis. For a discrete set of N samples occurring time δt apart, the autocorrelation function of the observable A at a given point i is found as

where δA(i) defines the deviation of the current observable from its mean

In our script we use ∂U/∂λ as the observable A, although any of the potential energy differences ΔU_i,j can be used as well.

The autocorrelation time τ is defined as the integral of C_A(i), and becomes noisy at long times, especially at more than half of the simulation time. To obtain reliable estimate of τ, the simulation time should, as a general rule of thumb, exceed 50τ. Once the correlation time is found, a set of independent samples is built up by picking every gth (where g = 1 + 2τ) sample out of the original set. A detailed derivation can be found elsewhere [8].

3.1 Time-reversed convergence plots reveal non-equilibrated regions

The trajectory snapshots to be analyzed must be sampled at equilibrium. To get rid of any non-equilibrated region of the trajectory its location should first be identified. Automatic identification of non-equilibrated regions remains a major research challenge, and in fact equilibration may have occurred prior to start of data collection depending on the equilibration protocol. So assessment of equilibration remains a major task for the researcher conducting free energy simulations, but automated tools can provide at least some help. A standard approach is to look at convergence of the free energy estimate as a function of simulation time, as in Fig. 5.

While convergence analysis is normally applied on data in the order in which it was collected, the same analysis can also be applied to time-reversed data, as in the work of Yang et al. [41]. Yang et al. focused on automatic detection of equilibration based on reverse cumulative averaging starting with the end point of the simulation. Our approach is slightly different, and we simply compare forward and reverse free energy estimates. That is, we might compare an estimate of the free energy change based on the first 10% of the data with an estimate based on the last 10% of the data. This results in two ΔG estimates for each observation time, one using data collected starting from the beginning of the trajectory and moving forwards, and the other using data collected starting from the end and moving backwards (Fig. 5).

Fig. 5 and 6 show convergence of free energy estimates obtained both from normal and time-reversed data. We include the forward estimate as this is the conventional way to examine the data and assess convergence, but we actually find the estimate with time-reversed data more helpful. In a favorable case, both sets of data should converge rapidly to within uncertainty of the final value (Fig. 5(a)). However, consider a hypothetical case where the free energy calculations were begun before the system is at equilibrium, and the first 40% of the data is in fact non-equilibrium. What would we expect in this case?

In the case where a substantial amount of data at the beginning of a set of simulations is not equilibrated, it is reasonable to assume that any free energy estimate based on this data would differ substantially from an estimate based on data from the equilibrated region. Thus as we examine free energy estimates from both normal and time-reversed data, both will reach the same final value but from opposite directions. We would expect that the time-reversed estimate will be steady (within uncertainty) around some value and then leave this value to reach a different final ΔG value as un-equilibrated data begins being included in the free energy estimate, while the forward estimate will exhibit some overall trend as the un-equilibrated data from the beginning of the simulation starts to carry less and less weight. This also means that forward and reverse ΔG estimates will tend to approach the final value from opposite directions-i.e. if the forward estimate ascends to the final value, the reverse estimate will descend to it. Thus the two free energy estimates have at least partial reflection symmetry around the line ΔG = ΔG _final.

To illustrate this, we created a hypothetical case of un-equilibrated data by taking our well-converged simulations of 3-methylindole and adding Gaussian noise to the first 40% of the time series at all lambda values. Fig. 5(b) shows this case where the ΔU_i,j distributions of the first 40% of the time series have been contaminated with the Gaussian noise centered at the original mean of each distribution but with a standard deviation of 5 k_BT units. This is a deliberately artificial example, but it provides a clear demonstration of how this analysis can be useful. This “un-equilibrated” data at the beginning of of the timeseries changes the behavior of the ΔG_reverse(t) function in the expected way – it remains stable essentially within uncertainty of a constant value (20 kcal/mol) over a plateau region extending to around t = 0.6t_total (i.e., until the reverse estimate begins to include the first 40% of the data). After this point, the non-equilibrated data adversely affects the reverse ΔG estimates and both the normal and time-reversed free energy estimates converge to the wrong ΔG value (the correct value is seen in Fig. 5(a)). This graph then suggests a simple solution: We can recover the correct free energy estimate if we simply recognize the region of non-equilibrated data and discard it. Fig. 6 shows the convergence plot obtained via this procedure, with the first 40% of the data removed from the analysis. As expected, the elimination of the non-equlibrium region from the analysis restores the well-behaved character of the ΔG_reverse(t) function, with all data points in both directions essentially lying within uncertainty of the final value.

Our discussion here has assumed that all of the data, after equilibration, is collected at equilibrium from the same distribution so that the forward and reverse free energy estimates agree within error after convergence. But what if there samples are instead collected from multiple metastable states, such as if the system undergoes a conformational transition? One can imagine a situation where the first and last halves of the data were sampled in two different metastable states. This could result in a convergence plot which displays the same behaviors as in the case of unequilibrated data. Thus, our analysis in this subsection refers to the case of adequate sampling, where all relevant metastable states are covered and there is sufficient number of transitions between them within each time block considered in the convergence plot. When this is the case, the deviation of the free energy estimates is entirely due to the presence of unequilibrated data. However, if the populations of different states differ drastically from block to block (i.e. if the transition time between metastable states is a reasonable fraction of the total simulation time) the technique described here will not discriminate between discrepancies caused by unequilibrated data and those caused by an unconverged free energy estimate. In other words, in the limit of adequate sampling, this technique can help to identify unequilibrated data. However, when sampling is not adequate, it will simply indicate a problem which can be due to slow convergence/inadequate sampling of metastable states (slow transitions between states), or to unequilibrated data.

3.2 Practical recommendations for using the convergence plot

In our view, then, the convergence plot is of great value for detecting potentially un-equilibrated data, provided the data to be analyzed was adequately sampled, i.e. with an adequate number of transitions between metastable states and reasonably correct populations in all time blocks. Whenever the time-reversed free energy estimate shows an extended plateau, then after some time begins to exhibit an overall trend in time (as in Fig. 5(b)), it suggests that the new data being included which leads to that trend may be non-equilibrated. Thus this should be considered a substantial warning sign, and data collected prior this point should perhaps be discarded.⁴ In contrast, if both forward and reverse estimates quickly approach the same value as in Fig. 5(a), and there is no plateau in the time-reversed data from which the free energy estimate departs as more data is included, then this measure suggests no concern.

3.3 The overlapping distribution method identifies regions with poor phase-space overlap

In addition to looking for consistency between free energy estimates from several methods discussed in Section 2.3, the Overlapping Distribution Method (ODM) [11,32]—introduced by Bennett [1] under the name of the Curve-Fitting Method— provides another useful technique for spotting trouble. It is a helpful tool for assessing consistency, when combined with another free energy estimator.

We start with the equation 19 from the Bennett paper [1]:

where P_i+1(ΔU) and P_i(ΔU) are the distributions of the potential energy differences between adjacent states obtained when sampling at state λ = i + 1 and λ = i, respectively. The free energy change between the states, ΔG_i,i+1, and corresponding potential energy difference, ΔU, are written in reduced units.

This equation can be rewritten as

where the new distribution functions g_i+1 and g_i were obtained by taking natural logarithm of both sides of Eq. 8 and splitting the ΔU term

where C is an arbitrary constant.

If the difference Δg_i,i+1 = g_i+1 − g_i is plotted versus ΔU, there should be a range of the ΔU values over which Δg_i,i+1 oscillates about the free energy estimate obtained by a standard technique such as BAR, provided there is an overlap between the two distributions and the sampling was sufficient.

This analysis, shown in Fig. 7(a), means that we can graphically inspect the difference Δg_i,i+1 over a range of the ΔU values and it should appear relatively constant. If not, it is a substantial warning sign. All the subplots of Fig. 7(a) exhibit good behavior: the Δg_i,i+1 function oscillates about the ΔG_BAR estimate (shown in purple, with a width corresponding to the uncertainty in ΔG_BAR). In contrast, the Δg_i,i+1 function depicted on the subplot 4–6 of Fig. 8(a) behaves abnormally. Figure 8 was obtained by processing exactly the same data but with the lambda state 5 left out from the analysis. This resulted in the deterioration of the distribution overlap between states 4 and 6, which the ODM helps highlight: the Δg_i,i+1 function is not continuous and does not exhibit saw-like behavior over a substantial interval of the ΔU_i,i–1 values as it does in Fig. 7(a). This tells us that the estimated free energy change for the 4–6 pair of states is getting less reliable, which is a problem.

To further highlight the problem, we exacerbate the overlap by throwing away yet another intermediate state, lambda state 6, from the analysis and re-examine the Δg_i,i+1 function (Fig. 10). This time, there is only one point on the graph for subplot 4–6, and the point is far from the unusually wide (because of the greater uncertainty) ΔG_BAR estimate, indicating clear and severe problems with overlap in this case.

In general, low overlap does not necessarily make the ΔG estimate completely wrong or substantially untrustworthy. However, it does substantially increase the corresponding variance, as only few samples contribute to the free energy estimate. Users themselves must decide what level of uncertainty in the free energy estimate can be tolerated, and whether the poor efficiency identified by this method might suggest restructuring the spacing or number of intermediate states before collecting additional data.

Figure 9 shows the free energies obtained from the analysis with different states omitted. The free energy change estimated between states 4 and 7 through the intermediate states 5 and 6 (−0.028 ± 0.332 kJ/mol) and 6 only (0.601 ± 0.651 kJ/mol) are within statistical noise. However, the 4–7 ΔG estimate without any intermediate states yields 3.986 ± 2.512 kJ/mol which is almost two standard deviations away from the original value of −0.028 ± 0.332 kJ/mol. On the case of extremely poor overlap (as in Fig. 10), the FEP-based methods tend to underestimate the variance, making free energy estimates untrustworthy.

Despite its ability to identify lambda regions with poor and good overlap, the overlapping distribution remains a qualitative method. If one needs to translate “poor overlap” and “good overlap” into concrete numbers, the overlap matrix, discussed in the next subsection, should be employed.

3.4 The overlap matrix is a quantitative estimator of the phase-space overlap

The overlap matrix is a helpful tool in finding the magnitude of the phase space overlap and we recommend using it as a consistency check whenever the free energy estimate relies on the overlap, as in the case of the FEP-based methods.

If we define the weight of each of the N samples x_n (collected from all K states) in the ith lambda state as:

then the covariance matrix Θ of the vector of reduced free energies βG can be written (From Eq. D6 of Ref. [35]) as

where ⁺ is the Moore-Penrose pseudoinverse, N is a diagonal matrix with its elements N_ii the number of samples collected in state i, and O is the overlap matrix, defined as:

We note that $e^{\beta G_i-\beta U_i(x_n)}$ in eq (11) is the probability p_i(x_n) of sample x_n occuring when simulating state i. The unnormalized probability of sample x_n is the Boltzmann weight $e^{-\beta U_i(x_n)}$ and since $e^{-\beta G_i}={\displaystyle \int e^{-\beta U_i(\mathrm{x})}}d\mathrm{x},e^{-\beta G_i}$ is the normalizing constant.

The overlap matrix is a K × K matrix with entries:

The N samples were collected with N₁ samples from the p₁(x) distribution, N₂ samples from the p₂(x) distribution, and so forth. This combination of the K distributions is known as a mixture distribution, and can be written mathematically as:

where $N_{\mathit{tot}}={\displaystyle {\sum }_{k=1}^K{N}_k}$.

We next note that the Monte Carlo estimate of the integral ∫ A(x)p(x)dx, where A is some function of x, is $\frac{1}{N}{\displaystyle {\sum }_{n=1}^{N}A(x_n)}$, where the samples x_n are drawn from the probability distribution p(x). This means that each element of the overlap matrix can be seen as a Monte Carlo estimate of the integral:

where the averages are over the probability distribution of samples in state j. O_ij can therefore be interpreted as the average probability of a sample generated in state j being observed in the ith state. This average is computed over samples collected from all of the states, not just the samples from state j. We can easily see that we must have Σ_i O_ij= 1. Because O_ij is a stochastic (or Markov) matrix, its eigenvalues are all real and positive, and the largest is 1. In fact, since it is a diagonal matrix (N) times a symmetric matrix (W^TW), then all the eigenvalues are real and positive.

We can then write (using the standard notation for the eigenvalue decomposition):

In this case, because it is a diagonal matrix, we can give a simple formula for the pseudoinverse; (Λ⁻¹ –I)⁺ is a diagonal matrix with one zero diagonal entry (corresponding to the largest eigenvalue, which is 1), and the other entries corresponding to the ith eigenvalue λ_i being λ_i/(1 – λ_i). The uncertainty in any free energy difference between states i and j is Θ_ij = Θ_ii + Θ_jj − 2Θ_ij which makes it difficult to write explicit formulas of the variance in terms of the overlap matrix. In general, uncertainties in free energy differences will be larger when the eigenvalues λ_i/(1 – λ_i) are large, which will occur when the eigenvalues are close to 1. Because of the factor of N⁻¹, any variance can be made arbitrarily low with enough samples. Clearly, the most efficient choices of λ points to simulate will be ones leading to smaller eigenvalues in the overlap matrix.

When will these eigenvalues be close to one? Again, it is difficult to completely generalize, but for stochastic matrices, i.e. matrices with the form of O with rows summing to 1, the smallest K – 1 eigenvalues approach one when the matrices can be nearly decomposed into independent block matrices. The absence of eigenvalues close to one indicates that the matrix is more connected.

3.5 Practical recommendations for using the overlap matrix

We cross-checked the overlap matrix with the overlapping distribution method (Fig. 7, 8, and 10) and arrived at several recommendations for trustworthy free energy calculations. First, trustworthy results should in general have at least a tridiagonal overlap matrix-that is, all pairs of adjacent states should have substantial overlap, and no element should be zero in the main diagonal and the first diagonals above and below the main one. Second, these nonzero elements should be appreciably different from zero. In our experience, with enough samples, values as low as 0.03 (Fig. 8) seems to be tolerable to yield a reliable free energy estimate (as long as the resulting error is sufficiently low), though obviously the number of samples required for an accurate estimate increases with decreasing overlap. Anything below that number should serve as a caution sign as the FEP-based methods tend to underestimate the variance when the phase space overlap is that low; in this case, not only will the estimated variance increase but the free energy estimate itself will likely be substantially incorrect, perhaps by far more than the estimated variance.

4.1 Tool execution and usage

Our tool is executed by the command python alchemical-analysis.py [options] with the list of options provided below. In general, plots and options are as described above, except when they relate to GROMACS in particular, in which case additional information is provided below.

Options include:

• –

  Simulation temperature

• –

  Directory with the input data

• –

  Datafile prefix and suffix

• –

  Time prior to which the data is to be discarded

• –

  Names of the free energy estimators to be used

• –

  The units the free energies are to be reported in

• –

  Number of decimal places the free energies are to be reported with

• –

  Graphical functionality (discussed above):

  • –

    The ΔG_i,i+1 vs. λ_i,i+1 bar plot as in Fig. 4(a)

  • –

    The $\langle \frac{\partial U}{\partial \mathrm{\lambda }}\rangle$ vs. λ plot as in Fig. 4(b)

  • –

    The ΔG vs. time plot, as in Fig. 5

  • –

    The overlap matrix calculation, as in Fig. 7

  • –

    The consistency check based on the overlap distribution method, as in Fig. 10

The script outputs a text file, results.txt, with a table of free energy differences computed by means of various methods for each pair of adjacent states, as well as overall totals. This file also shows the ΔG breakdown into its components discussed in Section 2.4. Full precision data is dumped to a Python pickle file, results.pickle, where it is stored in a form a class with multiple instances whose names are self-explanatory.