openmm / openmm

OpenMM is a toolkit for molecular simulation using high performance GPU code.
1.5k stars 523 forks source link

System creation failure for deprotonated lysine #462

Closed leeping closed 3 years ago

leeping commented 10 years ago

OpenMM fails to build the System for some unconventional amino acids.

pdb = PDBFile(sys.argv[1])
ff = ForceField('amber99sb.xml', 'tip3p.xml')
system = ff.createSystem(pdb.topology)

Traceback (most recent call last):
  File "../HIP/min.py", line 31, in <module>
    system = ff.createSystem(pdb.topology)
  File "/home/leeping/local/lib/python2.7/site-packages/simtk/openmm/app/forcefield.py", line 333, in createSystem
    raise ValueError('No template found for residue %d (%s).  %s' % (res.index+1, res.name, _findMatchErrors(self, res)))
ValueError: No template found for residue 1 (ACE).  The set of atoms matches ACE, but the bonds are different.  Perhaps the chain is missing a terminal group?

The PDB in question is here: https://www.dropbox.com/s/j75p39mrplwhe4p/ace-LYN-nme.pdb

The same error occurs for protonated aspartic acid (ASH) and glutamic acid (GLH), as well as deprotonated cysteine (CYM).

leeping commented 10 years ago

The topology building works if I explicitly specify the CONECT records in the pdb.

jchodera commented 9 years ago

It seems like this might still need some sort of fix in PDBFile, such as the addition of standard bond perception for deprotonated lysine. Maybe this can be marked bug?

Astrovicis commented 4 years ago

This is happening quite often for me. Eagerly awaiting a solution 😉

jchodera commented 4 years ago

@Astrovicis : Which force field(s) are you attempting to use, out of curiosity?

jchodera commented 4 years ago

It looks like amber14/protein.ff14SB.xml supports deprotonated lysine (LYN), but the residues.xml lacks a definition for this amino acid variant.

@peastman: What do you think about adding these protonation state variants to residues.xml? The PDB standard says that CONECT records are "only mandatory for HET groups (excluding water) and for other bonds not specified in the standard residue connectivity table". Since the definition of HET residues is vague, it could be interpreted as suggesting that standard amino acid variants / protonation states should be considered as not requiring CONECT residues.

@Astrovicis : As a workaround, you can add CONECT records to your PDB file.

peastman commented 4 years ago

What name does the residue have in the PDB file? If it has the standard name LYS, the PDB loader should handle it correctly even if one atom is missing. If the force field then has an appropriate template, it should work. On the other hand, if you've given it a nonstandard name, there's no way the PDB loader can know what it's supposed to be.

Since the definition of HET residues is vague, it could be interpreted as suggesting that standard amino acid variants / protonation states should be considered as not requiring CONECT residues.

The definition is pretty clear. The chemical component dictionary defines exactly what atoms are expected to be present in each standard residue. Missing atoms are tolerated, but the addition of any atom not present in the standard residue definition is not allowed.

jchodera commented 4 years ago

The definition is pretty clear. The chemical component dictionary defines exactly what atoms are expected to be present in each standard residue. Missing atoms are tolerated, but the addition of any atom not present in the standard residue definition is not allowed.

That's not quite correct. The chemical component dictionary for HIS, for example, includes atom names for both delta and epsilon protons and the C-terminal carboxylic acid hydrogen, but it does not include a third amino terminal proton needed if HIS were to appear at the N-terminus of a protein. And yet, despite this, we include this third amino-terminal proton in residues.xml.

jchodera commented 4 years ago

But @peastman has a good point: Replacing nonstandard amino acid names (LYN) with standard ones (LYS) in your PDB file will allow examples like this to work!