[gmx-users] Free Energy Calculations of Octanol in Octanol, What is the Issue?
Justin Lemkul
jalemkul at vt.edu
Wed Jun 27 23:32:11 CEST 2018
On 6/26/18 8:48 PM, Dallas Warren wrote:
> 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"?
That's what I see, yes.
-Justin
> 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.
>>
--
==================================================
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