\section{MBAR in alchemical free energy calculations}

In this test example, we compare MBAR with other standard methods of
computing free energies for alchemical transformations.  The other
methods for this comparison are thermodynamic integration (TI), the
Zwanzig formula in both the insertion and deletion directions (denoted
IEXP and DEXP respectively), and a pairwise application of Bennett
acceptance ratio method (BAR) between consecutive steps.

We test these methods on the problem of calculating the hydration free
energy of a solute by decoupling the interactions between the solute
and the solvent.~\cite{???}.  In the $\lambda=0$ state, both the
intermolecular Coulombic and Lennard-Jones energies are turned on; in
the $\lambda=1$ state, the intermolecular solute-solvent Columbic and
Lennard-Jones energies are turned off, but the intramolecular
solute-solvent interactions remain on.  The free energy of this
transformation is thus equivalent to the free energy of transfer from
the vapor state to the hydrated state, or chemical potential at
infinite dilution in water.

The solute molecule used in this test is 3-methylindole, the side
chain analog of tryptophan that results from truncating the tryptophan
side chain at that $\alpha$-carbon and replacing it with a hydrogen.
The molecule is both large and carries moderately large partial
charges, resulting in a challenging but biophysically relevant test
case.

Force field parameters for 3-methylindole were taken from OPLS-AA, and
were identical to the parameters used in~\cite{????}.  The solute were
simulated in cubic box with fixed size length 3.0073 nm, with 891
TIP3P water molecules.  All simulations were done with a preliminary
version of Gromacs 4.5, modified to give slightly different output for
free energy calculations for ease of analysis.

We simulated the solvation with molecular dynamics using the velocity
Verlet (first included in Gromacs 4.5) with modifications to the
Gromacs free energies targeted for inclusion in Gromacs 4.6.  A time
step of 2 ns was used, with linear center of mass removal every step.
Neighbor lists were recomputed every 10 steps, with an neighbor-list
radius of 1 nm.  Particle mesh Ewald electrostatics (PME) were used to
treat long range electrostatic interactions, with 6th order
interpolation, with a real space cutoff of 0.9 nm, an Ewald tolerance
of $10^{-6}$ with the FFT grid spacing at 0.1 nm.  The real space part
was tapered from 0.88 to 0.9 nm to improve energy conservation.
Lennard-Jones interactions were cutoff at 0.9 nm, with a tapered
switch from 0.8 to 0.9 nm.  Bond to hydrogens were constrained with
RATTLE~\cite{???} or settle in the case of the TIP3P
water.~\cite{???}. Initial velocities were chosen from a
Maxwell-Boltzmann distribution at 298 K.  Volume was fixed, and
temperature was controlled using Nose-Hoover chains with length 10(?)
with a time constant of 10 ps.  A linearly dependent $\lambda$
dependent analytical dispersion correction was added at each step.

We ran 38 4 ns simulations, each at a different $\lambda$ state in the
decoupling trajectory.  Values of $\dhdl$ and potential energy
differences between $\lambda$ states were printed every 200 steps (0.4
ps).  First, intermolecular Coulombic interactions were decoupled off
linearly.  Next, Lennard-Jones intermolecular interactions were turned
off with a ``soft core'' interaction as described in~\cite{????}  We
also ran 11 50 ns simulations of 3-methylindole with the same
decoupling pathway.  In this case values of $\dhdl$ and potential
energy differences between $\lambda$ states were printed every 500
steps (1 ps).  In both cases, the first 100 samples were discarded
(approximately 40 ps for the 38 step set, 100 ps for the 11 step set),
for a total of 150.5 ns for the 11 state set, and 549 ns for the
larger set, so 3.65 times as much simulation time is performed in the
11 ns step set, although since the data is collected less frequently
(0.5 ps vs 0.2 ps), the actual amount of data collected is only about
1.5 times greater.

