[gmx-users] problem with free energy calculations
Michael Brunsteiner
mbx0009 at yahoo.com
Tue Feb 6 06:31:14 CET 2007
Dear all,
I seem to have a problem with my free energy TI calculations
using Gromacs, I'd be grateful for any suggestions about what
could have gone wrong there...
I calculated the free energy of solvation of water
(basically its chemical potential). The literature value for TIP3P
water, which I used, is about -6.3 kcal/mol -
I get something like +2.5 (!) ...
my calculations are fairly standard, I use PME with
conservative parameters, soft core potentials,
do 21 windows, each one lasting 500 pico seconds, then
after discarding the first 100 ps of each dgdl*xvg file
i use the averages and do a numerical integration,
the works ...
I also use two thermostates (one for the water that disappears,
one for the rest)
I see that the fluctuations of my dU/dlambda values are fairly
large ... can it be that the soft core parameters i use (sc_alpha=0.7,
sc_power=1, sc_sigma=0.3) are inappropriate if, as i do here,
i turn on/off the VdW and Coulomb interactions simultanioulsy ??
mdp plus topology files and the results are included below,
thanks in advance for any suggestions!
cheers,
michael
================ md.mdp ==================================
title = water mu
cpp = /usr/bin/cpp
;
integrator = md
nsteps = 250000
dt = 0.002
;
nstxtcout = 20
nstenergy = 20
nst_log = 20
nstxout = 10000
nstvout = 10000
;
nstlist = 12
nstype = grid
pbc = xyz
nstcomm = 100
rlist = 1.2
;
coulombtype = PME
rcoulomb = 1.2
vdwtype = shift
rvdw = 1.2
;
Tcoupl = nose-hoover
tc_grps = w1 rest
ref_t = 300.0 300.0
tau_t = 0.1 0.1
gen_vel = yes
;
Pcoupl = no
;
pme_order = 4
ewald_rtol = 1e-05
ewald_geometry = 3d
epsilon_surface = 0
optimize_fft = no
fourier_nx = 64
fourier_ny = 64
fourier_nz = 64
;
constraints = hbonds
unconstrained-start = no
lincs_iter = 4
;
free_energy = yes
init_lambda = 0.00 ; (0.05, 0.10, etc in the other runs)
delta_lambda = 0.0
sc_alpha = 0.7
sc_power = 1
sc_sigma = 0.3
=== topol.top =====================================================
; Include forcefield parameters
#include "ffamber03.itp"
; a TIP3P water with zero charge and no VdW interactions in the B
; topology the dummy atom types dumow and dumhw are defined in the
; ff*.atp and ff*nb.itp files used here
[ moleculetype ]
; molname nrexcl
wdummy 1
[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB chargeB massB
1 amber99_42 1 SOL OW 1 0.000 16.00000 dumow 0.00 16.00000
2 amber99_27 1 SOL HW1 1 0.000 1.00800 dumhw 0.00 1.00800
3 amber99_27 1 SOL HW2 1 0.000 1.00800 dumhw 0.00 1.00800
[ settles ]
; OW funct doh dhh
1 1 0.09572 0.15139
[ exclusions ]
1 2 3
2 1 3
3 1 2
; the other waters
; Include water topology
#include "ffamber_tip3p.itp"
; the system
[ system ]
; Name
462 tip3p water
[ molecules ]
; Compound #mols
wdummy 1
SOL 461
=== the results (averages from dgdl*xvg files) ===
# lamda <dU/dl> rmsd
0.00 -4.85 14.04
0.05 -5.44 14.42
0.10 -5.66 14.40
0.15 -6.77 14.93
0.20 -7.98 17.11
0.25 -9.17 17.85
0.30 -11.09 19.52
0.35 -14.23 21.80
0.40 -18.69 27.28
0.45 -27.41 29.40
0.50 -27.75 22.81
0.55 -26.55 17.26
0.60 -20.99 11.85
0.65 -15.47 8.66
0.70 -10.12 6.01
0.75 -6.54 4.53
0.80 -2.75 3.30
0.85 -0.02 2.50
0.90 2.37 1.93
0.95 3.78 1.45
1.00 5.07 1.15
____________________________________________________________________________________
Get your own web address.
Have a HUGE year through Yahoo! Small Business.
http://smallbusiness.yahoo.com/domains/?p=BESTDEAL
More information about the gromacs.org_gmx-users
mailing list