[gmx-users] Alchemical free energy calculations - right choice of lambda points? What else might I do wrong?

[Julian Zachmann]

Dear Gromacs-Users,

I calculate relative binding energies using free energy perturbation (FEP)
and mbar. I have the binding data for 4 receptors and two ligands which are
almost identical. The only difference is the change from a methyl group to
a hydrogen. I alchemically convert a methyl group into one hydrogen atom
with 3 dummy atoms - inside the receptor and in water. First I have several
equilibration steps; the production run takes 4ns. Unfortunately my results
don't match the experimental results and I would like to ask for input
about what I might do wrong.

First the results:

Experimental | Theoretical (kcal) mbar calculated with pymbar
Receptor 1:    0.85 |  0.57
Receptor 2:    0.43 | -1.94
Receptor 3:    1.13 |   0.78
Receptor 4:    2.41 |   0.29

I am using the following lambda points:
;                                       0          1          2         3
       4           5          6          7         8          9         10
       11        12        13        14        15        16            17
      18         19
fep_lambdas             =  0.0000 0.2000 0.4000 0.5000 0.6000 0.7000 0.8000
0.9000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1.0000 1.0000 1.0000
vdw_lambdas             =  0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.2000 0.4000 0.6000 0.7000 0.8000 0.9000 0.9300 0.9600
0.9900 0.9950 1.0000

I think this might be a mistake as the graphic of free energy differences
<https://drive.google.com/open?id=0B2M9aqeJrxnYMWRXbGZ4dkxrb00> shows there
is a huge difference between lambda points 7-8 - when fep_lambdas turn from
0.9 to 1.0. I think of either introducing another lambda point such as
fep_lambdas = 0.95 or to split up fep_lambdas and converge bonded_lambdas,
coul_lambdas, mass_lambdas differently fast to 1.0. I think
temperature_lambdas and restraint_lambdas are not important for my system
so I my gradually convert them from 0.0 to 1.0 over all 20 lambda points. I
would be very thankful if you could advise me how it would be best to
converge the bonded-, coul-, mass_lambdas, and vdw_lambdas to 1.0.

Other guesses about what I might do wrong are in the mol.itp
<https://drive.google.com/file/d/0B2M9aqeJrxnYeGFBMl9XTmprZXM/view?usp=sharing>
 file:

I am not 100% sure, but I think it is correct to put the bond, angle, and
dihedral strengths for dummy atoms to zero.

[ bonds ]
;  ai    aj funct  r  k
(...)
   16    18     1  1.0920e-01  2.8225e+05  1.0920e-01  2.8225e+05
   13    14     1  1.0910e-01  2.8342e+05  1.0910e-01  2.8342e+05
   13    15     1  1.0910e-01  2.8342e+05  1.0910e-01  2.8342e+05
    8     9      1  1.0910e-01  2.8342e+05  1.0910e-01  0.0000
    8    10     1  1.0910e-01  2.8342e+05  1.0910e-01  0.0000
    8    11     1  1.0910e-01  2.8342e+05  1.0910e-01  0.0000
    7    12     1  1.0330e-01  3.0878e+05  1.0330e-01  3.0878e+05
    4     5      1  1.0910e-01  2.8342e+05  1.0910e-01  2.8342e+05
    4     6      1  1.0910e-01  2.8342e+05  1.0910e-01  2.8342e+05
(...)

[ angles ]
;  ai    aj    ak funct  theta   cth
(...)
   13    16    18     1  1.1005e+02  3.8802e+02  1.1005e+02  3.8802e+02
   12     7    13     1  1.1011e+02  3.8652e+02  1.1011e+02  3.8652e+02
   10     8    11     1  1.1074e+02  3.2669e+02  1.1074e+02  0.0000
    9     8    10     1  1.1074e+02  3.2669e+02  1.1074e+02  0.0000
    9     8    11     1  1.1074e+02  3.2669e+02  1.1074e+02  0.0000
    8     7    12     1  1.1011e+02  3.8652e+02  1.0811e+02  3.3907e+02
    7     8     9     1  1.0791e+02  4.1020e+02  1.0791e+02  0.0000
    7     8    10     1  1.0791e+02  4.1020e+02  1.0791e+02  0.0000
    7     8    11     1  1.0791e+02  4.1020e+02  1.0791e+02  0.0000
    7    13    14     1  1.0791e+02  4.1020e+02  1.0791e+02  4.1020e+02
    7    13    15     1  1.0791e+02  4.1020e+02  1.0791e+02  4.1020e+02
(...)

[ dihedrals ]
;i  j   k  l     func   C0  ...  C5
(...)
    12   7    13   15     3     0.65270     1.95811     0.00000    -2.61082
    0.00000     0.00000       0.65270     1.95811     0.00000    -2.61082
  0.00000     0.00000
    12   7    13   16     3     0.65270     1.95811     0.00000    -2.61082
    0.00000     0.00000       0.65270     1.95811     0.00000    -2.61082
  0.00000     0.00000
    11   8    7    12     3     0.65270     1.95811     0.00000    -2.61082
    0.00000     0.00000       0.00000     0.00000     0.00000     0.00000
  0.00000     0.00000
    11   8    7    13     3     0.65270     1.95811     0.00000    -2.61082
    0.00000     0.00000       0.00000     0.00000     0.00000     0.00000
  0.00000     0.00000
    10   8    7    12     3     0.65270     1.95811     0.00000    -2.61082
    0.00000     0.00000       0.00000     0.00000     0.00000     0.00000
  0.00000     0.00000
    10   8    7    13     3     0.65270     1.95811     0.00000    -2.61082
    0.00000     0.00000       0.00000     0.00000     0.00000     0.00000
  0.00000     0.00000
    9    8    7    12     3     0.65270     1.95811     0.00000    -2.61082
    0.00000     0.00000       0.00000     0.00000     0.00000     0.00000
  0.00000     0.00000
    9    8    7    13     3     0.65270     1.95811     0.00000    -2.61082
    0.00000     0.00000       0.00000     0.00000     0.00000     0.00000
  0.00000     0.00000
    8    7    13   14     3     0.65270     1.95811     0.00000    -2.61082
    0.00000     0.00000       0.65270     1.95811     0.00000    -2.61082
  0.00000     0.00000
    8    7    13   15     3     0.65270     1.95811     0.00000    -2.61082
    0.00000     0.00000       0.65270     1.95811     0.00000    -2.61082
  0.00000     0.00000
(...)


My feeling tells me, that mdp files
<https://drive.google.com/open?id=0B2M9aqeJrxnYdDdYYUJoTmk3MDA> are
correct. I sticked close to the tutorial
<http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/01_theory.html>
on
the Gromacs page.

I would be very thankful for any input you could give me!

Best wishes,
Julian