[gmx-users] Umbrella Sampling - Protein Ligand

[Steven Neumann]

Dear Gmx Users, Dear Justin,

I pulled my ligand away from my protein. Ligand was attached to lower part
of my protein, I pulled in Z coordinate it using:


; Run parameters
integrator = md ; leap-frog integrator
nsteps = 5000000 ; 2 * 5000000 = 10 ns
dt = 0.002 ; 2 fs
tinit = 0
nstcomm = 10
; Output control
nstxout = 50000 ; save coordinates every 100 ps
nstvout = 50000 ; save velocities every
nstfout = 5000
nstxtcout = 5000 ; every 10 ps
nstenergy = 5000
; Bond parameters
continuation = yes ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained
; Neighborsearching
ns_type = grid ; search neighboring grid cells
nstlist = 5 ; 10 fs
rlist = 0.9 ; short-range neighborlist cutoff (in nm)
rcoulomb = 0.9 ; short-range electrostatic cutoff (in nm)
rvdw = 0.9 ; short-range van der Waals cutoff (in nm)
ewald_rtol = 1e-5 ; relative strenght of the Ewald-shifted potential
rcoulomb
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.12 ; grid spacing for FFT
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
optimize_fft = yes
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc_grps = Protein_LIG Water_and_ions ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 298 298 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 1.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; assign velocities from Maxwell distribution
; These options remove COM motion of the system
; Pull code
pull = umbrella
pull_geometry = distance
pull_dim = N N Y
pull_start = yes
pull_ngroups = 1
pull_group0 = Protein
pull_group1 = LIG182
pull_init1 = 0
pull_rate1 = 0.0
pull_k1 = 200 ; kJ mol^-1 nm^-2
pull_nstxout = 1000 ; every 2 ps
pull_nstfout = 1000 ; every 2 ps

Following Justin's tutorial I used perl script to extract coordinate for
each window.
0       2.4595039
1       2.4745028
...
500    8.74

My ligand at the begining was at such distance as it was in the lower part
of the protein. Then I used 0.1 nm spacing at the begining (till 4 nm) and
0.2 nm later on.
And following equilibration in each window I run umbrella sampling for 10ns
in app 49 windows:
Run parameters
integrator = md ; leap-frog integrator
nsteps = 5000000 ; 2 * 5000000 = 10 ns
dt = 0.002 ; 2 fs
tinit = 0
nstcomm = 10
; Output control
nstxout = 50000 ; save coordinates every 100 ps
nstvout = 50000 ; save velocities every
nstfout = 5000
nstxtcout = 5000 ; every 10 ps
nstenergy = 5000
; Bond parameters
continuation = yes ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained
; Neighborsearching
ns_type = grid ; search neighboring grid cells
nstlist = 5 ; 10 fs
rlist = 0.9 ; short-range neighborlist cutoff (in nm)
rcoulomb = 0.9 ; short-range electrostatic cutoff (in nm)
rvdw = 0.9 ; short-range van der Waals cutoff (in nm)
ewald_rtol = 1e-5 ; relative strenght of the Ewald-shifted potential
rcoulomb
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.12 ; grid spacing for FFT
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
optimize_fft = yes
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc_grps = Protein_LIG Water_and_ions ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 298 298 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 1.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; assign velocities from Maxwell distribution
; These options remove COM motion of the system
; Pull code
pull = umbrella
pull_geometry = distance
pull_dim = N N Y
pull_start = yes
pull_ngroups = 1
pull_group0 = Protein
pull_group1 = LIG182
pull_init1 = 0
pull_rate1 = 0.0
pull_k1 = 200 ; kJ mol^-1 nm^-2
pull_nstxout = 1000 ; every 2 ps
pull_nstfout = 1000 ; every 2 ps

My PMF profile:
http://speedy.sh/zerqZ/profile.JPG

My histogram: http://speedy.sh/PyhnN/Histo.JPG

Why g_wham takes into account distances below 2.45 nm as the 1st structure
was at 2.45. If I get rid of the distances below 2.45 (those weird values
PMF values) I obtain beautiful profile:

http://speedy.sh/TUXGC/profile1.JPG

Please, explain!
Steven
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20120308/0a34f901/attachment.html>