[gmx-users] large free energy difference of mutating ATP to GTP in solution using FEP method in Gromacs-4.6.5
Michael Shirts
mrshirts at gmail.com
Sat May 24 19:06:31 CEST 2014
There are many ways for free energy calculations to go wrong. Look at
http://www.alchemistry.org for some discussions.
One thing I noticed (I didn't have a chance to look at everything) -
you only talk about one simulation. When doing classical simulations,
then you can't just mutate atoms, because classical simulations
neglect quantum mechanical terms. You can only look at differences in
alchemical simulations. For example, you can look at the differences
between solvation free energies by performing the mutation in both
water and vacuum. You have to think VERY carefully about the
thermodynamic end states that you are using.
On Sat, May 24, 2014 at 1:18 AM, dbaogen <dbaogen at gmail.com> wrote:
> 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
> --
> 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