For the 11 step set, the four Coulomb lambda values were 0.0, 0.2,
0.5, 1.0, and the eight Lennard-Jones lambda values were 0, 0.3, 0.5
0.65, 0.75, 0.85, 0.9, and 1.0 (the last Coulomb and the first
Lennard-Jones states are shared). For the 38 steps set, the 11 Coulomb
lambda values were 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
1.0, and the 28 Lennard-Jones lambda values were 0.0, 0.1, 0.2, 0.3,
0.35, 0.4, 0.45, 0.5, 0.55, 0.58, 0.6, 0.62, 0.64, 0.66, 0.68, 0.70,
0.72, 0.74, 0.76, 0.78, 0.80, 0.82, 0.84, 0.86, 0.88, 0.92, 0.96, and
1.0.

The data was subsampled by computing the autocorrelation time of the
total $\dhdl$, and identifying the subsampling time as the first time
the autocorrelation function crosses zero.  Other possibilities, such
as fitting to an exponential and subsampling at intervals of $2\tau$
are possible~\cite{???}In most cases, the correlation times were
between 1 and 2 samples for both sets, resulting in the majority of
the data being retained after subsampling.

Uncertainty estimates were calculated using both analytical error
estimates and a bootstrap error estimate.  For TI, the standard
deviation in the average of $\dhdl$ at each $\lambda$ value was
calculated, and this uncertainty was propagated through the numerical
integration, with the assumption of statistically uncorrelated error.
For DEXP, IEXP, and BAR, the uncertainty in the free energy between
each consecutive pair of states was calculated using the standard
analytical uncertainty estimates of each method, and then propagated
into the sum by again assuming the individual free energies are
uncorrelated.  We note that this is not necessarily a good
approximation; for BAR, for example, the free energy estimate between
states 1 and 2 depends on data collected from simulations at state 1
and 2, and the free energy estimate between states 2 and 3 depends on
data collected from simulations of state 2 and 3.  Therefore, the two
free energy differences both depend on the same set of simulation data
of state 2, and will not be entirely statistically independent.  For
MBAR, the uncertainty estimate between the initial and final state as
described in ~\cite{????} was used, which makes makes no assumptions
about correlation between states.

We also calculated the overall uncertainties for all five methods
using the bootstrap method.~\cite{????} Simulation data from each
state of each data set was resampled with replacement, and all free
energies were resampled from the collected data were recalculated,
with this resampling process being repeated 100 times.  The bootstrap
uncertainty of each method is calculated by taking the standard error
of the distribution of 100 free energies obtained by from the 100
resampled data sets.  With sufficient resamplings, this process gives
the variance in the empirical distribution of samples collected, which
is equal to the full distribution of data in the limit of large
samples of data collected.  Trials indicate that with 100 samples from
the bootstrap distribution, the bootstrap variance is converged to
within approximately 10\%, in line with general rules of thumb to use
between 50 and 200 bootstrap samples~\cite{????}

Comparisons of the 11 state and 38 state state set is shown in
Table~\ref{table.set11} and Table~\ref{table.set38} respectively, with
both the free energies between individual states and totals shown,
along with uncertainties computed from the analytical error estimates.
Inspection of the 38 state set indicates good agreement between BAR,
MBAR, DEXP, and IEXP for the 38 state test set, which in turn agree
with the BAR and MBAR estimates for the 11 state test set.  Based on
the agreement of BAR and MBAR with the 11 step and 38 step sets, as
well as the agreement with IEXP and DEXP with the 38 step set, we can
conclude that the free energies are 18.95 kJ/mol$\pm ~\approx$ 0.1
kJ/mol, where the uncertainty depends on the method and the data set.
However, for 11 states, TI significantly underestimates the true free
energy.  Even with 38 states, the free energy calculated with TI
appears to be off the consensus value by approximately 0.55 kJ/mol.
Inspection of the individual states indicates that most of this bias
is in the first few states of the Lennard-Jones decoupling; for
essentially all other states, TI and the other methods are within
statistical uncertainty.  With a different selection of states, it
seems likely that 38 states would be sufficient for TI; however, MBAR
and BAR seem much more robust different selection of values for
$\lambda$.  One might typically choose TI states based on lowering the
variance.  Although the variance in $\dhdl$ is largest in the
intervals between the 1st and 2nd and 2nd and 3rd Lennard-Jones
$\lambda$, it is not significantly larger than the variance from other
states.  For MBAR and BAR, we can switch from 11 states to 38 states
and get exactly the same answer to within statistical uncertainty.
Improved numerical integration methods, such as Gaussian quadrature
might in some cases lead to more efficient TI estimates.  especially
for the Coulombic portion which is more smooth. However, this
requires knowing the curvature of $\dhdl$ before the simulation which
may not always be possible.

