[gmx-users] Free Energy Calculations of Octanol in Octanol, What is the Issue?

[Dallas Warren]

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.

Catch ya,

Dr. Dallas Warren
Drug Delivery, Disposition and Dynamics
Monash Institute of Pharmaceutical Sciences, Monash University
381 Royal Parade, Parkville VIC 3052
dallas.warren at monash.edu
---------------------------------
When the only tool you own is a hammer, every problem begins to resemble a nail.