(Note: These tutorials are meant to provide
illustrative examples of how to use the AMBER software suite to carry out
simulations that can be run on a simple workstation in a reasonable period of
time. They do not necessarily provide the optimal choice of parameters or
methods for the particular application area.)
Copyright Andrew T. Fenley & Michael K. Gilson 2014
Section 4: Computing the Binding Enthalpy
We now describe how to process the potential energy values from the output files of the production simulations described in section 3. The output files are also accessible here: all_prod1_out.tar.gz.
The output files (*prod1.out) contain all the information we need for computing the binding enthalpy value. We will utilize common Unix/Linux shell commands to extract the necessary information. Below are the commandline commands used to estimate the average total potential energy from each production simulation.
grep EPtot water_prod1.out | head -n -2 | awk '{if(NR>5000)sum+=$9}END{printf "%14.6f\n",sum/(NR-5000)}'
grep EPtot b2_prod1.out | head -n -2 | awk '{if(NR>5000)sum+=$9}END{printf "%14.6f\n",sum/(NR-5000)}'
grep EPtot cb7_prod1.out | head -n -2 | awk '{if(NR>5000)sum+=$9}END{printf "%14.6f\n",sum/(NR-5000)}'
grep EPtot cb7_b2_prod1.out | head -n -2 | awk '{if(NR>5000)sum+=$9}END{printf "%14.6f\n",sum/(NR-5000)}'
Here are a few details about the above commands. The head -n -2 ignores the last two potential energy entries from a given output file. We ignore these two lines because (assuming the simulation finished without any errors) they correspond to the average value over the whole trajectory and the fluctuation of the potential energy and therefore should not be included in our average. We then use awk to skip the first 5000 lines of the output to act as an additional equilibration time which accounts for switching from the NPT to the NVT ensemble. awk also adds all the potential energy values and then prints the final average value. Below are the average values of the total potential energy for each system. These values should be within ~0.2 kcal/mol of anyone rerunning the simulations as described in this tutorial.
| Simulation | Potential Energy (kcal/mol) |
|---|---|
| Pure Water | -14312.4 |
| B2 | -14283.7 |
| CB7 | -15226.2 |
| CB7 and B2 | -15218.9 | Binding Enthalpy (kcal/mol) |
| ΔH | -21.4 |
Estimating the Uncertainty of the Binding Enthalpy Value:
There are two popular approaches for estimating the uncertainty in the mean potential energy per simulation, blocking[2] and statistical inefficiency[3] analysis. Either approach should yield very similar estimates on the standard deviation of the mean potential energies. The Supporting Information of Fenley et al. discusses the application of such error analysis to these host-guest systems in more detail.
We calculated the standard deviation of the mean potential energies to be about 0.13 kcal/mol for each of the simulations in this example. To compute the uncertainty in the absolute binding enthalpy (σH), we add the uncertainties of each simulation in quadrature and then take the square root:
With a standard deviation of the mean potential energy equal to 0.13 kcal/mol for all the simulations, the total uncertainty of the binding enthalpy is about 0.3 kcal/mol. Thus, the binding enthalpy of guest B2 with CB7 is predicted to be -21.4 ± 0.3 kcal/mol which agrees quite well with the estimate published in Fenley et al.
At this point, we are done with the computation part of the tutorial. However, we recommend reading the last section as it contains a summary of potential pitfalls one may encounter in computing absolute or even relative binding enthalpies via the direct method. Additional information is located in the Supporting Information of Fenley et al.
Click here to go back to section 3
References
[2] Flyvbjerg, H. and Petersen, H. G., "Error Estimates on Averages of Correlated Data," J. Chem. Phys., 1989, 91, 461-466
[3] Chodera, J. D., Swope, W. C., Pitera, J. W., Seok, C., and Dill, K. A., "Use of the weighted histogram analysis method for the analysis of simulated and parallel tempering simulations," J. Chem. Theory and Comput., 2007 3(1), 26-41
(Note: These tutorials are meant to provide
illustrative examples of how to use the AMBER software suite to carry out
simulations that can be run on a simple workstation in a reasonable period of
time. They do not necessarily provide the optimal choice of parameters or
methods for the particular application area.)
Copyright Andrew T. Fenley & Michael K. Gilson 2014