[gmx-users] LIE precision and accuracy

[David van der Spoel]

On 12/17/09 9:24 PM, Marc Charendoff wrote:
> 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.
Cut-off should not be used either :).

You can estimate the error rather than std. dev. using g_analyze -ee


> 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
>


-- 
David.
________________________________________________________________________
David van der Spoel, PhD, Professor of Biology
Dept. of Cell and Molecular Biology, Uppsala University.
Husargatan 3, Box 596,  	75124 Uppsala, Sweden
phone:	46 18 471 4205		fax: 46 18 511 755
spoel at xray.bmc.uu.se	spoel at gromacs.org   http://xray.bmc.uu.se/~spoel
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++