Add hydrogens to caps (#320)

pdbfixer/pdbfixer.py

@@ -389,11 +389,13 @@ class PDBFixer(object):
         # Load the templates.
         # Load the templates.
 
 
         self.templates = {}
         self.templates = {}
+        self._standardTemplates = set()
         templatesPath = os.path.join(os.path.dirname(__file__), 'templates')
         templatesPath = os.path.join(os.path.dirname(__file__), 'templates')
         for file in os.listdir(templatesPath):
         for file in os.listdir(templatesPath):
             templatePdb = app.PDBFile(os.path.join(templatesPath, file))
             templatePdb = app.PDBFile(os.path.join(templatesPath, file))
             name = next(templatePdb.topology.residues()).name
             name = next(templatePdb.topology.residues()).name
             self.templates[name] = Template(templatePdb.topology, templatePdb.positions)
             self.templates[name] = Template(templatePdb.topology, templatePdb.positions)
+            self._standardTemplates.add(name)
 
 
     def _initializeFromPDB(self, file):
     def _initializeFromPDB(self, file):
         """Initialize this object by reading a PDB file."""
         """Initialize this object by reading a PDB file."""
@@ -1413,7 +1415,7 @@ class PDBFixer(object):
 
 
     def _describeVariant(self, residue, definitions):
     def _describeVariant(self, residue, definitions):
         """Build the variant description to pass to addHydrogens() for a residue."""
         """Build the variant description to pass to addHydrogens() for a residue."""
-        if residue.name not in app.PDBFile._standardResidues and self._getTemplate(residue.name) is not None:
+        if residue.name not in self._standardTemplates and self._getTemplate(residue.name) is not None:
             # The user has registered a template for this residue.  Use the hydrogens from it.
             # The user has registered a template for this residue.  Use the hydrogens from it.
             template = self._getTemplate(residue.name)
             template = self._getTemplate(residue.name)
             atoms = [(atom.name, atom.element.symbol.upper(), terminal) for atom, terminal in zip(template.topology.atoms(), template.terminal)]
             atoms = [(atom.name, atom.element.symbol.upper(), terminal) for atom, terminal in zip(template.topology.atoms(), template.terminal)]

pdbfixer/tests/data/alanine-dipeptide.pdb

+REMARK   1 CREATED WITH OPENMM 8.2, 2024-12-03
+HETATM    1  CH3 ACE A   1       2.000   2.090   0.000  1.00  0.00           C  
+HETATM    2  C   ACE A   1       3.427   2.641  -0.000  1.00  0.00           C  
+HETATM    3  O   ACE A   1       4.391   1.877  -0.000  1.00  0.00           O  
+ATOM      4  N   ALA A   2       3.555   3.970  -0.000  1.00  0.00           N  
+ATOM      5  CA  ALA A   2       4.853   4.614  -0.000  1.00  0.00           C  
+ATOM      6  CB  ALA A   2       5.661   4.221  -1.232  1.00  0.00           C  
+ATOM      7  C   ALA A   2       4.713   6.129   0.000  1.00  0.00           C  
+ATOM      8  O   ALA A   2       3.601   6.653   0.000  1.00  0.00           O  
+HETATM    9  N   NME A   3       5.846   6.835   0.000  1.00  0.00           N  
+HETATM   10  C   NME A   3       5.846   8.284   0.000  1.00  0.00           C  
+TER      11      NME A   3
+CONECT    1    2
+CONECT    2    1    3    4
+CONECT    3    2
+CONECT    4    2
+CONECT    7    9
+CONECT    9    7   10
+CONECT   10    9
+END

pdbfixer/tests/test_add_hydrogens.py

@@ -42,3 +42,10 @@ def test_registered_template():
         count = sum(1 for atom in residue.atoms() if atom.element.symbol == 'H')
         count = sum(1 for atom in residue.atoms() if atom.element.symbol == 'H')
         if residue.name == 'CAS':
         if residue.name == 'CAS':
             assert count == 9
             assert count == 9
+
+def test_end_caps():
+    """Test adding hydrogens to a chain capped with ACE and NME."""
+    fixer = pdbfixer.PDBFixer(filename=(Path(__file__).parent / "data" / "alanine-dipeptide.pdb").as_posix())
+    fixer.addMissingHydrogens()
+    forcefield = app.ForceField('amber14/protein.ff14SB.xml')
+    forcefield.createSystem(fixer.topology)
\ No newline at end of file