Merge pull request #1823 from swails/feature/flexibleConstraints

Flexible constraints

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

@@ -223,7 +223,7 @@ class ForceField(object):
             trees.append(tree)
 
         # Process includes.
-        
+
         for parentFile, tree in zip(files, trees):
             if isinstance(parentFile, str):
                 parentDir = os.path.dirname(parentFile)
@@ -1035,6 +1035,8 @@ class ForceField(object):
             not terminated properly.  This option can create ambiguities where multiple
             templates match the same residue.  If that happens, use the residueTemplates
             argument to specify which one to use.
+        flexibleConstraints : boolean=False
+            If True, parameters for constrained degrees of freedom will be added to the System
         args
             Arbitrary additional keyword arguments may also be specified.
             This allows extra parameters to be specified that are specific to
@@ -1140,9 +1142,9 @@ class ForceField(object):
             if not unit.is_quantity(hydrogenMass):
                 hydrogenMass *= unit.dalton
             for atom1, atom2 in topology.bonds():
-                if atom1.element == elem.hydrogen:
+                if atom1.element is elem.hydrogen:
                     (atom1, atom2) = (atom2, atom1)
-                if atom2.element == elem.hydrogen and atom1.element not in (elem.hydrogen, None):
+                if atom2.element is elem.hydrogen and atom1.element not in (elem.hydrogen, None):
                     transferMass = hydrogenMass-sys.getParticleMass(atom2.index)
                     sys.setParticleMass(atom2.index, hydrogenMass)
                     sys.setParticleMass(atom1.index, sys.getParticleMass(atom1.index)-transferMass)
@@ -1208,7 +1210,7 @@ class ForceField(object):
             for bond in data.bonds:
                 atom1 = data.atoms[bond.atom1]
                 atom2 = data.atoms[bond.atom2]
-                bond.isConstrained = atom1.name.startswith('H') or atom2.name.startswith('H')
+                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]
@@ -1224,11 +1226,11 @@ class ForceField(object):
                 atom2 = data.atoms[angle[1]]
                 atom3 = data.atoms[angle[2]]
                 numH = 0
-                if atom1.name.startswith('H'):
+                if atom1.element is elem.hydrogen:
                     numH += 1
-                if atom3.name.startswith('H'):
+                if atom3.element is elem.hydrogen:
                     numH += 1
-                data.isAngleConstrained.append(numH == 2 or (numH == 1 and atom2.name.startswith('O')))
+                data.isAngleConstrained.append(numH == 2 or (numH == 1 and atom2.element is elem.oxygen))
         else:
             data.isAngleConstrained = len(data.angles)*[False]
         if rigidWater:
@@ -1897,8 +1899,11 @@ class HarmonicBondGenerator(object):
                     bond.length = self.length[i]
                     if bond.isConstrained:
                         data.addConstraint(sys, bond.atom1, bond.atom2, self.length[i])
-                    elif self.k[i] != 0:
-                        force.addBond(bond.atom1, bond.atom2, self.length[i], self.k[i])
+                    if self.k[i] != 0:
+                        # flexibleConstraints allows us to add parameters even if the DOF is
+                        # constrained
+                        if not bond.isConstrained or args.get('flexibleConstraints', False):
+                            force.addBond(bond.atom1, bond.atom2, self.length[i], self.k[i])
                     break
 
 parsers["HarmonicBondForce"] = HarmonicBondGenerator.parseElement
@@ -1978,8 +1983,9 @@ class HarmonicAngleGenerator(object):
                             if l1 is not None and l2 is not None:
                                 length = sqrt(l1*l1 + l2*l2 - 2*l1*l2*cos(self.angle[i]))
                                 data.addConstraint(sys, angle[0], angle[2], length)
-                    elif self.k[i] != 0:
-                        force.addAngle(angle[0], angle[1], angle[2], self.angle[i], self.k[i])
+                    if self.k[i] != 0:
+                        if not isConstrained or args.get('flexibleConstraints', False):
+                            force.addAngle(angle[0], angle[1], angle[2], self.angle[i], self.k[i])
                     break
 
 parsers["HarmonicAngleForce"] = HarmonicAngleGenerator.parseElement
@@ -3160,8 +3166,9 @@ class AmoebaBondGenerator(object):
                     bond.length = self.length[i]
                     if bond.isConstrained:
                         data.addConstraint(sys, bond.atom1, bond.atom2, self.length[i])
-                    elif self.k[i] != 0:
-                        force.addBond(bond.atom1, bond.atom2, self.length[i], self.k[i])
+                    if self.k[i] != 0:
+                        if not bond.isConstrained or args.get('flexibleConstraints', False):
+                            force.addBond(bond.atom1, bond.atom2, self.length[i], self.k[i])
                     break
 
 parsers["AmoebaBondForce"] = AmoebaBondGenerator.parseElement
@@ -3288,6 +3295,8 @@ class AmoebaAngleGenerator(object):
         force.setAmoebaGlobalAnglePentic(self.pentic)
         force.setAmoebaGlobalAngleSextic(self.sextic)
 
