Modeller.addHydrogens() allows the user to provide an XML file with additional residue definitions

wrappers/python/simtk/openmm/app/modeller.py

@@ -487,7 +487,30 @@ class Modeller(object):
             self.variants = variants
             self.terminal = terminal
     
-    def addHydrogens(self, forcefield, pH=7.0, variants=None):
+    def _loadHydrogenDefinitions(self, file):
+        residueHydrogens = {}
+        tree = etree.parse(file)
+        infinity = float('Inf')
+        for residue in tree.getroot().findall('Residue'):
+            resName = residue.attrib['name']
+            data = Modeller._ResidueData(resName)
+            residueHydrogens[resName] = data
+            for variant in residue.findall('Variant'):
+                data.variants.append(variant.attrib['name'])
+            for hydrogen in residue.findall('H'):
+                maxph = infinity
+                if 'maxph' in hydrogen.attrib:
+                    maxph = float(hydrogen.attrib['maxph'])
+                atomVariants = None
+                if 'variant' in hydrogen.attrib:
+                    atomVariants = hydrogen.attrib['variant'].split(',')
+                terminal = None
+                if 'terminal' in hydrogen.attrib:
+                    terminal = hydrogen.attrib['terminal']
+                data.hydrogens.append(Modeller._Hydrogen(hydrogen.attrib['name'], hydrogen.attrib['parent'], maxph, atomVariants, terminal))
+        return residueHydrogens
+    
+    def addHydrogens(self, forcefield, pH=7.0, variants=None, hydrogenDefinitions=None):
         """Add missing hydrogens to the model.
         
         Some residues can exist in multiple forms depending on the pH and properties of the local environment.  These
@@ -531,6 +554,9 @@ class Modeller(object):
          - variants (list=None) an optional list of variants to use.  If this is specified, its length must equal the number
            of residues in the model.  variants[i] is the name of the variant to use for residue i (indexed starting at 0).
            If an element is None, the standard rules will be followed to select a variant for that residue.
+         - hydrogenDefinitions (File=None) an optional XML file containing definitions for hydrogens to add.  The built in
+           hydrogens.xml file contains definitions for standard amino acids and nucleotides.  This parameter can be used to
+           provide additional definitions for other residue types.
         Returns: a list of what variant was actually selected for each residue, in the same format as the variants parameter
         """
         # Check the list of variants.
@@ -546,26 +572,12 @@ class Modeller(object):
         # Load the residue specifications.
         
         if Modeller._residueHydrogens is None:
-            Modeller._residueHydrogens = {}
-            tree = etree.parse(os.path.join(os.path.dirname(__file__), 'data', 'hydrogens.xml'))
-            infinity = float('Inf')
-            for residue in tree.getroot().findall('Residue'):
-                resName = residue.attrib['name']
-                data = Modeller._ResidueData(resName)
-                Modeller._residueHydrogens[resName] = data
-                for variant in residue.findall('Variant'):
-                    data.variants.append(variant.attrib['name'])
-                for hydrogen in residue.findall('H'):
-                    maxph = infinity
-                    if 'maxph' in hydrogen.attrib:
-                        maxph = float(hydrogen.attrib['maxph'])
-                    atomVariants = None
-                    if 'variant' in hydrogen.attrib:
-                        atomVariants = hydrogen.attrib['variant'].split(',')
-                    terminal = None
-                    if 'terminal' in hydrogen.attrib:
-                        terminal = hydrogen.attrib['terminal']
-                    data.hydrogens.append(Modeller._Hydrogen(hydrogen.attrib['name'], hydrogen.attrib['parent'], maxph, atomVariants, terminal))
+            Modeller._residueHydrogens = self._loadHydrogenDefinitions(os.path.join(os.path.dirname(__file__), 'data', 'hydrogens.xml'))
+        if hydrogenDefinitions is None:
+            residueHydrogens = Modeller._residueHydrogens
+        else:
+            residueHydrogens = Modeller._residueHydrogens.copy()
+            residueHydrogens.update(self._loadHydrogenDefinitions(hydrogenDefinitions))
 
         # Make a list of atoms bonded to each atom.
         
@@ -601,10 +613,10 @@ class Modeller(object):
                 newResidue = newTopology.addResidue(residue.name, newChain)
                 isNTerminal = (residue == chain._residues[0])
                 isCTerminal = (residue == chain._residues[-1])
-                if residue.name in Modeller._residueHydrogens:
+                if residue.name in residueHydrogens:
                     # Add hydrogens.  First select which variant to use.
                     
-                    spec = Modeller._residueHydrogens[residue.name]
+                    spec = residueHydrogens[residue.name]
                     variant = variants[residue.index]
                     if variant is None:
                         if residue.name == 'CYS':