Additionally, MBAR and BAR do not require any analytical derivatives
being computed, only reevaluation of the potential energies at
additional $\lambda$ states, and then only at a small fraction of the
overall states simulated.  This can be done very efficiently either
during the simulation or using postprocesing. 

To directly test the utility of analytical and bootstrap estimates of
uncertainty, we divided the each of the 11 state 50 ns trajectories
into 40 equal sub-trajectories, and individually computed the free
energies of solvation for each complete set of sub-trajectories.  We
then computed a direct estimate of the uncertainty by computing the
standard deviation of this population of 40 shorter sets of
simulations (see Table~\ref{table.replicate40}).  For each of these 40
set, we also computed bootstrap errors by resampling each of the 40
simulation sets, as well as computing the analytical error estimate.  

Using Table~\ref{table.replicate40}, we can examine the bias in each
method by comparing the average free energies over all 40 sets with
the free energies obtained from the full set.  We see that DEXP and
IEXP are have significant bias.  Deletion exponential averaging (DEXP)
has significant bias, with a significant difference between the full
data set and 40 subset values.  Insertion exponential averaging (IEXP)
is significantly improved over deletion exponential averaging,
corresponding, and has relatively small bias compared to DEXP. It is
well established that going from regions of high entropy to regions of
low entropy yields lower statistical variance ~\cite{Kofke}, which is
reflected in the superiority of IEXP to DEXP.  However, there is still
moderate bias (0.3 kJ/mol).  BAR and MBAR exhibit very little bias
with the averages over 40 subsets being almost equal to the free
energy calculated from the full data set. TI results are the same,
indicating that the averaging of $\dhdl$ values has no bias, but the
bias from the numerical integration is obviously quite significant.

We then also compare the analytical estimate and direct uncertainties
for the four methods.  In all cases, the direct uncertainties are
greater than the predicted or bootstrap uncertainties by about 20\%.
Since for directly generated data from harmonic oscillators, the
direct uncertainties for BAR and MBAR match very well with the
predicted uncertainties, this may indicate that here are underlying
correlations in the data that are not taken into account by the
subsampling performed in this example.  However, it does mean that the
uncertainties are roughly in line with estimates. 

The BAR analytical estimate underestimates the direct uncertainty to a
greater extent than MBAR does, as does the analytical error estimate
for TI, while both IEXP and DEXP have significantly higher
uncertainties than MBAR and BAR.  BAR, MBAR, and TI all have similar
direct uncertainty, while IEXP and DEXP have successively greater
uncertainties.

Interestingly, analytically estimated errors match bootstrap errors
very precisely for MBAR, BAR, IEXP, and DEXP, while TI does not.
Analytical error estimate for MBAR and BAR are extremely robust,
yielding virtually the same answer despite resampling, which is a
useful property of the method. 

%BAR applied pairwise between neighboring states and MBAR give
%essentially identical average values (19.004 kJ/mol vs. 18.986 kJ/mol)
%indicating very little bias, as well almost identical directly
%calculated uncertainties over the set of 40 free energies (0.531 vs
%0.523 kJ/mol). However, the predicted uncertainties of BAR
%underestimate the true uncertainties by a larger amount than MBAR
%(0.353 kJ/mol predicted and 0.531 kJ/mol actual for BAR, versus 0.401
%predicted vs 0.523 actual for MBAR).  Additionally, average bootstrap
%uncertainties are different from both the predicted and the actual
%uncertainties for BAR, whereas the predicted and bootstrap uncertainty
%agree well for MBAR.

