GROMACS version: 2025.3
GROMACS modification: Yes (patched with PLUMED 2.9.4)
Hello, I am new to this forum and hope that this is the correct place to post my question!
I want to simulate a protein in a certain protonation state in GROMACS. I have protonated the protein with:
pdb2pqr --ff=AMBER --titration-state-method=propka --with-ph=6.8.
Now when I try to initialize the parametrization with:
gmx_mpi pdb2gmx -o processed.gro -p topol.top -ff amber99sb-star-ildnp -water tip3p
I get the following error:
Fatal error:
Atom HB3 in residue MET 1 was not found in rtp entry NMET with 19 atoms
while sorting atoms.
For a hydrogen, this can be a different protonation state, or it
might have had a different number in the PDB file and was rebuilt
(it might for instance have been H3, and we only expected H1 & H2).
Note that hydrogens might have been added to the entry for the N-terminus.
Remove this hydrogen or choose a different protonation state to solve it.
Option -ignh will ignore all hydrogens in the input.
My problem: using -ignh does not make any sense to me because I want to simulate exactly this protonation state. How can I tell GROMACS that the third H on that residue has to be there (and similarly for all other added H’s)? Or what is the default way to deal with this?
Using -ignh is the appropriate solution. pdb2gmx will rebuild the correct H atoms with correct names. Your issue is that your Cβ has atoms HB2 and HB3 bonded to it (very common) but GROMACS expects HB1 and HB2. Rather than renaming everything (error-prone and tedious), let pdb2gmx rebuild the appropriate H atoms for you. It won’t change protonation state.
Thank you for your answer! I ran the paramtrization with -ignh and compared the two files (the protonated pdb and the gro), and the processed gro file has 1 hydrogen less than the original pdb that I protonated at pH6.8. I also ran the parametrization with the original, unprotonated file, which resulted in the same amount of hydrogens as the protonated gro. This means that -ignh does affect the protonation state and will undo the effects of running pdb2pqr, right (even if that apparently is just one H in this case)? I also tried it with H++ Is there a way to keep the protonation state and simulate the protonated protein? Why do pdb2pqr and H++ even generate these weirdly named files that GROMACS cannot read?
pdb2gmx assumes a pH of 7 and that all residues have canonical protonation states. If your species is protonated differently, then you need to individually assign protonation states with the corresponding command-line options that allow you to tune these for every titratable residue.
Welcome to the Wild West of incompatible formats. Wait until you get into non-standard species, where most programs call them “UNK” and don’t even write distinct atom names, just elements…
Thank you. I solved this problem by protonating my structure at a pH of 7 and a pH of 6.8, then using a python script to count the protons per residue and rename every titratable residue with a differing number of protons between the two in the original pdb, which I then parametrized and let GROMACS build the protons.. Not the cleanest solution but at least it worked.