[gmx-users] Free Energy Calculations of Octanol in Octanol, What is the Issue?
Dallas Warren
dallas.warren at monash.edu
Wed Jun 27 02:48:57 CEST 2018
So I take it then (and to confirm this) that it is an issue with the output
from gmx bar to the histogram file, rather than results from the lambda
state simulations? But the results it has presented in the results table
and final value for the free energy change is "correct"?
Previously I'd only looked at the final result, but thought much better to
make sure all the underlying data is correct and valid too.
The number of set of values seem consistent, 21 lambda states, 21 dH/dl
sets between states, plus 2 extra (which I'd assumed were overall dH/dl for
coulomb and vdw terms.
Good to know close values David, thank you.
Mark, we are getting a programmer on the team here, so hopefully we will be
help out with things like this that are important to us in the near future.
On Tue, 26 Jun. 2018, 10:48 pm Justin Lemkul, <jalemkul at vt.edu> wrote:
>
>
> 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
>
> ==================================================
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>
More information about the gromacs.org_gmx-users
mailing list