\begin{table}\renewcommand{\arraystretch}{0.75}
\caption{Free energy of 3-methylindole with 38 total states, compared with 5 methods, including analytical error estimates}
\begin{tabular}{lccccc}
         &    TI (kJ/mol)      &   DEXP (kJ/mol)       &    IEXP (kJ/mol)     &   BAR (kJ/mol)        &     MBAR (kJ/mol) \\
    0:  &  13.122 $\pm$  0.010 &  13.078 $\pm$  0.020  & 13.095 $\pm$  0.024  &  13.091 $\pm$  0.011  &  13.090 $\pm$  0.010  \\
    1:  &  13.804 $\pm$  0.015 &  13.587 $\pm$  0.034  & 13.540 $\pm$  0.054  &  13.586 $\pm$  0.015  &  13.584 $\pm$  0.014  \\
    2:  &  12.419 $\pm$  0.022 &  11.966 $\pm$  0.166  & 12.044 $\pm$  0.141  &  12.049 $\pm$  0.023  &  12.044 $\pm$  0.020  \\
    3:  &  -9.937 $\pm$  0.027 &  -4.486 $\pm$  0.419  & -5.179 $\pm$  0.051  &  -5.165 $\pm$  0.026  &  -5.175 $\pm$  0.023  \\
    4:  &   2.845 $\pm$  0.018 &   3.498 $\pm$  0.202  &  3.365 $\pm$  0.031  &   3.382 $\pm$  0.019  &   3.382 $\pm$  0.018  \\
    5:  &  -0.188 $\pm$  0.020 &   1.030 $\pm$  0.173  &  0.395 $\pm$  0.039  &   0.372 $\pm$  0.021  &   0.371 $\pm$  0.021  \\
    6:  &  -4.652 $\pm$  0.018 &  -4.227 $\pm$  0.199  & -4.366 $\pm$  0.044  &  -4.368 $\pm$  0.021  &  -4.368 $\pm$  0.021  \\
    7:  &  -7.696 $\pm$  0.016 &  -8.741 $\pm$  0.050  & -8.757 $\pm$  0.054  &  -8.748 $\pm$  0.018  &  -8.748 $\pm$  0.016  \\
    8:  &  -3.158 $\pm$  0.004 &  -3.175 $\pm$  0.006  & -3.176 $\pm$  0.007  &  -3.177 $\pm$  0.004  &  -3.176 $\pm$  0.003  \\
    9:  &  -2.114 $\pm$  0.005 &  -2.013 $\pm$  0.007  & -2.013 $\pm$  0.008  &  -2.014 $\pm$  0.005  &  -2.014 $\pm$  0.004  \\
\hline  
TOTAL:  &  14.446 $\pm$  0.053 &  20.516 $\pm$  0.563  & 18.947 $\pm$  0.182  &  19.008 $\pm$  0.056  &  18.990 $\pm$  0.063  \\
\hline  
\end{tabular}
\label{table.set11}
\end{table}


\begin{table}\renewcommand{\arraystretch}{0.75}
\caption{Free energy of 3-methylindole with 38 total states, compared
  with 5 methods, including analytical error estimates.}
