[gmx-users] Free Energy Calculations of Octanol in Octanol, What is the Issue?
Justin Lemkul
jalemkul at vt.edu
Tue Jun 26 14:48:03 CEST 2018
On 6/25/18 10:34 PM, Dallas Warren wrote:
> I am trying to identify where an issue is with my free energy
> calculations of an octanol molecule in octanol (and whether it is
> actually an issue or not).
>
> The lambda states used are in 0.1 increments from 0 to 1, with coulomb
> interactions turned off first, then van der Waals. Settings within the
> mdp file shown below.
>
> ; Free Energy
> free_energy = yes ; free energy calculation
> init-lambda-state = 0
> coul_lambdas = 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
> 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
> vdw_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
>
> After performing each of those 21 simulations, I run the gmx bar
> script on the outputs. The final results from that analysis is:
>
> Temperature: 298 K
>
> Detailed results in kT (see help for explanation):
>
> lam_A lam_B DG +/- s_A +/- s_B +/- stdev +/-
> 0 1 2.57 0.07 0.29 0.06 0.35 0.07 0.78 0.05
> 1 2 1.95 0.08 0.27 0.10 0.32 0.10 0.95 0.03
> 2 3 1.06 0.10 0.57 0.04 0.49 0.05 0.98 0.03
> 3 4 0.43 0.06 0.15 0.05 0.09 0.05 0.56 0.03
> 4 5 0.24 0.04 0.09 0.04 0.07 0.03 0.38 0.03
> 5 6 0.13 0.02 0.04 0.02 0.03 0.02 0.26 0.02
> 6 7 0.08 0.01 0.02 0.01 0.02 0.01 0.19 0.01
> 7 8 0.05 0.00 0.01 0.01 0.01 0.01 0.16 0.01
> 8 9 0.03 0.00 0.01 0.00 0.01 0.00 0.14 0.01
> 9 10 0.01 0.00 0.01 0.00 0.01 0.00 0.14 0.01
> 10 11 2.30 0.01 0.25 0.01 0.34 0.01 0.77 0.01
> 11 12 2.24 0.01 0.31 0.01 0.43 0.01 0.82 0.00
> 12 13 2.13 0.01 0.32 0.01 0.46 0.01 0.87 0.01
> 13 14 1.98 0.02 0.41 0.01 0.60 0.02 0.96 0.01
> 14 15 1.77 0.03 0.47 0.01 0.70 0.02 1.06 0.01
> 15 16 1.48 0.01 0.60 0.02 0.93 0.03 1.21 0.01
> 16 17 0.87 0.02 0.92 0.01 1.53 0.01 1.54 0.01
> 17 18 -0.55 0.05 1.68 0.05 2.86 0.06 2.29 0.03
> 18 19 -2.11 0.03 1.70 0.04 1.84 0.04 2.12 0.03
> 19 20 -1.20 0.02 0.64 0.01 0.59 0.01 1.11 0.01
>
>
> Final results in kJ/mol:
>
> point 0 - 1, DG 6.37 +/- 0.17
> point 1 - 2, DG 4.84 +/- 0.19
> point 2 - 3, DG 2.63 +/- 0.26
> point 3 - 4, DG 1.06 +/- 0.14
> point 4 - 5, DG 0.60 +/- 0.09
> point 5 - 6, DG 0.32 +/- 0.04
> point 6 - 7, DG 0.19 +/- 0.02
> point 7 - 8, DG 0.12 +/- 0.01
> point 8 - 9, DG 0.08 +/- 0.00
> point 9 - 10, DG 0.02 +/- 0.00
> point 10 - 11, DG 5.71 +/- 0.02
> point 11 - 12, DG 5.55 +/- 0.03
> point 12 - 13, DG 5.28 +/- 0.03
> point 13 - 14, DG 4.91 +/- 0.05
> point 14 - 15, DG 4.40 +/- 0.07
> point 15 - 16, DG 3.67 +/- 0.03
> point 16 - 17, DG 2.16 +/- 0.04
> point 17 - 18, DG -1.37 +/- 0.12
> point 18 - 19, DG -5.22 +/- 0.08
> point 19 - 20, DG -2.97 +/- 0.05
>
> total 0 - 20, DG 38.35 +/- 0.45
>
> A combined plot of the bar.xvg and barint.xvg data can be found at
> https://twitter.com/dr_dbw/status/1011427248807632896
>
> Everything seems reasonable from that point.
>
> However, when look into the histogram.xvg file generated, in some of
> the calculations between states something weird is going on, and I'm
> trying to identify what the actually issue is, if it is a particular
> lambda state etc, and having a bit of difficulty with that.
>
> Looking at the second set of values in histogram.xvg, which from what
> I can tell is the dH/dl for the vdw can be found at
> https://twitter.com/dr_dbw/status/1011429960173490176 That looks as
> expected.
>
> Now if look at the first set of values in histogram.xvg, which is
> dH/dl for coulomb?, a graph can be found at
> https://twitter.com/dr_dbw/status/1011430685695864834 , there are
> obviously some extreme values, zooming in some some of the sets look
> okay, but obviously some are not
> https://twitter.com/dr_dbw/status/1011432844684451841
>
> I'm at bit of a loss from here to work out what is going on.
>
> To try and work out what is going on, I have plotted the dVcoul/dl (
> https://twitter.com/dr_dbw/status/1011434677905731584 ) and dVvdw/dl (
> https://twitter.com/dr_dbw/status/1011434327253516288 ) values for
> each of the simulations, to see if something stands out there as being
> wrong.
I ran into this exact same thing yesterday as I was re-writing one of my
tutorials for the 2018 version. The histogram file produced by gmx bar
is just a mess of different data sets, many of which have the same exact
legend, and including some that just correspond to the number of data
points, for which I do not see any purpose.
FWIW, I don't think there's anything obviously wrong with your results,
but the histogram output is basically unintelligible. This would be a
great idea for a feature request in Redmine, to write different
histogram files for different data types and to remove unnecessary
"histograms" with only a single value.
-Justin
--
==================================================
Justin A. Lemkul, Ph.D.
Assistant Professor
Virginia Tech Department of Biochemistry
303 Engel Hall
340 West Campus Dr.
Blacksburg, VA 24061
jalemkul at vt.edu | (540) 231-3129
http://www.thelemkullab.com
==================================================
More information about the gromacs.org_gmx-users
mailing list