+        DEG_TO_RAD = math.pi / 180
+
         for angleDict in angleList:
             angle = angleDict['angle']
             isConstrained = angleDict['isConstrained']
@@ -3302,8 +3311,8 @@ class AmoebaAngleGenerator(object):
                 if (type1 in types1 and type2 in types2 and type3 in types3) or (type1 in types3 and type2 in types2 and type3 in types1):
                     if isConstrained and self.k[i] != 0.0:
                         angleDict['idealAngle'] = self.angle[i][0]
-                        addAngleConstraint(angle, self.angle[i][0]*math.pi/180.0, data, sys)
-                    elif self.k[i] != 0:
+                        addAngleConstraint(angle, self.angle[i][0]*DEG_TO_RAD, data, sys)
+                    if self.k[i] != 0 and (not isConstrained or args.get('flexibleConstraints', False)):
                         lenAngle = len(self.angle[i])
                         if (lenAngle > 1):
                             # get k-index by counting number of non-angle hydrogens on the central atom
@@ -3372,7 +3381,7 @@ class AmoebaAngleGenerator(object):
                     angleDict['idealAngle'] = self.angle[i][0]
                     if (isConstrained and self.k[i] != 0.0):
                         addAngleConstraint(angle, self.angle[i][0]*math.pi/180.0, data, sys)
-                    else:
+                    if self.k[i] != 0.0 and (not isConstrained or args.get('flexibleConstraints', False)):
                         force.addAngle(angle[0], angle[1], angle[2], angle[3], self.angle[i][0], self.k[i])
                     break
 
@@ -5431,7 +5440,7 @@ class AmoebaUreyBradleyGenerator(object):
             force = existing[0]
 
         for (angle, isConstrained) in zip(data.angles, data.isAngleConstrained):
-            if (isConstrained):
+            if (isConstrained and not args.get('flexibleConstraints', False)):
                 continue
             type1 = data.atomType[data.atoms[angle[0]]]
             type2 = data.atomType[data.atoms[angle[1]]]

wrappers/python/tests/TestForceField.py

@@ -95,6 +95,37 @@ class TestForceField(unittest.TestCase):
                 validateConstraints(self, topology, system,
                                     constraints_value, rigidWater_value)
 
+    def test_flexibleConstraints(self):
+        """ Test the flexibleConstraints keyword """
+        topology = self.pdb1.topology
+        system1 = self.forcefield1.createSystem(topology, constraints=HAngles,
+                                                rigidWater=True)
+        system2 = self.forcefield1.createSystem(topology, constraints=HAngles,
+                                                rigidWater=True, flexibleConstraints=True)
+        system3 = self.forcefield1.createSystem(topology, constraints=None, rigidWater=False)
+        validateConstraints(self, topology, system1, HAngles, True)
+        # validateConstraints fails for system2 since by definition atom pairs can be in both bond
+        # and constraint lists. So just check that the number of constraints is the same for both
+        # system1 and system2
+        self.assertEqual(system1.getNumConstraints(), system2.getNumConstraints())
+        for force in system1.getForces():
+            if isinstance(force, HarmonicBondForce):
+                bf1 = force
+            elif isinstance(force, HarmonicAngleForce):
+                af1 = force
+        for force in system2.getForces():
+            if isinstance(force, HarmonicBondForce):
+                bf2 = force
+            elif isinstance(force, HarmonicAngleForce):
+                af2 = force
+        for force in system3.getForces():
+            if isinstance(force, HarmonicAngleForce):
+                af3 = force
+        # Make sure we picked up extra bond terms with flexibleConstraints
+        self.assertGreater(bf2.getNumBonds(), bf1.getNumBonds())
+        # Make sure flexibleConstraints yields just as many angles as no constraints
+        self.assertEqual(af2.getNumAngles(), af3.getNumAngles())
+
     def test_ImplicitSolvent(self):
         """Test the four types of implicit solvents using the implicitSolvent
         parameter.
@@ -687,7 +718,7 @@ class TestForceField(unittest.TestCase):
     def test_NBFix(self):
         """Test using LennardJonesGenerator to implement NBFix terms."""
         # Create a chain of five atoms.
-        
+
         top = Topology()
         chain = top.addChain()
         res = top.addResidue('RES', chain)
@@ -701,9 +732,9 @@ class TestForceField(unittest.TestCase):
         top.addBond(atoms[1], atoms[2])
         top.addBond(atoms[2], atoms[3])
         top.addBond(atoms[3], atoms[4])
-        
+
         # Create the force field and system.
-        
+
         xml = """
 <ForceField>
  <AtomTypes>
@@ -738,9 +769,9 @@ class TestForceField(unittest.TestCase):
 </ForceField> """
         ff = ForceField(StringIO(xml))
         system = ff.createSystem(top)
-        
+
         # Check that it produces the correct energy.
-        
+
         integrator = VerletIntegrator(0.001)
         context = Context(system, integrator, Platform.getPlatform(0))
         positions = [Vec3(i, 0, 0) for i in range(5)]*nanometers