\begin{tabular}{lccccc}
       &      TI (kJ/mol)      &     DEXP (kJ/mol)     &    IEXP (kJ/mol)     &     BAR (kJ/mol)     &    MBAR (kJ/mol)     \\ 
    0: &    6.978 $\pm$  0.013 &   6.985 $\pm$  0.021  &  6.969 $\pm$  0.022  &  6.975 $\pm$  0.014  &  6.982 $\pm$  0.012  \\
    1: &    6.077 $\pm$  0.014 &   6.067 $\pm$  0.020  &  6.088 $\pm$  0.022  &  6.070 $\pm$  0.014  &  6.080 $\pm$  0.010  \\
    2: &    5.268 $\pm$  0.014 &   5.238 $\pm$  0.020  &  5.265 $\pm$  0.024  &  5.258 $\pm$  0.014  &  5.246 $\pm$  0.009  \\
    3: &    4.520 $\pm$  0.014 &   4.539 $\pm$  0.021  &  4.489 $\pm$  0.021  &  4.510 $\pm$  0.014  &  4.504 $\pm$  0.008  \\
    4: &    3.864 $\pm$  0.012 &   3.853 $\pm$  0.017  &  3.875 $\pm$  0.018  &  3.857 $\pm$  0.012  &  3.860 $\pm$  0.007  \\
    5: &    3.313 $\pm$  0.011 &   3.306 $\pm$  0.015  &  3.302 $\pm$  0.018  &  3.307 $\pm$  0.011  &  3.305 $\pm$  0.006  \\
    6: &    2.827 $\pm$  0.011 &   2.820 $\pm$  0.016  &  2.827 $\pm$  0.017  &  2.824 $\pm$  0.011  &  2.819 $\pm$  0.006  \\
    7: &    2.393 $\pm$  0.011 &   2.392 $\pm$  0.015  &  2.382 $\pm$  0.016  &  2.390 $\pm$  0.011  &  2.384 $\pm$  0.005  \\
    8: &    1.980 $\pm$  0.011 &   1.994 $\pm$  0.015  &  1.975 $\pm$  0.017  &  1.979 $\pm$  0.011  &  1.984 $\pm$  0.005  \\
    9: &    1.583 $\pm$  0.012 &   1.586 $\pm$  0.016  &  1.591 $\pm$  0.019  &  1.584 $\pm$  0.012  &  1.601 $\pm$  0.006  \\
   10: &   -5.332 $\pm$  0.025 &  -4.986 $\pm$  0.054  & -5.038 $\pm$  0.025  & -5.026 $\pm$  0.020  & -5.027 $\pm$  0.011  \\
   11: &   -1.205 $\pm$  0.016 &  -0.977 $\pm$  0.038  & -1.050 $\pm$  0.021  & -1.044 $\pm$  0.016  & -1.050 $\pm$  0.010  \\
   12: &    0.787 $\pm$  0.017 &   0.863 $\pm$  0.052  &  0.881 $\pm$  0.023  &  0.882 $\pm$  0.017  &  0.896 $\pm$  0.011  \\
   13: &    0.780 $\pm$  0.009 &   0.778 $\pm$  0.018  &  0.789 $\pm$  0.012  &  0.790 $\pm$  0.009  &  0.804 $\pm$  0.006  \\
   14: &    0.889 $\pm$  0.010 &   0.878 $\pm$  0.020  &  0.912 $\pm$  0.013  &  0.900 $\pm$  0.010  &  0.898 $\pm$  0.006  \\
   15: &    0.888 $\pm$  0.011 &   0.878 $\pm$  0.028  &  0.888 $\pm$  0.014  &  0.896 $\pm$  0.011  &  0.894 $\pm$  0.006  \\
   16: &    0.780 $\pm$  0.012 &   0.778 $\pm$  0.025  &  0.803 $\pm$  0.016  &  0.791 $\pm$  0.012  &  0.790 $\pm$  0.007  \\
   17: &    0.595 $\pm$  0.014 &   0.578 $\pm$  0.027  &  0.606 $\pm$  0.019  &  0.604 $\pm$  0.014  &  0.572 $\pm$  0.008  \\
   18: &    0.198 $\pm$  0.009 &   0.235 $\pm$  0.015  &  0.186 $\pm$  0.012  &  0.200 $\pm$  0.009  &  0.178 $\pm$  0.005  \\
   19: &    0.026 $\pm$  0.007 &   0.030 $\pm$  0.011  &  0.022 $\pm$  0.009  &  0.026 $\pm$  0.007  &  0.024 $\pm$  0.003  \\
   20: &   -0.073 $\pm$  0.007 &  -0.066 $\pm$  0.011  & -0.075 $\pm$  0.009  & -0.072 $\pm$  0.007  & -0.071 $\pm$  0.004  \\
   21: &   -0.195 $\pm$  0.007 &  -0.182 $\pm$  0.011  & -0.205 $\pm$  0.010  & -0.193 $\pm$  0.007  & -0.187 $\pm$  0.004  \\
   22: &   -0.344 $\pm$  0.008 &  -0.334 $\pm$  0.013  & -0.336 $\pm$  0.011  & -0.341 $\pm$  0.008  & -0.330 $\pm$  0.004  \\
   23: &   -0.499 $\pm$  0.009 &  -0.522 $\pm$  0.015  & -0.482 $\pm$  0.012  & -0.495 $\pm$  0.009  & -0.507 $\pm$  0.005  \\
   24: &   -0.715 $\pm$  0.009 &  -0.680 $\pm$  0.015  & -0.734 $\pm$  0.013  & -0.713 $\pm$  0.009  & -0.724 $\pm$  0.005  \\
   25: &   -1.006 $\pm$  0.010 &  -0.987 $\pm$  0.016  & -1.012 $\pm$  0.013  & -1.004 $\pm$  0.010  & -0.982 $\pm$  0.006  \\
   26: &   -1.286 $\pm$  0.010 &  -1.310 $\pm$  0.017  & -1.271 $\pm$  0.014  & -1.287 $\pm$  0.010  & -1.270 $\pm$  0.006  \\
   27: &   -1.539 $\pm$  0.010 &  -1.544 $\pm$  0.016  & -1.531 $\pm$  0.015  & -1.544 $\pm$  0.010  & -1.550 $\pm$  0.006  \\
   28: &   -1.745 $\pm$  0.009 &  -1.761 $\pm$  0.014  & -1.760 $\pm$  0.013  & -1.755 $\pm$  0.009  & -1.762 $\pm$  0.005  \\
   29: &   -1.860 $\pm$  0.008 &  -1.856 $\pm$  0.012  & -1.885 $\pm$  0.011  & -1.872 $\pm$  0.008  & -1.863 $\pm$  0.004  \\
   30: &   -1.846 $\pm$  0.007 &  -1.861 $\pm$  0.009  & -1.841 $\pm$  0.010  & -1.853 $\pm$  0.007  & -1.844 $\pm$  0.004  \\
   31: &   -1.714 $\pm$  0.006 &  -1.732 $\pm$  0.008  & -1.706 $\pm$  0.009  & -1.719 $\pm$  0.006  & -1.728 $\pm$  0.003  \\
   32: &   -1.535 $\pm$  0.005 &  -1.533 $\pm$  0.007  & -1.547 $\pm$  0.007  & -1.539 $\pm$  0.005  & -1.548 $\pm$  0.002  \\
   33: &   -1.329 $\pm$  0.004 &  -1.330 $\pm$  0.005  & -1.330 $\pm$  0.006  & -1.331 $\pm$  0.004  & -1.332 $\pm$  0.002  \\
   34: &   -1.951 $\pm$  0.006 &  -1.955 $\pm$  0.009  & -1.950 $\pm$  0.009  & -1.951 $\pm$  0.006  & -1.956 $\pm$  0.003  \\
   35: &   -1.022 $\pm$  0.005 &  -1.007 $\pm$  0.006  & -1.027 $\pm$  0.007  & -1.016 $\pm$  0.004  & -1.015 $\pm$  0.003  \\
   36: &   -0.154 $\pm$  0.003 &  -0.149 $\pm$  0.005  & -0.141 $\pm$  0.005  & -0.146 $\pm$  0.003  & -0.141 $\pm$  0.002  \\
