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 ending with the dreaded LINCS or XTC warning as the system explodes. 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), but have also run with smaller box sizes 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. I am guessing it is some improper parameters within the charmm36.ff files, which I am currently looking into, but hesitant to alter given the complexity and effect it may have on proper parametrization.
Altering protonation states with pdb2gmx manually is quite tedious so I built a python script, protonation_test.py copied below to automate the process as well as automate mdrun for em/NVT/NPT equilibration by using the subprocess module. The .py file runs on windows, but might not work properly on mac/linux without some alterations or you can just run the gmx calls normally without using the .py file. If you do run the .py file, start it in a dir with 1ubq.pdb and the respective mdp files which 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.
protonation_test.py file with deprotonated 42ARG (see comments in code below)
#!/usr/bin/env python
import subprocess
#expands coordinate file to 10.0 nm^3 box and removes waters
subprocess.run(['gmx', 'editconf', '-f', '1ubq.pdb', '-o', 'prot.gro', '-box', '10.0'], shell=True)
subprocess.run(['echo', 'q', '|', 'gmx', 'make_ndx', '-f', 'prot.gro', '-o', '1ubq.ndx'], shell=True)
subprocess.run(['echo', '1', '|', 'gmx', 'editconf', '-f', 'prot.gro', '-o', 'prot.gro', '-n', '1ubq.ndx'], shell=True)
#list of ubiq amino acid residues from PDB:1UBQ
proton_res = ['6LYSN', '11LYSN', '27LYS', '29LYS', '33LYS', '48LYS', \
'63LYS', '42ARG', '54ARG', '72ARG', '74ARG', '21ASP', '32ASP', '39ASP', \
'52ASP', '58ASP', '16GLU', '18GLU', '24GLU', '34GLU', '51GLU', '64GLU', \
'68HISE', 'Nterm', 'Cterm']
#dict of protonation states to be fed into pdb2gmx
#flip any of the ARG numbers from 1 (protonated ARG) -> 0 (deprotonated ARGN)
#currently 42ARG is set to be deprotonated as it is set to 0
#if 42ARG set to 1, the simulation will run with no errors to completion
proton_map = {'6LYSN': 1, '11LYSN': 1, '27LYS': 1, '29LYS': 1, '33LYS': 1, \
'48LYS': 1, '63LYS': 1, '42ARG': 0, '54ARG': 1, '72ARG': 1, '74ARG': 1, '21ASP': 0, \
'32ASP': 0, '39ASP': 0, '52ASP': 0, '58ASP': 0, '16GLU': 0, '18GLU': 0, '24GLU': 0, \
'34GLU': 0, '51GLU': 0, '64GLU': 0, '68HISE': 1, 'Nterm': 1, 'Cterm': 0}
#calling pdb2gmx and feeding values from dict
proton_map_vals = []
for i in range(len(proton_res)):
proton_map_vals.append(str(proton_map[proton_res[i]]))
proc = subprocess.Popen(['gmx','pdb2gmx' ,'-f', 'prot.gro', '-o', 'prot.gro', '-p', 'topol.top', '-ff', 'charmm36', \
'-water', 'none', '-ter', '-lys', '-arg', '-asp', '-glu', '-his', '-ignh'], stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
for input in proton_map_vals:
proc.stdin.write(f"{input}\n".encode())
proc.stdin.flush()
output, errors = proc.communicate()
print(output)
#em equilibriation
subprocess.run(['echo','q','|','gmx','make_ndx','-f','prot.gro','-o','prot.ndx'], shell=True)
subprocess.run(['gmx','grompp','-f','em.mdp','-c','prot.gro','-r','prot.gro','-p','topol.top','-n','prot.ndx','-o','em.tpr','-maxwarn','10'], shell=True)
subprocess.run(['gmx','mdrun','-v','-deffnm','em'], shell=True)
#running protein in vacuum so need to expand to max box size
#see Konermann, L. et al., Methods 144, 104–112 (2018) if curious
subprocess.run(['gmx', 'editconf', '-f', 'em.gro', '-o', 'expanded.gro', '-box', '999'], shell=True)
#nvt equilibriation
subprocess.run(['echo','q','|','gmx','make_ndx','-f','expanded.gro','-o','prot.ndx'], shell=True)
subprocess.run(['gmx','grompp','-f','nvt.mdp','-c','expanded.gro','-r','expanded.gro','-p','topol.top','-n','prot.ndx','-o','nvt.tpr','-maxwarn','10'], shell=True)
subprocess.run(['gmx','mdrun','-v','-deffnm','nvt'], shell=True)
#npt equilibriation
subprocess.run(['echo','q','|','gmx','make_ndx','-f','nvt.gro','-o','prot.ndx'], shell=True)
subprocess.run(['gmx','grompp','-f','npt.mdp','-c','nvt.gro','-r','nvt.gro','-p','topol.top','-n','prot.ndx','-o','npt.tpr','-maxwarn','10'], shell=True)
subprocess.run(['gmx','mdrun','-v','-deffnm','npt'], shell=True)
#prodrun
subprocess.run(['echo','q','|','gmx','make_ndx','-f','npt.gro','-o','prot.ndx'], shell=True)
subprocess.run(['gmx','grompp','-f','prodrun.mdp','-c','npt.gro','-r','npt.gro','-p','topol.top','-n','prot.ndx','-o','prodrun.tpr','-maxwarn','10'], shell=True)
subprocess.run(['gmx','mdrun','-v','-deffnm','prodrun'], shell=True)
#deletes tpr/log files so mdruns don't rerun old files if grompp fails
#this is with windows syntax, change del to rm if on mac/linux
subprocess.run(['del', '*log'], shell=True)
subprocess.run(['del', '*tpr'], shell=True)
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