[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