Merge pull request #2529 from peastman/rigidwater

Improved behavior of rigidWater option

wrappers/python/simtk/openmm/app/data/amoeba2009.xml

@@ -2983,7 +2983,7 @@
    <ExternalBond from="14" />
    <ExternalBond from="0" />
  </Residue>
-  <Residue name="HOH">
+  <Residue name="HOH" rigidWater="false">
    <Atom name="H1" type="403" />
    <Atom name="H2" type="403" />
    <Atom name="O" type="402" />

wrappers/python/simtk/openmm/app/data/amoeba2013.xml

@@ -1756,7 +1756,7 @@
    <ExternalBond from="14" />
    <ExternalBond from="0" />
  </Residue>
-  <Residue name="HOH">
+  <Residue name="HOH" rigidWater="false">
    <Atom name="H1" type="248" />
    <Atom name="H2" type="248" />
    <Atom name="O" type="247" />

wrappers/python/simtk/openmm/app/data/iamoeba.xml

@@ -4,7 +4,7 @@
   <Type name="381" class="74" element="H" mass="1.008"/>
  </AtomTypes>
  <Residues>
-  <Residue name="HOH">
+  <Residue name="HOH" rigidWater="false">
    <Atom name="H1" type="381"/>
    <Atom name="H2" type="381"/>
    <Atom name="O" type="380"/>

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

@@ -264,6 +264,8 @@ class ForceField(object):
                     template = ForceField._TemplateData(resName)
                     if 'override' in residue.attrib:
                         template.overrideLevel = int(residue.attrib['override'])
+                    if 'rigidWater' in residue.attrib:
+                        template.rigidWater = (residue.attrib['rigidWater'].lower() == 'true')
                     atomIndices = template.atomIndices
                     for ia, atom in enumerate(residue.findall('Atom')):
                         params = {}
@@ -594,6 +596,7 @@ class ForceField(object):
             self.bonds = []
             self.externalBonds = []
             self.overrideLevel = 0
+            self.rigidWater = True
 
         def getAtomIndexByName(self, atom_name):
             """Look up an atom index by atom name, providing a helpful error message if not found."""
@@ -1051,7 +1054,7 @@ class ForceField(object):
         return [templates, unique_unmatched_residues]
 
     def createSystem(self, topology, nonbondedMethod=NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
-                     constraints=None, rigidWater=True, removeCMMotion=True, hydrogenMass=None, residueTemplates=dict(),
+                     constraints=None, rigidWater=None, removeCMMotion=True, hydrogenMass=None, residueTemplates=dict(),
                      ignoreExternalBonds=False, switchDistance=None, flexibleConstraints=False, **args):
         """Construct an OpenMM System representing a Topology with this force field.
 
@@ -1067,9 +1070,10 @@ class ForceField(object):
         constraints : object=None
             Specifies which bonds and angles should be implemented with constraints.
             Allowed values are None, HBonds, AllBonds, or HAngles.
-        rigidWater : boolean=True
+        rigidWater : boolean=None
             If true, water molecules will be fully rigid regardless of the value
-            passed for the constraints argument
+            passed for the constraints argument.  If None (the default), it uses the
+            default behavior for this force field's water model.
         removeCMMotion : boolean=True
             If true, a CMMotionRemover will be added to the System
         hydrogenMass : mass=None
@@ -1109,6 +1113,7 @@ class ForceField(object):
         data.atoms = list(topology.atoms())
         for atom in data.atoms:
             data.excludeAtomWith.append([])
+        rigidResidue = [False]*topology.getNumResidues()
 
         # Make a list of all bonds
 
@@ -1146,6 +1151,13 @@ class ForceField(object):
                     unmatchedResidues.append(res)
                 else:
                     data.recordMatchedAtomParameters(res, template, matches)
+                if res.name == 'HOH':
+                    # Determine whether this should be a rigid water.
+
+                    if rigidWater is None and template is not None:
+                        rigidResidue[res.index] = template.rigidWater
+                    elif rigidWater:
+                        rigidResidue[res.index] = True
 
         # Try to apply patches to find matches for any unmatched residues.
 
@@ -1269,11 +1281,10 @@ class ForceField(object):
                 atom1 = data.atoms[bond.atom1]
                 atom2 = data.atoms[bond.atom2]
                 bond.isConstrained = atom1.element is elem.hydrogen or atom2.element is elem.hydrogen
-        if rigidWater:
         for bond in data.bonds:
             atom1 = data.atoms[bond.atom1]
             atom2 = data.atoms[bond.atom2]
-                if atom1.residue.name == 'HOH' and atom2.residue.name == 'HOH':
+            if rigidResidue[atom1.residue.index] and rigidResidue[atom2.residue.index]:
                 bond.isConstrained = True
 
         # Identify angles that should be implemented with constraints
@@ -1291,13 +1302,12 @@ class ForceField(object):
                 data.isAngleConstrained.append(numH == 2 or (numH == 1 and atom2.element is elem.oxygen))
         else:
             data.isAngleConstrained = len(data.angles)*[False]
-        if rigidWater:
         for i in range(len(data.angles)):
             angle = data.angles[i]
             atom1 = data.atoms[angle[0]]
             atom2 = data.atoms[angle[1]]
             atom3 = data.atoms[angle[2]]
-                if atom1.residue.name == 'HOH' and atom2.residue.name == 'HOH' and atom3.residue.name == 'HOH':
+            if rigidResidue[atom1.residue.index] and rigidResidue[atom2.residue.index] and rigidResidue[atom3.residue.index]:
                 data.isAngleConstrained[i] = True
 
         # Add virtual sites

wrappers/python/tests/TestForceField.py

@@ -165,12 +165,12 @@ class TestForceField(unittest.TestCase):
 
         topology = self.pdb1.topology
         for constraints_value in [None, HBonds, AllBonds, HAngles]:
-            for rigidWater_value in [True, False]:
+            for rigidWater_value in [True, False, None]:
                 system = self.forcefield1.createSystem(topology,
                                                        constraints=constraints_value,
                                                        rigidWater=rigidWater_value)
                 validateConstraints(self, topology, system,
-                                    constraints_value, rigidWater_value)
+                                    constraints_value, rigidWater_value != False)
 
     def test_flexibleConstraints(self):
         """ Test the flexibleConstraints keyword """
@@ -1062,6 +1062,15 @@ class AmoebaTestForceField(unittest.TestCase):
         self.assertAlmostEqual(constraints[(0,2)], hoDist)
         self.assertAlmostEqual(constraints[(1,2)], hohDist)
         
+        # Check that all values of rigidWater are interpreted correctly.
+        
+        numWaters = 215
+        self.assertEqual(3*numWaters, system.getNumConstraints())
+        system = self.forcefield1.createSystem(self.pdb1.topology, rigidWater=False)
+        self.assertEqual(0, system.getNumConstraints())
+        system = self.forcefield1.createSystem(self.pdb1.topology, rigidWater=None)
+        self.assertEqual(0, system.getNumConstraints())
+
     def test_Forces(self):
         """Compute forces and compare them to ones generated with a previous version of OpenMM to ensure they haven't changed."""