[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