[gmx-users] md_pull code in umbrella sampling

[Christopher Neale]

It's not US if you have a non-zero pull-rate1 that gets applied during your simulation, although with your pull-geometry=distance I don't think that pull-rate is even used (therefore I think that part is technically correct, although your .mdp file is confusing on this issue).

Careful that X and Y are going to be free to diffuse introducing what will likely be a convergence nightmare for separation distances at which the drug and protein are just coming into contact. Specifically, while you may set up the system with the drug and protein having the same COM in X and Y, that will change over time (more or less, depending on what define = -DPOSRES_Protein does for your system). This likely isn't a problem for large separation distances (although there will be some convergence issues here due to the need to sample over all X/Y at some larger displacements for which coulomb interaction are still non-negligable). Neither is it likely to be a problem for umbrellas in which the protein constrains your drug within a favourable binding pocket. However, in umbrellas where the drug and protein just come into contact, you will probably start with contact, but then need to obtain the proper boltzmann populations for interactions that depend on X/Y. What's worse, if you end up predominantly sampling conformational space that overlaps along your order parameter (Z) but not in some other orthogonal degree of freedom (e.g., X/Y), then your PMF will be simply wrong (that will get fixed with extensive sampling, but the amount of sampling that you need might be unobtainably large).

Personally, I would solve this by adding a second pull group with a flat-bottom potential that acts only on X and Y to keep the drug within a cylinder that extends along Z away from your protein. However, you would need to modify gromacs source code for this (I have some posts on the uses lists over the years about how to do this). An alternative solution is to define an order parameter that acts in all 3 dimensions, although you may need to compute the free energy of releasing this restraint (a) in the binding pocket and (b) at large separations where your PMF flattens out to zero.

Careful also about gen_vel = no (i.e., be sure to load in velocities if you are not generating them)

Chris.
________________________________________
From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of Arunima Shilpi <writetoash28 at gmail.com>
Sent: 03 March 2014 09:34
To: gmx-users at gromacs.org
Subject: [gmx-users] md_pull code in umbrella sampling

Dear Sir

I am trying to calculate the potential mean force (PMF) between protein
and Ligand. I have applied position restrain to protein. I have applied
reference pulling group to Protein and Pulling force has been applied to
ligand. The pul force is applied along Z-axis. I had query as to whether I
am proceeding in the right direction. Here I have provided the detail of
the content of md_pull.mdp

md_pull.mdp

title       = Umbrella pulling simulation
define      = -DPOSRES_Protein
; Run parameters
integrator  = md
dt          = 0.002
tinit       = 0
nsteps      = 250000    ; 500 ps
nstcomm     = 10
; Output parameters
nstxout     = 5000      ; every 10 ps
nstvout     = 5000
nstfout     = 500
nstxtcout   = 500       ; every 1 ps
nstenergy   = 500
; Bond parameters
constraint_algorithm    = lincs
constraints             = all-bonds
continuation            = yes       ; continuing from NPT
; Single-range cutoff scheme
nstlist     = 5
ns_type     = grid
rlist       = 1.4
rcoulomb    = 1.4
rvdw        = 1.4
; PME electrostatics parameters
coulombtype     = PME
fourierspacing  = 0.12
fourier_nx      = 0
fourier_ny      = 0
fourier_nz      = 0
pme_order       = 4
ewald_rtol      = 1e-5
optimize_fft    = yes
; Berendsen temperature coupling is on in two groups
Tcoupl      = Nose-Hoover
tc_grps     = Protein   Non-Protein
tau_t       = 0.5       0.5
ref_t       = 310       310
; Pressure coupling is on
Pcoupl          = Parrinello-Rahman
pcoupltype      = isotropic
tau_p           = 1.0
compressibility = 4.5e-5
ref_p           = 1.0
refcoord_scaling = com
; Generate velocities is off
gen_vel     = no
; Periodic boundary conditions are on in all directions
pbc     = xyz
; Long-range dispersion correction
DispCorr    = EnerPres
; Pull code
pull            = umbrella
pull_geometry   = distance  ; simple distance increase
pull_dim        = N N Y
pull_start      = yes       ; define initial COM distance > 0
pull_ngroups    = 1
pull_group0     = Protein
pull_group1     = DRG
pull_rate1      = 0.01      ; 0.01 nm per ps = 10 nm per ns
pull_k1         = 1000      ; kJ mol^-1 nm^-2

Regards


Arunima
--
Gromacs Users mailing list

* Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.