Segmentation fault during slow SMD pulling through a constricted pore

[Priyanka_yadav]

GROMACS version: 2023.3
GROMACS modification: Yes/No

Hi everyone,

I am performing a steered MD (SMD) simulation in GROMACS to pull a Na⁺ ion through a protein channel.

I use the pore-lining Cα atoms as the reference group, pull-coord1-geometry = cylinder, and pull-cylinder-r = 0.6 nm.

I can successfully pull the ion over a distance of 12 nm in 200 ps using a relatively fast pulling rate. However, I wanted to perform a smoother pulling simulation, so I reduced the pulling speed to 2 Å/ns (0.0002 nm/ps), which requires about 60 ns to pull the ion the same distance.

With the slower pulling rate, the simulation runs initially, but when the ion reaches a constricted region of the pore, it eventually terminates with a segmentation fault. There are no LINCS warnings or other obvious GROMACS fatal errors before the crash.

These are all parameters

integrator = md
dt = 0.002
nsteps = 30000000
nstxout-compressed = 5000
nstxout = 0
nstvout = 0
nstfout = 0
nstcalcenergy = 100
nstenergy = 100
nstlog = 100
;
cutoff-scheme = Verlet
nstlist = 20
rlist = 1.2
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 1.0
rvdw = 1.2
coulombtype = PME
rcoulomb = 1.2
;
tcoupl = v-rescale
tc_grps = SOLU MEMB SOLV
tau_t = 1.0 1.0 1.0
ref_t = 303.15 303.15 303.15
;
pcoupl = C-rescale
pcoupltype = semiisotropic
tau_p = 5.0
compressibility = 4.5e-5 4.5e-5
ref_p = 1.0 1.0
;
constraints = h-bonds
constraint_algorithm = LINCS
continuation = yes
;
nstcomm = 100
comm_mode = linear
comm_grps = SOLU_MEMB SOLV

;SMD

pull = yes
pull-ngroups = 2
pull-ncoords = 1

pull-group1-name = PORE ; Ca atoms of pore lining residues
pull-group2-name = NA_PULL : serial 175028
pull-coord1-groups = 1 2

pull-coord1-type = umbrella
pull-coord1-geometry = cylinder
pull-cylinder-r = 0.6

pull-coord1-dim = N N Y
pull-coord1-vec = 0 0 -1
pull-coord1-start = yes ; Define initial COM distance >0
pull-print-com = yes
pull-print-ref-value = yes
pull-coord1-rate = 0.0002 ; nm/ps = 0.2 nm/ns
pull-coord1-k = 1000 ; kJ mol^-1 nm^-2
pull-pbc-ref-prev-step-com = yes

pull-nstxout = 100
pull-nstfout = 100
~

I have also tried 500 kJ mol⁻¹ nm⁻², still the same issue.

My questions are:

• Is using the pore-lining Cα atoms as the reference group with pull-coord1-geometry = cylinder and pull-cylinder-r = 0.6 nm appropriate for SMD?
• Could the very slow pulling rate (0.0002 nm/ps) be causing the instability, or is the segmentation fault likely due to another issue?
• Since my goal is only to generate umbrella sampling windows, would it be better to use a faster pulling rate (e.g., complete the 12 nm pull in 500 ps–1 ns)?

Any suggestions would be greatly appreciated.

[hess]

You should never get an segmentation fault without any warnings before that.

How much does the location of the ion deviate from the reference close the time of the crash?

I am wondering if your rather low force constant could be the source of the issue.

[Priyanka_yadav]

Dear Hess,
Thank you for your response.

• The simulation consistently terminates with Segmentation fault (core dumped) whenever the Na⁺ ion encounters a significant barrier while entering the channel. This can occur even before reaching the selectivity filter, for example due to interactions with the extracellular loop region.
• I restarted the simulation 2–3 times, but each time it crashes within a few seconds after the ion encounters such a barrier and does not proceed further.
• I selected the Na⁺ ion that was initially located at the extracellular (top) mouth of the channel and pulled it along the pore axis toward the intracellular side.
• The ion enters the pore but, upon interacting with the loop region or the selectivity filter, it begins to fluctuate around that region instead of progressing further. Shortly afterward, the simulation terminates with a segmentation fault.
• I am attaching the pulling distance (pullx.xvg) and pulling force (pullf.xvg) plots near the time of the crash.
• I am also attaching snapshots of the system. The Na⁺ ion is shown in yellow, the selectivity filter residues are shown in red, and the pore-lining Cα atoms used as the reference group are shown in blue. The snapshots correspond to the initial configuration and the configuration near the crash (~15 ns).

[hess]

I see nothing out of the normal there.

I now realize that you are using geometry=cylinder with a very small radius. Maybe the causes issues. Normally one would use a region of the protein as reference group with geometry=direction.

[Priyanka_yadav]

Dear Hess,

Thank you for your suggestions.

Based on your recommendations, I modified the pulling setup as follows:

• pull-coord1-geometry = direction. I also defined pull-group1-pbcatom to ensure that the pulling reference remained consistent and the Na⁺ ion followed the channel axis instead of diffusing away from the pore.
• pull-group1-pbcatom = 15612 (atom near the center of the protein used for PBC reconstruction).
• pull-coord1-k = 1000 kJ mol⁻¹ nm⁻².
• pull-pbc-ref-prev-step-com = yes.

I tested the simulation for 10 ns, and so far it is running smoothly. The Na⁺ ion is able to pass through the channel without the segmentation fault that I was encountering previously.

Thank you very much for your time and helpful suggestions.

[Priyanka_yadav]

Dear Hess,

I have a follow-up question regarding the same SMD setup. The simulation worked as expected for a shorter run (10 ns). However, when I extended the simulation to 120 ns and reduced the pulling rate to:

pull-coord1-rate = 0.0001 ; nm/ps
pull-coord1-k    = 1000   ; kJ mol^-1 nm^-2
pull-pbc-ref-prev-step-com  = yes

During the slow pulling simulation, the Na⁺ ion diffuses away from the pore entrance into the bulk solvent, causing the pulling restraint to drag the ion outside the channel instead of through the pore.