Bug #2480

pdb2gmx does not protonate correctly for united-atom Gromos

Added by Mikko Huhtala 3 months ago. Updated about 1 month ago.

preprocessing (pdb2gmx,grompp)
Target version:
Affected version - extra info:
Affected version:


pdb2gmx builds ASPH residues without the HD2 hydrogen for Gromos 54a7. I changed manually the residue name of Asp95 in the PDB file 3GUU from ASP to ASPH and ran pdb2gmx. The result was a ASPH residue name both in the topology and the structure file, but both without the HD2 hydrogen, i.e. ASP structure, not ASPH. Changing His366 to HISA worked correctly, giving a correct HISA residue in the output. Edited PDB file included (3GUU A chain + crystal water, excluding heteroatoms, modified Asp95 and His366).

Version: Gromacs 2018, single precision, Linux x86_64.

Discussed on the mailing list:

3guu_A.pdb (336 KB) 3guu_A.pdb Mikko Huhtala, 04/16/2018 06:14 PM


#1 Updated by Mikko Huhtala 2 months ago

Tested 2018.1, the issue is still there.

#2 Updated by Mikko Huhtala 2 months ago

Realized that the issue cuts both ways: renaming ASP to ASPH in the input will build a ASP topology, and having an actual ASPH with the hydrogen in the input is not recognized. I get the following error:

"Fatal error:
Atom HD2 in residue ASPH 61 was not found in rtp entry ASP with 9 atoms
while sorting atoms."

ASPH has 10 atoms, not 9. For some reason pdb2gmx tries to build an ASP residue with 9 atoms.

This is when selecting GROMOS96 54a7 in Gromacs 2018.1.

#3 Updated by Mark Abraham 2 months ago

What happens when you use gmx pdb2gmx -asp, which is the intended way for you to specify protonation states of such residues?

#4 Updated by Mikko Huhtala 2 months ago

That will work, but avoiding that was my whole point starting the thread on the mailing list. The manual entering of protonation states is a pain, if you need to e.g. input protonation calculated with another tool, or prepare multiple versions of the protein with residue substitutions etc. You want text input that's readable and scriptable.

#5 Updated by Mikko Huhtala about 1 month ago

Work intervened, so I'm adding this comment late. I made another post weeks ago on the mailing list, and Justin Lemkul confirmed that pdb2gmx does behave inconsistently with different force fields. ASH gets built correctly for Amber and Charmm, but ASPH does not for Gromos and OPLS-AA. See

Justin Lemkul commented earlier that building the protonated states according to the force field residue name in the input PDB file is the expected behavior. I can't actually find an explicit statement of that anywhere in Gromacs documentation, but that may be just me. If the user is always expected to enter the protonation states manually, I think there should be at least a mention of it somewhere ("no guarantees on what happens to ASH/ASPH etc." or something to that effect). In any case, I think the current interactive interface is a very bad idea anyway - entering long lists of residues by hand. If you make one mistake, you have to start from the beginning. And the residues are not in the order of sequence, but by type, so using the output of a external protonation calculation tool manually becomes very inconvenient. E.g. my case happened to be a 430-residue protein with 6 non-default protonation states across three different residue types (total > 40 residues). Preparing topologies for a wild type and a couple of mutants means going through the error-prone manual keying three times.

#6 Updated by Mark Abraham about 1 month ago

  • Target version set to 2018.2

We wholeheartedly agree that pdb2gmx needs an overhaul. I'll see what I can do on the minor points for 2018.2, but the structural issues will need more far-reaching re-design of system preparation workflows.

#7 Updated by Mark Abraham about 1 month ago

  • Category set to preprocessing (pdb2gmx,grompp)
  • Assignee set to Mark Abraham
  • Target version changed from 2018.2 to 2018.3

Will reconsider for 2018.3

Also available in: Atom PDF