How to pull an ion towards binding site in Umbrella sampling

User discussions
1

GROMACS version: 2022.4
GROMACS modification: No
Hi,
I’m struggling with a very naive problem may be. My system is an ion channel in a bilayer membrane. Now I’m trying pull an ion along positive Z-axis (aligned with pore) from the solvent to the pore using COM of pore.
My basic understanding so far says we always pull from a bound state (like ligand in an active site) and away from it. But is it possible to pull things inside a cleft or pore too using Umbrella sampling method ?
I’m trying to work with a very basic pull code which looks like this:

title       = Umbrella pulling simulation
;
define      = -DPOSRES -DPOSRES_FC_BB=50.0 -DPOSRES_FC_SC=0.0 -DPOSRES_FC_LIPID=0.0 -DDIHRES -DDIHRES_FC=0.0
; 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
cutoff-scheme   = Verlet
nstlist         = 20
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       = 1.0	1.0
ref_t       = 277.15    277.15
; 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                    = yes
pull_ncoords            = 1         ; only one reaction coordinate
pull_ngroups            = 2         ; two groups defining one reaction coordinate
pull_group1_name        = Pore
pull_group2_name        = X
pull_coord1_type        = umbrella  ; harmonic potential
pull_coord1_geometry    = distance  ; simple distance increase
pull_coord1_dim         = N N Y
pull_coord1_groups	= 1 2
pull_coord1_start	= yes       ; define initial COM distance > 0
pull_coord1_rate        = 0.01      ; 0.01 nm per ps = 10 nm per ns
pull_coord1_k           = 1000      ; kJ mol^-1 nm^-2

Where X is the ion I’m trying to pull and Pore is the channel-pore where I want the ion to end up after the pull.
When I grompp this (using:gmx grompp -f md_pull.mdp -c step7_production_volt_10.gro -p topol.top -r step7_production_volt_10.gro -n index.ndx -t step7_production_volt_10.cpt -o pull.tpr) I’m getting the following error:

ERROR 1 [file md_pull.mdp]:
  When the maximum distance from a pull group reference atom to other atoms
  in the group is larger than 0.5 times half the box size a centrally
  placed atom should be chosen as pbcatom. Pull group 1 is larger than that
  and does not have a specific atom selected as reference atom.

Pull group  natoms  pbc atom  distance at start  reference at t=0
       1       197     14834
       2         1         0       2.273 nm          2.273 nm

Please suggest how this can be achieved to generate a PMF.

Thank you.

2

Dear Abhysek,

This means that group 1 is very large and some atoms in that group are distant from each other more than 0.5 times the box size. This is a problem for GROMACS when it tries to calculate distances, since PBC are on and the true distance is not well defined anymore.

Try to visualize group 1 with some tool like VMD and try to find an atom that is ca. at the centre of your group 1. Take that atom number from the .gro/.pdb/your structure file and then add the option

pull-group1-pbcatom = #number_of_your_central_atom

to your mdp file. Take a look at the mdp pull options in GROMACS manual here.

When you do not specify this option, I think that GROMACS just searches for the central atom in the group from the list, i.e., if you have a group of 10 atoms gromacs uses the fifth, which does not correspond necessarily with the geometrical centre of your group. This usually is not a problem, but when groups are larger than half the box then the choice could be not ideal to calculate distances with PBC, and so you have to specify a better atom for GROMACS to use.

Also, I think it may be relevant for you to consider adding also a cylindrical constrain to keep the ion within the channel (not necessarily, but just a hint).

Hope this helps