Results

Gromacs provides an inbuilt utility gmx bar to analyze thermodynamic integration results, using the Bennett acceptance ratio (BAR) method to compute the free energy difference between adjacent lambda values, which are then summed to give \( \Delta G \).  Alternatively, we can examine the output .xvg files containing time series for \( dH/d \lambda \); the time average gives the integrand, with its error bar estimated from the variance and number of points as \( (\sigma^2/N)^{1/2} \).

Figure 4 shows a typical time series for \( dH/d \lambda\) (here, for ethanol at \( \lambda = 0\).  The fluctuations are large and rapid, but the time series is long enough to give a reliable average value.  Figure 5 shows the resulting integrand with error bars.  The integrand varies rapidly at small lambda, which is why we took more closely spaced lambda values there.

Figure 4. dH/dλ of fully coupled ethanol solute (λ = 0).

Figure 5. Integrand of dH/dλ for ethanol from fully coupled state (λ =0) to fully uncoupled state (λ = 1).

The table below shows the solvation free energy \( \Delta G\)  (in kJ/mol) for each alcohol, from the integrand method, gmx bar method, and experiment.

The error bars reflect only statistical error, not systematic error from inaccuracies in the simulation potential itself.  The results from gmx bar and direct integration of \( \langle dH/d \lambda \rangle \) are consistent.  (The small discrepancies between these analyses reflect errors in the integrand near \( \lambda = 0.25 \), where it bends sharply, and more points should be added.)

As expected, solvation becomes steadily less favorable (\( \Delta G\) becomes less negative) from methanol to ethanol to propanol, reflecting the influence of the longer alkyl tail.  The simulation results are systematically too small in magnitude relative to experiment, by about 3 kJ/mol.  This likely reflects an error in the simulation potential; for example, the partial charges on the -OH group may be too small, so that it is less well solvated by water.  Previous work has shown that solvation free energies from simulations using the OPLS-AA potentials and TIP3P water model deviate from experiment, and depend on the choice of water model [1,2].

The following tarball includes Gromacs run and results files and Mathematica analysis notebook used for all alcohols.

henry.tar

References

1. MR Shirts, JW Pitera, WC Swope, VS Pande; Extremely precise free energy calculations of amino acid side chain analogs: Comparison of common molecular mechanics force fields for proteins; The Journal of chemical physics 119(11), 5740-5761
2. T. Wescott, L. R. Fisher, S. Hanna; Use of thermodynamic integration to calculate the hydration free energies of n-alkanes; The Journal of chemical physics 116, 2365 (2002)