Merge pull request #505 from peastman/forcefield

Added a test case to make sure future changes to ForceField do not change forces

wrappers/python/tests/TestForceField.py

@@ -132,6 +132,23 @@ class TestForceField(unittest.TestCase):
         totalMass1 = sum([system1.getParticleMass(i) for i in range(system1.getNumParticles())]).value_in_unit(amu)
         totalMass2 = sum([system2.getParticleMass(i) for i in range(system2.getNumParticles())]).value_in_unit(amu)
         self.assertAlmostEqual(totalMass1, totalMass2)
+    
+    def test_Forces(self):
+        """Compute forces and compare them to ones generated with a previous version of OpenMM to ensure they haven't changed."""
+        
+        pdb = PDBFile('systems/lysozyme-implicit.pdb')
+        system = self.forcefield2.createSystem(pdb.topology)
+        integrator = VerletIntegrator(0.001)
+        context = Context(system, integrator)
+        context.setPositions(pdb.positions)
+        state1 = context.getState(getForces=True)
+        state2 = XmlSerializer.deserialize(open('systems/lysozyme-implicit-forces.xml').read())
+        numDifferences = 0
+        for f1, f2, in zip(state1.getForces().value_in_unit(kilojoules_per_mole/nanometer), state2.getForces().value_in_unit(kilojoules_per_mole/nanometer)):
+            diff = norm(f1-f2)
+            if diff > 0.1 and diff/norm(f1) > 1e-3:
+                numDifferences += 1
+        self.assertTrue(numDifferences < system.getNumParticles()/20) # Tolerate occasional differences from numerical error
 
 class AmoebaTestForceField(unittest.TestCase):
     """Test the ForceField.createSystem() method with the AMOEBA forcefield."""