[gmx-users] Umbrella Sampling - Ligand Protein

[Justin A. Lemkul]

Steven Neumann wrote:
> 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!
> 

The way you're thinking about distance is not consistent.  Again, this is a 
hazard of trying to map my tutorial onto your problem.  You say you have a 
ligand bound to the "lower part" of your protein, and then you're pulling in the 
z-direction.  The COM distance (as measured by g_dist and extracted using my 
script) is not equivalent to the distance along the reaction coordinate, if that 
reaction coordinate is only one dimension.  In the tutorial, it was.  Here, it 
is not, hence the massive sampling defects that you're observing and 
considerable redundancy in many of your windows.

Check the output of grompp for the actual restraint distances that mdrun will 
interpret.  They are printed to the screen.

-Justin

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================