[gmx-users] large free energy difference of mutating ATP to GTP in solution using FEP method in Gromacs-4.6.5
dbaogen
dbaogen at gmail.com
Sat May 24 07:28:53 CEST 2014
Dear all,
I want to calculate the free energy difference of mutating ATP to GTP in solution using FEP method. Firstly, the hybrid topology and structure files for A (ATP) and B (GTP) state using dummy atoms were constructed. Secondly, the system is running for 10 ns to reach an equilibrium state. And then the structure at 10 ns is as the starting structure to carry out FEP calculation. In the course of FEP, the coulomb interaction was firstly changed, and then the VDW interactions. Total 32 lambda points are set in the mdp file shown in the following:
integrator = sd
nsteps = 100000
dt = 0.002
nstenergy = 1000
nstlog = 5000
nstcalcenergy = 100
nstcomm = 1
cutoff-scheme = group
rlist = 1.2
dispcorr = EnerPres
vdw-type = switch
;cut-off lengths
rvdw = 1.1
rvdw-switch = 1
; Coulomb interactions
coulombtype = pme
rcoulomb = 1.2
fourierspacing = 0.12
; Constraints
constraints = all-bonds
; set temperature to 310K
tcoupl = v-rescale
tc-grps = system
tau-t = 1.0
ref-t = 310
; pressure control
pcoupl = Parrinello-Rahman
ref-p = 1
compressibility = 4.5e-5
tau-p = 0.5
; and set the free energy parameters
free-energy = yes
sc-power = 1
sc-sigma = 0.3
sc-alpha = 0.5
sc-coul = no
sc-r-power = 6
; we still want the molecule to interact with itself at lambda=0
couple-intramol = no
couple-lambda0 = vdw-q
couple-lambda1 = vdw-q
; $LAMBDA$ changed from 0 to 32
init-lambda-state = $LAMBDA$
nstdhdl = 100
calc-lambda-neighbors = 1
fep-lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.01 0.03 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0
;change electrostatic and then LJ interaction
coul-lambdas = 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.00 1.00 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
vdw-lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.01 0.03 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0
atom part of hybrid topology file is :
[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB chargeB
1 O3 1 RAG O1G 1 -0.95260 16.000000
2 P 1 RAG PG 2 1.26500 31.000000
3 O3 1 RAG O2G 3 -0.95260 16.000000
4 O3 1 RAG O3G 4 -0.95260 16.000000
5 OS 1 RAG O3B 5 -0.53220 16.000000
6 P 1 RAG PB 6 1.38520 31.000000
7 O2 1 RAG O1B 7 -0.88940 16.000000
8 O2 1 RAG O2B 8 -0.88940 16.000000
9 OS 1 RAG O3A 9 -0.56890 16.000000
10 P 1 RAG PA 10 1.25320 31.000000
11 O2 1 RAG O1A 11 -0.87990 16.000000
12 O2 1 RAG O2A 12 -0.87990 16.000000
13 OS 1 RAG O5* 13 -0.59870 16.000000
14 CT 1 RAG C5* 14 0.05580 12.000000
15 H1 1 RAG H50 15 0.06790 1.008000
16 H1 1 RAG H51 16 0.06790 1.008000
17 CT 1 RAG C4* 17 0.10650 12.000000
18 H1 1 RAG H40 18 0.11740 1.008000
19 OS 1 RAG O4* 19 -0.35480 16.000000
20 CT 1 RAG C1* 20 0.03940 12.000000 CT 0.01910 12.000000
21 H2 1 RAG H10 21 0.20070 1.008000 H2 0.20060 1.008000
22 N* 1 RAG N9 22 -0.02510 14.000000 N* 0.04920 14.000000
23 CK 1 RAG C8 23 0.20060 12.000000 CK 0.13740 12.000000
24 H5 1 RAG H80 24 0.15530 1.008000 H5 0.16400 1.008000
25 NB 1 RAG N7 25 -0.60730 14.000000 NB -0.57090 14.000000
26 CB 1 RAG C5 26 0.05150 12.000000 CB 0.17440 12.000000
27 CA 1 RAG C6 27 0.70090 12.000000 C 0.47700 12.000000
28 N2 1 RAG N6 28 -0.90190 14.000000 O -0.55970 16.000000
29 H 1 RAG H60 29 0.41150 1.008000 DUM_1 0.000000 1.0080
30 H 1 RAG H61 30 0.41150 1.008000 DUM_1 0.000000 1.0080
31 NC 1 RAG N1 31 -0.76150 14.000000 NA -0.47870 14.000000
32 CQ 1 RAG C2 32 0.58750 12.000000 CA 0.76570 12.000000
33 H5 1 RAG H2 33 0.04730 1.008000 N2 -0.96720 14.000000
34 NC 1 RAG N3 34 -0.69970 14.000000 NC -0.63230 14.000000
35 CB 1 RAG C4 35 0.30530 12.000000 CB 0.12220 12.000000
36 CT 1 RAG C3* 36 0.20220 12.000000
37 H1 1 RAG H30 37 0.06150 1.008000
38 OH 1 RAG O3* 38 -0.65410 16.000000
39 HO 1 RAG H3' 39 0.43760 1.008000
40 CT 1 RAG C2* 40 0.06700 12.000000
41 H1 1 RAG H20 41 0.09720 1.008000
42 OH 1 RAG O2* 42 -0.61390 16.000000
43 HO 1 RAG H2' 43 0.41860 1.008000
; add dummy atoms
44 DUM_1 1 RAG DH1 44 0.000000 1.008000 H 0.34240 1.008000
45 DUM_1 1 RAG DH21 45 0.000000 1.008000 H 0.43640 1.008000
46 DUM_1 1 RAG DH22 46 0.000000 1.008000 H 0.43640 1.008000
After the FEP simulation is done, g_bar program is used to obtain the free energy difference. The value is about -260 kJ/mol which is too large. Because the difference between ATP and GTP is only base part. In addition, if the other parameters are kept and only change sc-coul = yes and couple-intramol = yes, the free energy difference is -160 kJ/mol. So far, I have not found the solution. Could anyone give me some suggestions about this? whether the mdp file or top for FEP calculation is wrong? The equilibrium hybrid structure and topology files are enclosed. Thanks for your great help!
Best
Duan Baogen
Beijing Computational Science Research Center
dbaogen
More information about the gromacs.org_gmx-users
mailing list