Protonation/disulfide issues with pdb2gmx?

GROMACS version: 2021.5
GROMACS modification: No
Here post your question

Dear GROMACS users,
I would like to ask for some advice on file preparation and pdb2gmx. I plan to run simulations in higher than neutral pH, so I had to determine the correct protonation state of my protein. I used pdb2pqr/propka for this, which returns the structure in Amber format (having LYN, HID, CYX, etc. residues). GROMACS uses a different format for hydrogens, so when I process the file with pdb2gmx I have to use the -ignh flag, but GROMACS adds the correct number of hydrogens to all residues.

My input contains also cysteines forming disulfide bridges (CYX) and they are correctly protonated by GROMACS, both the *.gro and topology files contain the right number of hydrogens, and I get no error messages about having CYX in the structure. However, despite the correct protonation (indicating the presence of a disulfide bridge) when CYX is in the input, GROMACS does not recognize disulfide bridges and the simulation is performed as they weren’t there, resulting in erratic dynamics.

So I have two questions: First, I wonder whether there is a way of making GROMACS “CYX-aware”, or is it a known issue/bug, and cysteines need to be corrected manually? Second, I wonder whether in the case of other residues of my input (LYN, HID), which are also correctly protonated by GROMACS, are their hydrogen bonds treated correctly during simulations?

Thank you.

PS.: I am using Amber (ildn) force field.

Dear @ga01

To fix this problem, try to follow these steps:

1- Go to where the GROMACS installation folder is.
2- Find the “specbond.dat” file and open it.
3- Add the following information to the specbond.dat file:

CYX SG 1 CYX SG 1 0.2 CYS2 CYS2

This modification will fix your issue.

Best regards,

Hassan