\hline
TOTAL: &   18.398 $\pm$  0.067 &  19.025 $\pm$  0.126  & 18.927 $\pm$  0.093  &  18.942 $\pm$  0.064 &  18.938 $\pm$  0.080   \\
\hline
\end{tabular}
\label{table.set38} 
\end{table}

% loading time for 11 states = 74 s
% loading time for 38 states = 122 s
\begin{table}\renewcommand{\arraystretch}{0.75}
\begin{tabular}{lccccc}
\hline
               &   TI              &    DEXP             &     IEXP             &    BAR           & MBAR               \\ 
\hline
                         & \multicolumn{4}{c}{11 states} &                                                                  \\
\hline
Value          &  14.446           &   20.515            &   18.946             &  19.008          & 18.990             \\
Predicted Unc. &   0.053           &    0.563            &    0.182             &   0.056          &  0.063             \\
Bootstrap Unc. &   0.073           &    0.544            &    0.156             &   0.070          &  0.067             \\
\hline
                         & \multicolumn{4}{c}{38 states} &                                                                  \\
\hline
Value          &  18.398           &    19.025           &   18.926             &  18.942         &  18.938             \\
Predicted Unc. &   0.067           &     0.126           &    0.093             &   0.064         &   0.080             \\
Bootstrap Unc. &   0.077           &     0.115           &    0.077             &   0.076         &   0.065             \\
\hline
                  & \multicolumn{4}{c}{40 fold replicate of 11 step system} &                                           \\
