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 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 = 10
> 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 1RAGO1G 1 -0.95260 16.00
> 2 P 1RAG PG 21.26500 31.00
> 3 O3 1RAGO2G 3 -0.95260 16.00
> 4 O3 1RAGO3G 4 -0.95260 16.00
> 5 OS 1RAGO3B 5 -0.53220 16.00
> 6 P 1RAG PB 61.38520 31.00
> 7 O2 1RAGO1B 7 -0.88940 16.00
> 8 O2 1RAGO2B 8 -0.88940 16.00
> 9 OS 1RAGO3A 9 -0.56890 16.00
> 10 P 1RAG PA 101.25320 31.00
> 11 O2 1RAGO1A 11 -0.87990 16.00
> 12 O2 1RAGO2A 12 -0.87990 16.00
> 13 OS 1RAGO5* 13 -0.59870 16.00
> 14 CT 1RAGC5* 140.05580 12.00
> 15 H1 1RAGH50 150.06790 1.008000
> 16 H1 1RAGH51 160.06790 1.008000
> 17 CT 1RAGC4* 170.10650 12.00
> 18 H1 1RAGH40 180.11740 1.008000
> 19 OS 1RAGO4* 19 -0.35480 16.00
> 20 CT 1RAGC1* 200.03940 12.00 CT
> 0.01910 12.00
> 21 H2 1RAGH10 210.20070 1.008000 H2
> 0.20060