Don't add extra bond between water hydrogens (#4219)

wrappers/python/openmm/app/amberprmtopfile.py

@@ -160,7 +160,13 @@ class AmberPrmtopFile(object):
 
         atoms = list(top.atoms())
         for bond in prmtop.getBondsWithH():
-            top.addBond(atoms[bond[0]], atoms[bond[1]])
+            a1 = atoms[bond[0]]
+            a2 = atoms[bond[1]]
+            if a1.residue.name == 'HOH' and a1.element == elem.hydrogen and a2.element == elem.hydrogen:
+                # Don't add the "bond" Amber lists between the two hydrogens of a water molecule.
+                pass
+            else:
+                top.addBond(a1, a2)
         for bond in prmtop.getBondsNoH():
             top.addBond(atoms[bond[0]], atoms[bond[1]])
 

wrappers/python/tests/TestAmberPrmtopFile.py

@@ -78,12 +78,14 @@ class TestAmberPrmtopFile(unittest.TestCase):
         """Test all eight options for the constraints and rigidWater parameters."""
 
         topology = prmtop1.topology
-        for constraints_value in [None, HBonds, AllBonds, HAngles]:
-            for rigidWater_value in [True, False]:
-                system = prmtop1.createSystem(constraints=constraints_value,
-                                              rigidWater=rigidWater_value)
-                validateConstraints(self, topology, system,
-                                    constraints_value, rigidWater_value)
+        for constraints in [None, HBonds, AllBonds, HAngles]:
+            for rigidWater in [True, False]:
+                system = prmtop1.createSystem(constraints=constraints, rigidWater=rigidWater)
+                if constraints != None:
+                    # Amber adds an extra "bond" between water hydrogens, so any constraint
+                    # method except None is equivalent to rigidWater=True.
+                    rigidWater = True
+                validateConstraints(self, topology, system, constraints, rigidWater)
 
     def test_ImplicitSolvent(self):
         """Test the four types of implicit solvents using the implicitSolvent
@@ -460,5 +462,16 @@ class TestAmberPrmtopFile(unittest.TestCase):
             # If a constraint was added to a massless particle, this will throw an exception.
             context = Context(system, integrator, Platform.getPlatformByName('Reference'))
 
+    def testWaterBonds(self):
+        """Test that water molecules have the right set of bonds"""
+        top = prmtop1.topology
+        for residue in top.residues():
+            if residue.name == 'HOH':
+                bonds = list(residue.bonds())
+                self.assertEqual(2, len(bonds))
+                for a1, a2 in bonds:
+                    self.assertTrue(a1.element == elem.oxygen or a2.element == elem.oxygen)
+                    self.assertTrue(a1.element == elem.hydrogen or a2.element == elem.hydrogen)
+
 if __name__ == '__main__':
     unittest.main()

wrappers/python/tests/TestForceField.py

@@ -439,7 +439,6 @@ class TestForceField(unittest.TestCase):
    <Atom name="H2" type="tip3p-H" charge="0.417"/>
    <Bond from="0" to="1"/>
    <Bond from="0" to="2"/>
-   <Bond from="1" to="2"/>
   </Residue>
   <Residue name="HOH2">
    <Atom name="O" type="tip3p-O" charge="0.834"/>
@@ -447,7 +446,6 @@ class TestForceField(unittest.TestCase):
    <Atom name="H2" type="tip3p-H" charge="-0.417"/>
    <Bond from="0" to="1"/>
    <Bond from="0" to="2"/>
-   <Bond from="1" to="2"/>
   </Residue>
  </Residues>
  <NonbondedForce coulomb14scale="0.833333" lj14scale="0.5">