\hline
Av. Value      &  14.447           &  21.687            &   18.771             & 19.004           & 18.986              \\   
Direct    Unc. &   0.538           &   2.272            &    1.038             &  0.531           &  0.532              \\
Predicted Unc. &   0.338$\pm$0.003 &   1.512$\pm$0.483 &    0.875$\pm$ 0.204  &  0.353$\pm$0.002 &  0.401$\pm$0.002     \\
Bootstrap Unc. &   0.440$\pm$0.032 &   1.551$\pm$0.605 &    0.873$\pm$ 0.229  &  0.422$\pm$0.035 &  0.400$\pm$0.032     \\
\hline
\end{tabular}
\label{table.replicate40}

\end{table}

Overall, MBAR appears to be the best method for computing free
energies of larger, biophysically relevant solutes.  It is
significantly more statistically efficient and less bias-prone than
IEXP and DEXP, and unlike TI, does not suffer problems with bias with
even with relatively large intervals.  With a given amount of
simulation time, the results are relatively insensitive to the
intervals chosen, as seen by arriving at the same answer with 11 and
38 state sets.  

Using MBAR, relatively small numbers of states are required to
calculate accurate free energies, even for transformations generally
considered large, such as the solvation of 3-methylindole.  The values
calculated with BAR are almost identical to those computed with MBAR.
However, the uncertainties computed with MBAR more accurately describe
the uncertainties obtained with a bootstrap variance, or computed
directly from repeated samplings.  BAR has a tendency to underestimate
the variance because of assumptions in independence in the error
propagation.  We note that as overlap between states becomes greater,
MBAR should provide greater efficiency; but at the levels of overlap
sufficiently low to give high accuracy, there does not appear to be a
substantial efficiency gain for BAR over MBAR.  The advantages are
instead in the calculation of the uncertainties.

We note that MBAR does take more time than any of the other methods.
For example, on a MacBook 2.1 GHz Core 2 Duo with 4 GB RAM, MBAR for
the entire 38 state dataset took 131 seconds, while BAR took 4.3
second, with both exponential averaging methods and TI taking under a
second.  For the 11 state data set, MBAR took 56 seconds, while BAR
took 7 seconds, and again TI and exponential averaging taking under a
second.  However, this time is dwarfed by the time required to
actually generate the hundreds of nanoseconds of simulation.
Additionally, for all methods, it took 74 seconds to load the 11 state
data set and 122 second to load the 38 state data set using the
example scripts released with the distribution, so the actual time
required for analysis of MBAR was never more than twice that required
for any other method.

Finally, we note some caveats in the statement that MBAR is
necessarily always the best method for determining free energies.  For
solvation of rigid solutes, Gibbs ensemble based methods may be more
efficient, but for, larger more flexible molecules, or for complicated
processes such as protein-ligand binding, alchemical methods are
required.  In certain cases, TI can take advantage of smooth
$\expect{\dhdl}$ curves with properly chosen quadrature values, and
can easily be both low bias, and more accurate, but MBAR has very good
performance over a wide range of cases.

\subsubsection{Entropy-Enthapy Decomposition}
Frequently, researchers are interested in the decomposition of Gibbs
or Helmholtz free energy differences into entropy and ethalpy
differences, where $\Delta G = \Delta H - T\Delta S$, or for canonical
systems, $\Delta A = \Delta U - T\Delta S$.  The most common methods
are for calculating this decomposition are by calculating the free
energy at multiple temperatures and approximating the derivative
$\Delta S = \frac{Delta G}{\Delta T}$, and by an endpoint difference
method, calculating the average internal energy $\expect{U}$ or
enthalpy $\expect{U+PV}$ at two endpoints and subtracting to find the
difference between states.  Both methods are problematic, since the
finite difference method results. 

We note that MBAR allows a better estimator of the average potential
energy at the endstates by using all of the states sampled, not just
the endstate.  Consequently, the variances are about 30-40\% lower in
this particular case for the potential average determined with MBAR
than with simply the endpoints.

