addHydrogens() allows specifying exactly what hydrogens to add (#4585)

* addHydrogens() allows specifying exactly what hydrogens to add

* Prevent CI from using numpy 2.0

devtools/ci/gh-actions/conda-envs/build-M1-arm64.yml

@@ -10,7 +10,7 @@ dependencies:
 - python
 - cython
 - swig
-- numpy
+- numpy <2.0
 - doxygen 1.9.1
 # test
 - pytest

devtools/ci/gh-actions/conda-envs/build-macos-latest.yml

@@ -10,7 +10,7 @@ dependencies:
 - python
 - cython
 - swig
-- numpy
+- numpy <2.0
 - doxygen 1.8.14
 # test
 - pytest

devtools/ci/gh-actions/conda-envs/build-ubuntu-latest.yml

@@ -12,7 +12,7 @@ dependencies:
 - python
 - cython
 - swig
-- numpy
+- numpy <2.0
 - ocl-icd-system
 - doxygen 1.8.14
 # test

devtools/ci/gh-actions/conda-envs/build-ubuntu-pypy.yml

@@ -12,7 +12,7 @@ dependencies:
 - pypy
 - cython
 - swig
-- numpy
+- numpy <2.0
 - ocl-icd-system
 - doxygen 1.8.14
 # test

devtools/ci/gh-actions/conda-envs/build-windows-latest.yml

@@ -12,7 +12,7 @@ dependencies:
 - python
 - cython
 - swig
-- numpy
+- numpy <2.0
 - doxygen 1.8.14
 - khronos-opencl-icd-loader
 # test

devtools/ci/gh-actions/conda-envs/docs.yml

@@ -8,7 +8,7 @@ dependencies:
 # host
 - python
 - pip
-- numpy
+- numpy <2.0
 - cython
 - swig
 - doxygen 1.8.14

wrappers/python/openmm/app/modeller.py

@@ -736,6 +736,15 @@ class Modeller(object):
                     terminal = hydrogen.attrib['terminal']
                 data.hydrogens.append(Modeller._Hydrogen(hydrogen.attrib['name'], hydrogen.attrib['parent'], maxph, atomVariants, terminal))
 
+    @staticmethod
+    def _loadStandardHydrogenDefinitions():
+        """Load the definitions of hydrogens for standard residues.  Normally there is no need to call this directly.
+        It is automatically called by addHydrogens().  If the definitions have already been loaded, this returns without
+        doing anything."""
+        if not Modeller._hasLoadedStandardHydrogens:
+            Modeller.loadHydrogenDefinitions(os.path.join(os.path.dirname(__file__), 'data', 'hydrogens.xml'))
+            Modeller._hasLoadedStandardHydrogens = True
+
     def addHydrogens(self, forcefield=None, pH=7.0, variants=None, platform=None, residueTemplates=dict()):
         """Add missing hydrogens to the model.
 
@@ -798,7 +807,11 @@ class Modeller(object):
             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.
+            select a variant for that residue.  Alternatively, an element may specify
+            exactly which hydrogens to add.  In that case, variants[i] should be
+            a list of tuples [(name1, parent1), (name2, parent2), ...].  Each
+            tuple specifies the name of a hydrogen and the name of the parent atom
+            it should be bonded to.
         platform : Platform=None
             the Platform to use when computing the hydrogen atom positions.  If
             this is None, the default Platform will be used.
@@ -827,8 +840,7 @@ class Modeller(object):
 
         # Load the residue specifications.
 
-        if not Modeller._hasLoadedStandardHydrogens:
-            Modeller.loadHydrogenDefinitions(os.path.join(os.path.dirname(__file__), 'data', 'hydrogens.xml'))
+        Modeller._loadStandardHydrogenDefinitions()
 
         # Make a list of atoms bonded to each atom.
 
@@ -867,10 +879,9 @@ class Modeller(object):
                     newResidueTemplates[newResidue] = residueTemplates[residue]
                 isNTerminal = (residue == chain._residues[0])
                 isCTerminal = (residue == chain._residues[-1])
-                if residue.name in Modeller._residueHydrogens:
+                if residue.name in Modeller._residueHydrogens or isinstance(variants[residue.index], list):
                     # Add hydrogens.  First select which variant to use.
 
-                    spec = Modeller._residueHydrogens[residue.name]
                     variant = variants[residue.index]
                     if variant is None:
                         if residue.name == 'CYS':
@@ -930,8 +941,14 @@ class Modeller(object):
                                     variant = 'HID'
                         elif residue.name == 'HIS':
                             variant = 'HIP'
-                    if variant is not None and variant not in spec.variants:
-                        raise ValueError('Illegal variant for %s residue: %s' % (residue.name, variant))
+                    if isinstance(variant, list):
+                        spec = Modeller._ResidueData(residue.name)
+                        infinity = float('Inf')
+                        spec.hydrogens = [Modeller._Hydrogen(name, parent, infinity, None, None) for name, parent in variant]
+                    else:
+                        spec = Modeller._residueHydrogens[residue.name]
+                        if variant is not None and variant not in spec.variants:
+                            raise ValueError('Illegal variant for %s residue: %s' % (residue.name, variant))
                     actualVariants[residue.index] = variant
                     removeExtraHydrogens = (variants[residue.index] is not None)
 

wrappers/python/tests/TestModeller.py

 from collections import defaultdict
 import unittest
+import random
 
 from validateModeller import *
 from openmm.app import *
@@ -980,6 +981,33 @@ class TestModeller(unittest.TestCase):
             names1 = sorted([a.name for a in res1.atoms()])
             names2 = sorted([a.name for a in res2.atoms()])
             self.assertEqual(names1, names2)
+        # Reset the loaded definitions so we don't affect other tests.
+        Modeller._residueHydrogens = {}
+        Modeller._hasLoadedStandardHydrogens = False
+
+    def test_addSpecificHydrogens(self):
+        """Test specifying exactly which hydrogens to add."""
+        pdb = PDBFile('systems/glycopeptide.pdb')
+        variants = [None]*pdb.topology.getNumResidues()
+        for residue in pdb.topology.residues():
+            if residue.name != 'ALA':
+                var = []
+                for atom1, atom2 in residue.bonds():
+                    if atom1.element == element.hydrogen:
+                        var.append((atom1.name, atom2.name))
+                    elif atom2.element == element.hydrogen:
+                        var.append((atom2.name, atom1.name))
+                variants[residue.index] = var
+        modeller = Modeller(pdb.topology, pdb.positions)
+        hydrogens = [a for a in modeller.topology.atoms() if a.element == element.hydrogen and random.random() < 0.7]
+        modeller.delete(hydrogens)
+        self.assertTrue(modeller.topology.getNumAtoms() < pdb.topology.getNumAtoms())
+        modeller.addHydrogens(variants=variants)
+        self.assertEqual(modeller.topology.getNumAtoms(), pdb.topology.getNumAtoms())
+        for res1, res2 in zip(pdb.topology.residues(), modeller.topology.residues()):
+            names1 = sorted([a.name for a in res1.atoms()])
+            names2 = sorted([a.name for a in res2.atoms()])
+            self.assertEqual(names1, names2)
 
     def test_removeExtraHydrogens(self):
         """Test that addHydrogens() can remove hydrogens that shouldn't be there. """