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 dbao...@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 = 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 1.008000
22 N* 1RAG N9 22 -0.02510
