In spite of the fact that I have never used g_lie, I have used a
variety of different techniques to calculate binding free energies,
and I would be astounded if it was possible to get the precision that
you desire with only 1 ns of sampling per state. Are you sure that the
reported imprecision is different from the actual imprecision obtained
via 1 ns of sampling? Generally, one evaluates the convergence by
block averaging or some other convergence indicator. If you do find a
measure by which your precision is reported to be very good, then I
would suspect that your sampling is trapped and that your free
energies are not accurate in any event. In short, I think that your
comment about "rendering the energy difference meaningless" is
actually quite astute.
Chris.
-- original message --
Hello,
I am in the middle of using Gromacs 3.3.3 to evaulate
interaction energies between a protein and two separate ligands to
compare their binding affinities. In doing this I have run each
solvated complex through a 1 ns MD simulation and each solvated ligand
through a 1 ns MD simulation. After the trajectories were complete, I
ran g_lie on the trajectories for each case. Next, I subtracted the
average ligand LIE from the average complex LIE for each case to get
the interaction energy . Data is as follows:
Ligand 1 Ligand 2
LIE (kcal/mol) -5.7 -8.9
St dev. 6.7 4.7
While the average energies themselves are reasonable, the st dev's
seem to me to be too large - rendering the energy difference
meaningless. I calculated the st. dev's based on the trajectories'
g_lie values (in this case 1000 data points for a 1 ns run). I read in
the tutorials that PME condition should not be used with g_lie, so I
switched to cut-offs for my electrostatic forces. Is there anything
else I might be missing? Again, I know the st. dev's are just too high
to make g_lie useful so I blame myself for the initial setup. My .mdp
file is below.
Thanks, Marc
; mdrun for C4AHL on unix parameters
title = trp_drg MD
cpp = /usr/bin/cpp ; location of cpp on linux
constraints = all-bonds
integrator = md
dt = 0.002 ; ps!
nsteps = 500000 ; total 1000 ps
nstcomm = 1
nstxout = 500 ; output coords every 1 ps
nstvout = 0
nstfout = 0
nstenergy = 500
nstlist = 10
ns_type = grid
pbc = xyz
rlist = 0.9
coulombtype = cut-off
rcoulomb = 1.2
vdwtype = cut-off
rvdw = 1.4
;fourierspacing = 0.12
;fourier_nx = 0
;fourier_ny = 0
;fourier_nz = 0
;pme_order = 6
;ewald_rtol = 1e-5
;optimize_fft = yes
; Berendsen temperature coupling is on in five groups
Tcoupl = berendsen
tau_t = 0.1 0.1
tc-grps = protein non-protein
ref_t = 300 300
; Use Energy group monitoring
energygrps = protein C6L OHY SOL Na+
; Pressure coupling is on
Pcoupl = berendsen
pcoupltype = isotropic
tau_p = 1.0
compressibility = 4.5e-5
ref_p = 1.0
; Generate velocities is on at 300 K
gen_vel = yes
gen_temp = 300.0
gen_seed = -1
free_energy = yes
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