We also note that because MBAR gives the covariances between $f$ and
$\expect{u}$, as well as the variances in $f$ and $\expect{u}$, we can
propagate the direct uncertainty in $s$, rather than relying on an
assumption of independence variances of $f$ and $u$, which is clearly
untrue, as both calculations use the same data. Specifically, for 
$s = f-u$, we have $\var{s} = \var{f-u} = \var{f} + \var{u} - 2\Cov{f,u}$
Although useful, we do note that in many cases, the total uncertainty
in the $s=\expect{u}-f$ is dominated by the uncertainty in
$\expect{u}$, and the correlation between $f$ and $u$ is relatively
unmportant in the final calculation.
\begin{table}\renewcommand{\arraystretch}{0.75}
\begin{tabular}{lcccc}
                &$\Delta G$                   &            $\Delta H$        &            $T\Delta S$  \\
                       \multicolumn{4}{c}{11 states}                                      \\
Endpoint method &18.990 $\pm$  (0.063,0.068)  &  156.859 $\pm$  (1.47,1.44)  &   137.869 $\pm$  (1.44,1.43) \\
MBAR            &18.990 $\pm$  (0.063,0.068)  &  136.918 $\pm$  (1.05,1.00)  &   117.928 $\pm$  (1.05,1.00) \\
                        \multicolumn{4}{c}{38 states}                                    \\
Endpoint method &18.938 $\pm$  (0.080,    )   &  131.970 $\pm$  (3.61,4.85)  &   113.032 $\pm$  (3.61,4.86) \\  
MBAR            &18.938 $\pm$  (0.080,    )   &  126.847 $\pm$  (2.10,2.30)  &   107.909 $\pm$  (2.10,2.31) \\
\end{tabular}
\caption{\label{table.bootstrap}Free energy, entropy, and enthalpy differences for
  3-methylindole in kJ/mol.  In both cases, free energies are determined by
  MBAR.  Potential energy difference expectations determined by the
  endpoint method and by MBAR, with entropies determined by
  subtraction.  The first variances use analytical error estimates;
  the second variances in the list for each calculation are bootstrap
  variances over 40 bootstrap samples.}
\end{table}

% MRS: This really bugs me -- why are the endpoint and MBAR free
% energies so far apart? Why are the answers from MBAR for 11 states
% and 38 states so far apart, but are not reflected in the
% uncertainties?  This is exactly the same data.  Analytical variance
% and bootstrap calculations are comparable.  The bias of both methods
% should be very low.
% Trying to understand this better.  Below, I compute: 
%
%\sum A_n[SN[k]:SN[k+1] * W_n[SN[k]:SN[k+1],l] / W_n[SN[k]:SN[k+1],l] 
% Where SN[k] = sum(N_k[0:k].
%
% The 'weight' column is W_n[SN[k]:SN[k+1],l].  
% 
%For 10, contrib from  6: -11692.00019, weight:    0.00009
%For 10, contrib from  7: -11687.01181, weight:    0.02636
%For 10, contrib from  8: -11689.26618, weight:    0.23466
%For 10, contrib from  9: -11683.57593, weight:    0.31940
%For 10, contrib from 10: -11672.67210, weight:    0.41949
%
%For a final answer of: -11680.36433
%
%The direct endpoint average is: -11672.44030
%

%So clearly, the reason they don't agree is because of the fact that
%the other nearby points are not actually giving the right average.
%Is this a big deal?  What does it mean?  I don't really know.  There
%really shouldn't be that much difference.  The analytical error
%estimate for the overall expectation is about 0.5 kT, so the
%different averages shouldn't be off by much more than that if they
%all have similar weights!


As examples that can be modified by users to compute free energies
with MBAR, the first 400 ps steps of each file in the the 38 step data
set are included in
examples/alchemical-free-energy/data/3-methylindole-11steps, as are
the first 500 ps of the 11 step set in
examples/alchemical-free-energy/data/3-methylindole-38steps.

TI, IEXP, DEXP, pairwise BAR, and MBAR can be computed from these
truncated data sets using the script alchemical-gromacs.py included in
examples/alchemical-free-energy.

The full data sets are available as separate downloads, as are the
gromacs input files from simtk.org.


% LocalWords:  MBAR alchemical TI Zwanzig IEXP DEXP solute Coulombic Lennard AA
% LocalWords:  Columbic intramolecular Gromacs FFT subsampled subsampling kJ Av
% LocalWords:  resampled resampling resamplings postprocesing Unc solutes py
% LocalWords:  gromacs simtk org MacBook GHz Duo
