[gmx-users] Free Energy Calculations of Octanol in Octanol, What is the Issue?
Mark Abraham
mark.j.abraham at gmail.com
Tue Jun 26 15:13:49 CEST 2018
Hi,
gmx bar got a last-minute rewrite before GROMACS 4.6 (some implementations
of some of the free-energy algorithms changed) and I'm sure it's still in
need of love. However, the core team can't do all the things we'd like to
do, so if there are people who want to have working alternative
implementations of BAR analysis code, then we'd particularly love to have
inputs and trusted results (e.g. from 4.5.x g_bar, or some other tool) so
that we can use those to frame tests.
Mark
On Tue, Jun 26, 2018 at 2: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 re
> <https://maps.google.com/?q=files+for+different+data+types+and+to+re&entry=gmail&source=g>move
> 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