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-requ...@gromacs.org.