Preserve bond type and order in all `Modeller` operations (#4323)

* Py wrapper: preserve bond data in Modeller operations

* regression tests for add(), delete(), deleteWater(), addSolvent()

* Tests for addHydrogens and addExtraParticles

* Cosmetic

* A single regression test spanning a whole workflow

* Remove now-redundant tests

* Test also deleteWater and addHydrogens. Remove excessive assertions.

wrappers/python/openmm/app/modeller.py

@@ -178,7 +178,7 @@ class Modeller(object):
         for bond in self.topology.bonds():
             if bond[0] in newAtoms and bond[1] in newAtoms:
                 if bond not in deleteSet and (bond[1], bond[0]) not in deleteSet:
-                    newTopology.addBond(newAtoms[bond[0]], newAtoms[bond[1]])
+                    newTopology.addBond(newAtoms[bond[0]], newAtoms[bond[1]], bond.type, bond.order)
         self.topology = newTopology
         self.positions = newPositions
 
@@ -253,7 +253,7 @@ class Modeller(object):
                         newPositions.append(deepcopy(self.positions[atom.index]))
         for bond in self.topology.bonds():
             if bond[0] in newAtoms and bond[1] in newAtoms:
-                newTopology.addBond(newAtoms[bond[0]], newAtoms[bond[1]])
+                newTopology.addBond(newAtoms[bond[0]], newAtoms[bond[1]], bond.type, bond.order)
         self.topology = newTopology
         self.positions = newPositions
 
@@ -549,7 +549,7 @@ class Modeller(object):
                     newAtoms[atom] = newAtom
                     newPositions.append(deepcopy(self.positions[atom.index]))
         for bond in self.topology.bonds():
-            newTopology.addBond(newAtoms[bond[0]], newAtoms[bond[1]])
+            newTopology.addBond(newAtoms[bond[0]], newAtoms[bond[1]], bond.type, bond.order)
 
         # Sort the solute atoms into cells for fast lookup.
 
@@ -1000,7 +1000,7 @@ class Modeller(object):
                         newPositions.append(deepcopy(self.positions[atom.index]))
         for bond in self.topology.bonds():
             if bond[0] in newAtoms and bond[1] in newAtoms:
-                newTopology.addBond(newAtoms[bond[0]], newAtoms[bond[1]])
+                newTopology.addBond(newAtoms[bond[0]], newAtoms[bond[1]], bond.type, bond.order)
 
         # The hydrogens were added at random positions.  Now perform an energy minimization to fix them up.
 
@@ -1224,7 +1224,7 @@ class Modeller(object):
 
         for bond in self.topology.bonds():
             if bond[0] in newAtoms and bond[1] in newAtoms:
-                newTopology.addBond(newAtoms[bond[0]], newAtoms[bond[1]])
+                newTopology.addBond(newAtoms[bond[0]], newAtoms[bond[1]], bond.type, bond.order)
 
         if len(missingPositions) > 0:
             # There were particles whose position we couldn't identify before, since they were neither virtual sites nor Drude particles.

wrappers/python/tests/TestModeller.py

@@ -1200,6 +1200,75 @@ class TestModeller(unittest.TestCase):
         for i in range(3):
             self.assertTrue(newSize[i] >= originalSize[i]+1.1*nanometers)
 
+    def test_bondTypeAndOrderPreserved(self):
+        """ Check that bond type and order are preserved across multiple operations. 
+
+        Regression test for issue #4112 and similar behaviors.
+        """
+
+        # Given: an isolated molecule
+        pdb = PDBFile("systems/alanine-dipeptide-implicit.pdb")
+        topology, positions = pdb.topology, pdb.positions
+        topology.setUnitCellDimensions(Vec3(3.5, 3.5, 3.5) * nanometers)
+        # with some bonds carrying type and order information
+        for bond in topology.bonds():
+            if ((bond.atom1.element, bond.atom2.element) in [
+                (element.carbon, element.oxygen), (element.oxygen, element.carbon)
+            ]):
+                bond.type = Double
+                bond.order = 2.0
+        modeller = Modeller(topology, positions)
+
+        # When (1): add solvent
+        forcefield = ForceField("amber10.xml", "tip3p.xml")
+        modeller.addSolvent(forcefield, "tip3p")
+        # sanity check: water was added
+        self.assertTrue(any(r.name == "HOH" for r in modeller.topology.residues()))
+
+        # When (2): convert water (no sites added)
+        modeller.convertWater("spce")
+
+        # When (3): convert water with addExtraParticles
+        new_forcefield = ForceField('amber10.xml', 'tip4pew.xml')
+        modeller.addExtraParticles(new_forcefield)
+        # sanity check: extra sites were added
+        self.assertEqual(
+            set([len(list(res.atoms())) for res in modeller.topology.residues() if res.name == "HOH"]),
+            {4}
+        )
+
+        # When (4): delete water (with deleteWater) and hydrogens (with delete)
+        modeller.deleteWater()
+        hydrogens = [a for a in modeller.topology.atoms() if a.element == element.hydrogen]
+        modeller.delete(hydrogens)
+        # sanity check: all gone
+        self.assertFalse(any(a.element == element.hydrogen for a in modeller.topology.atoms()))
+        self.assertFalse(any(r.name == "HOH" for r in modeller.topology.residues()))
+
+        # When (5): add back hydrogens
+        modeller.addHydrogens()
+        # sanity check: hydrogens are back
+        self.assertTrue(any(a.element == element.hydrogen for a in modeller.topology.atoms()))
+
+        # Then (intermediate): bond info have been retained throughout the workflow
+        self.assertIn((Double, 2.0), [(b.type, b.order) for b in modeller.topology.bonds()])
+
+        # When (6): add a modeller (which also bears some bond info)
+        to_add = PDBFile('systems/methanol-box.pdb')
+        topology_to_add = to_add.topology
+        positions_to_add = to_add.positions
+        # add a dummy bond to the "to_add" system to check that it also is preserved
+        atom0, atom1 = (atom for i, atom in enumerate(topology_to_add.atoms()) if i < 2)
+        topology_to_add.addBond(atom0, atom1, Single, 1.0)
+        modeller.add(topology_to_add, positions_to_add)
+
+        # Then: bond info are retained for both the old and the new system
+        all_bond_extra_data = [(b.type, b.order) for b in modeller.topology.bonds()]
+        self.assertEqual(
+            set(all_bond_extra_data),
+            # None and Double from topology; other Nones and Single from topology_to_add
+            {(None, None), (Single, 1.0), (Double, 2.0)}
+        )
 
     def assertVecAlmostEqual(self, p1, p2, tol=1e-7):
         scale = max(1.0, norm(p1),)