GROMACS version: 2022.4
GROMACS modification: No
When altering the protonation state of a positively charged ARG to a neutral ARGN using topologies from pdb2gmx with the -arg flag enabled with the charmm36 FF, there seems to be massive mdrun instability the simulation exploding. grompp has no complaints and the instability is centered around the ARGN atoms HH11 and HH21.
The system I am testing on is a coordinate file from PDB:1UBQ with no alterations and waters removed. I am running in a 999nm^3 box to simulate vacuum conditions (see Konermann, L. et al., Methods 144, 104–112 (2018) if curious) although I have tried running with a smaller box size/with water to no avail. I have also listed all the files used in my test below. This failure happens both on native windows and linux (HPC) builds with/without GPU acceleration.
Through testing as far as I can tell this is not specific to the coordinate file I am using (PDB:1UBQ), it is not specific to any one ARG residue (I have tested each ARG residue in 1ubq.pdb), and it is not specific to the charmm36 version I am using (I have tested charmm36-jul2017.ff, charmm36-mar2019.ff, and charmm36-jul2022.ff). All other protonable amino acids including LYS, GLU, ASP, HIS, and C/N-termini can handle altered protonation states fine so it seems to be specific to ARG. Altering ARG → ARGN works fine with other FF’s like OPLS or GROMOS, so it also seem to be specific to CHARMM36.
My respective mdp files and gmx calls are also listed below. When run, the NVT equilibriation always fails around step 600-800. The trajectory file is quite interesting as the deprotonated ARG behaves normally until HH11 and HH12 pull together seemingly randomly at step ~600 before oscillating to explosion. If you run the prodrun.mdp file first you get similar results. I have tried editing nearly every parameter of these two files to no help. If run with only protonated ARG the calls below will run to completion with no complaints.
gmx calls
#set box size to 10 nm^3 and remove waters
gmx editconf -f 1ubq.pdb -o prot.gro -box 10.0
echo q | gmx make_ndx -f prot.gro -o 1ubq.ndx
echo 1 | gmx editconf -f pro.gro -o prot.gro -n 1ubq.ndx
#gmx pdb2gmx call, selecting 1 (protonated ARG form) for all ARG will run fine
#if select 0 for deprotonated ARG no matter which one, the run will fail once it reaches mdrun
gmx pdb2gmx -f prot.gro -o prot.gro -p topol.top -ff charmm36 \
-water none -ter -lys -arg -asp -gly -his -ignh
#em equilibriation will run fine even with ARGN
echo q | gmx make_ndx -f prot.gro -o prot.ndx
gmx grompp -f em.mdp -c prot.gro -r prot.gro -p topol.top -n prot.ndx -o em.tpr -maxwarn 10
gmx mdrun -v -deffnm em
#running protein in vacuum so need to expand to max box size
#see Konermann, L. et al., Methods 144, 104–112 (2018) if curious
gmx editconf -f em.gro -o expanded.gro -box 999
#nvt equilibriation, will fail here at step ~600 with ARGN
echo q | gmx make_ndx -f expanded.gro -o prot.ndx
gmx grompp -f nvt.mdp -c expanded.gro -r expanded.gro -p topol.top -n prot.ndx -o nvt.tpr -maxwarn 10
gmx mdrun -v -deffnm nvt
#npt equilibriation
echo q | gmx make_ndx -f nvt.gro -o prot.ndx
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -p topol.top -n prot.ndx -o npt.tpr -maxwarn 10
gmx mdrun -v -deffnm npt
#prodrun, will fail similarly to nvt if run first
echo q | gmx make_ndx -f npt.gro -o prot.ndx
gmx grompp -f prodrun.mdp -c npt.gro -r npt.gro -p topol.top -n prot.ndx -o prodrun.tpr -maxwarn 10
gmx mdrun -v -deffnm prodrun
em.mdp
integrator = steep
emtol = 0.01
emstep = 0.01
nstxout-compressed = 100
nsteps = 10000
nstlist = 10
cutoff-scheme = Verlet
rlist = 1.2
coulombtype = cut-off
rcoulomb = 1.2
vdwtype = cut-off
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2
pbc = xyz
DispCorr = no
nvt.mdp file
integrator = md
nsteps = 10000
dt = 0.001
nstxout = 0
nstvout = 0
nstenergy = 0
nstlog = 0
nstxout-compressed = 1
continuation = no
constraint_algorithm = lincs
constraints = h-bonds
lincs_iter = 1
lincs_order = 4
comm-mode = linear
cutoff-scheme = Verlet
nstlist = 20
rlist = 333
vdwtype = cut-off
rvdw = 333
DispCorr = no
coulombtype = cut-off
rcoulomb = 333
tcoupl = nose-Hoover
nh-chain-length = 1
tc-grps = system
tau_t = 2.0
ref_t = 370
pcoupl = no
gen_vel = yes
pbc = xyz
gen_temp = 370
gen_seed = -1
npt.mdp file
integrator = md
nsteps = 10000
dt = 0.001
nstxout = 0
nstvout = 0
nstenergy = 0
nstlog = 0
nstxout-compressed = 0
continuation = yes
constraint_algorithm = lincs
constraints = h-bonds
lincs_iter = 1
lincs_order = 4
comm-mode = linear
cutoff-scheme = Verlet
nstlist = 20
rlist = 333
vdwtype = cut-off
rvdw = 333
DispCorr = no
coulombtype = cut-off
rcoulomb = 333
tcoupl = Nose-Hoover
nh-chain-length = 1
tc-grps = System
tau_t = 2.0
ref_t = 370
pcoupl = Parrinello-Rahman
pcoupltype = isotropic
tau_p = 4.0
ref_p = 0.0
compressibility = 4.5e-5
refcoord_scaling = com
pbc = xyz
gen_vel = no
prodrun.mdp
integrator = md
nsteps = 10000
dt = 0.001
nstenergy = 0
nstcalcenergy = 0
nstlog = 0
nstxout-compressed = 0
continuation = yes
constraint_algorithm = lincs
constraints = h-bonds
lincs_iter = 1
lincs_order = 4
comm-mode = linear
cutoff-scheme = Verlet
nstlist = 20
rlist = 333
vdwtype = cut-off
rvdw = 333
DispCorr = no
coulombtype = cut-off
rcoulomb = 333
tcoupl = Nose-Hoover
nh-chain-length = 1
tc-grps = System
tau_t = 2.0
ref_t = 370
pcoupl = no
pbc = xyz
gen_vel = no
Any help would be greatly appreciated.
Thanks,
Michael Cordes
Baylor University