# Analyzing Simulation Results

Free Energy Fundamentals |
---|

Methods of Free Energy Simulations |

Free Energy How-to's |
---|

Once you have setup your transformation pathway, defined necessary intermediate states, run your equilibrium simulation, and acquired the correct data from the runs, you may now analyze the data to get free energies. Each of the free energy techniques discussed in the theory section of these fundamentals talks about how you get the free energies, but will also be recapped here. This page also talks about how you can get uncertainty estimates with the bootstrap method.

## Contents

# Calculating Free Energies

Shown below is a brief summary of the information needed and the method to calculating free energy with the various techniques. If you have read the pages for the free energy techniques, this will be a recap.

- EXP is a very straightforward calculation. You will need the [math]\Delta U[/math] from a given pair of states. From there, you can calculate the free energy from the effective sample size found by correcting for correlation/subsampling. Variances will be additive.
- TI needs [math]\frac{dU}{d\lambda} [/math] and the average at each of [math]K[/math] states needs calculated. Since its more common to have discrete [math]\lambda[/math] states, you will need to choose an appropriate weighting method to correctly calculate the free energy from the uncorrelated/subsampled data. By-state Variance will not be additive, but this is simple to account for.
- WHAM requires binning all your results then calculating [math]\Delta U[/math] from all states. It is highly recommended that you take advantage of the tools already out there for calculating WHAM instead of writing your own. There is no direct way to calculate variance, and so methods such as bootstrap sampling are needed.
- BAR requires the [math]\Delta U[/math] from two states. An iterative, often numeric, solution is needed to find the free energy, but this is relatively straightforward and can be done for each pair of states. The variance does not have as simple a relationship as with TI, and methods such as bootstrap sampling are recommended. There is a Python implementation of BAR available with PyMBAR.
- MBAR also requires [math]\Delta U[/math], but it needs it for all combination of states. Again, an iterative set of solutions is needed and it is not recommended that users attempt to code MBAR themselves. Instead, please take advantage of a Python implementation available called PyMBAR. Error estimates are directly available with MBAR, however, it can still benefit from methods such as bootstrap sampling.

# Bootstrap Sampling

Bootstrap sampling is a statistical tool in when we can get decent estimates for the uncertainty with very limited data;^{[1]} this method works with all techniques presented in the fundamentals section. We carry out bootstrap sampling as the time taken to do bootstrap is often substantially less than the time it would take to generate new samples.

If we assume we have some function, [math]F(x)[/math] that is computed from [math]N[/math] data samples, [math]x_1,x_2\ldots,x_N[/math], and we have that [math]\lim_{N\to\infty} F(x) = \overline{F}[/math] where [math]\overline{F}[/math] is a constant. As an example, say that [math]F(x)[/math] is the simple average of all [math]x[/math]

- [math]\displaystyle F(x)=\frac{1}{N}\sum\limits_{i=1}^N x_i[/math],

or the average of some function [math] X(x)[/math] such that

- [math]\displaystyle F(x)=\frac{1}{N}\sum\limits_{i=1}^N X(x_i)[/math];

we are not limited to such simple choices as the functions could be as complicated the MBAR or WHAM free energy estimators. To calculate the bootstrap variance, carry out the following procedure.

**Pick with replacement**[math]n[/math] samples from the complete list of samples [math]\vec{x}\{x_1,\ldots,x_N\}[/math] to create a new set of samples [math]\vec{x}_i[/math]. Since you are picking randomly with replacement, there are bound to be repeated samples. This method is called sampling from the*empirical distribution,*that is, sampling from the distribution we measured, rather than true distribution. For example, if [math]\vec{x} = \{1,2,6,4,3\}[/math], a possible set of [math]\vec{x}_i[/math] might be [math]\vec{x}_i = \{6,1,1,4,4\}[/math]; note that [math]\vec{x}_i[/math] is*the same size*as the original [math]\vec{x}[/math] set.**Compute [math]\hat{F}_i = F(x_i)[/math].**That is, find the estimate for your function based on the bootstrap sample taken from the empirical distribution. If we were finding simple averages, we would calculate the average of the bootstrap sample. If we were calculating free energies with MBAR, we would generate a bootstrap sample at*each*[math]K[/math] state and estimate the full free energy at each state based on the bootstrap samples.**Repeat steps 1 and 2**for [math]M[/math] number of times. You will need at least 50-200 times to obtain accurate variances^{[1]}. If the calculation in step 2 is cheap, then [math]M[/math] can be very large to get even better estimates of the variance; relative uncertainty scales as [math]M^{-1/2}[/math] and can take more than 1000 steps to get variance within 1% convergence.**Compute the variance**[math]\mathrm{var}\left(\hat{F}\right)[/math], which is the variance over the [math]M[/math] bootstrap values. The standard deviation is then just the square root of the bootstrap variance.

The values of [math]\hat{F_i}[/math] are the bootstrap distribution of the function [math]F(x)[/math]. In the large sample limit, the bootstrap distribution will converge to the true distribution of [math]\overline{F}[/math] under most conditions, but will still be a good estimate for moderate values of [math]N[/math]. If N is too low though, the bootstrap distribution will not have good agreement with the true distribution.

Bootstrap sampling can be used with *any* statistical estimator, regardless of how complicated. There is the extra overhead of calculating [math]F(x)[/math] [math]M[/math] times, but this is often negligible compared to how long it would take to collect [math]M[/math] times the number of data. For MBAR, the time required for bootstrap is only 5-10min, where as TI is only seconds. There is also a variant of bootstrap called "moving block bootstrapping" which accounts for Simulation Information Gathering#Correlation time correlations *without* subsampling by taking random block of length [math]\tau[/math].^{[1]}

## Example of Bootstrap Method

Lets take a simple example where we will sample values of [math]x[/math] from a 1D Gaussian with mean 0 and and variance 1 ([math]x[/math]~[math]N(0,1)[/math]). We will now compute EXP with [math]k_BT=1[/math] for this sampling by

- [math]F(x_i) = -\ln \left[N^{-1}\sum_{i=1}^N \exp(-x_i)\right][/math].

We then compute the distribution of free energies obtained by [math]N[/math] points for 100,000 bootstrap samples, and 100,000 independent samples (see figure on right). For small [math]N[/math], the distributions are not close, but nearly converge for larger [math]N[/math]. Even at N=200, the variances are near identical. To emphasize the point, the "Repeated Samples" distribution uses data from [math]1000\times 100,000[/math] points collected, where the "Bootstrap Samples" needs only *one* set of 1000 samples run through the bootstrap process 100,000 times. The same amount of analysis is carried out in both cases, but the bootstrap samples only require 0.001% of the data collection. One catch to the bootstrap method that should be noted is that rare events requiring more than 1000 samples would not have been observed, so it still takes careful analysis to properly use bootstrap sampling.

# References

- ↑
^{1.0}^{1.1}^{1.2}Efron, B., and Tibshirani, R. J.*An Introduction to the Bootstrap*; Chapman and Hall/CRC:Boca Raton, FL, 1993. - Find at Cite-U-Like