[gmx-users] Alchemical free energy calculations - right choice of lambda points? What else might I do wrong?
hannes.loeffler at stfc.ac.uk
hannes.loeffler at stfc.ac.uk
Thu Jul 30 08:29:23 CEST 2015
Hi,
you wouid _not_ set bonded terms to zero. That would mean that you have essentially a non-bonded end state allowing the atom to be "everywhere" (reading Stefan Boresch' papers from 1999 on this may be helpful). The dihedrals are debatable (a collaborator of mine thinks that having them zero may help in replicate methods or so). Our strategy with FESetup (http://www.hecbiosim.ac.uk/fesetup) is to copy those terms over from the state without the dummies.
The mass-lambdas are best kept at one end state as discussed previously to avoid any bad interactions with constraints (you say you have X-C-H3 to X-H-DU3).
Cheers,
Hannes.
________________________________________
From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se [gromacs.org_gmx-users-bounces at maillist.sys.kth.se] on behalf of Julian Zachmann [FrankJulian.Zachmann at uab.cat]
Sent: 29 July 2015 14:10
To: gmx-users at gromacs.org
Subject: [gmx-users] Alchemical free energy calculations - right choice of lambda points? What else might I do wrong?
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
--
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