Merge pull request #1618 from jchodera/fix-mts

Fix precision issue in MTSIntegrator

wrappers/python/simtk/openmm/mtsintegrator.py

@@ -98,10 +98,11 @@ class MTSIntegrator(CustomIntegrator):
         for i in range(stepsPerParentStep):
             self.addComputePerDof("v", "v+0.5*(dt/"+str(substeps)+")*f"+str(group)+"/m")
             if len(groups) == 1:
-                self.addComputePerDof("x1", "x")
                 self.addComputePerDof("x", "x+(dt/"+str(substeps)+")*v")
+                self.addComputePerDof("x1", "x")
                 self.addConstrainPositions();
-                self.addComputePerDof("v", "(x-x1)/(dt/"+str(substeps)+")");
+                self.addComputePerDof("v", "v+(x-x1)/(dt/"+str(substeps)+")");
+                self.addConstrainVelocities()
             else:
                 self._createSubsteps(substeps, groups[1:])
             self.addComputePerDof("v", "v+0.5*(dt/"+str(substeps)+")*f"+str(group)+"/m")

wrappers/python/tests/TestIntegrators.py

+import unittest
+import tempfile
+from datetime import datetime, timedelta
+from simtk.openmm import *
+from simtk.openmm.app import *
+from simtk.unit import *
+import math, random
+
+class TestIntegrators(unittest.TestCase):
+    """Test Python Integrator classes"""
+
+    def testMTSIntegratorExplicit(self):
+        """Test the MTS integrator on an explicit solvent system"""
+        # Create a periodic solvated system with PME
+        pdb = PDBFile('systems/alanine-dipeptide-explicit.pdb')
+        ff = ForceField('amber99sbildn.xml', 'tip3p.xml')
+        system = ff.createSystem(pdb.topology, cutoffMethod=PME)
+
+        # Split forces into groups
+        for force in system.getForces():
+            if force.__class__.__name__ == 'NonbondedForce':
+                force.setForceGroup(1)
+                force.setReciprocalSpaceForceGroup(2)
+            else:
+                force.setForceGroup(0)
+
+        # Create an integrator
+        integrator = MTSIntegrator(4*femtoseconds, [(2,1), (1,2), (0,8)])
+
+        # Run a few steps of dynamics
+        context = Context(system, integrator)
+        context.setPositions(pdb.positions)
+        integrator.step(10)
+
+        # Ensure energy is well-behaved.
+        state = context.getState(getEnergy=True)
+        if not (state.getPotentialEnergy() / kilojoules_per_mole < 0.0):
+            raise Exception('Potential energy of alanine dipeptide system with MTS integrator is blowing up: %s' % str(state.getPotentialEnergy()))
+
+    def testMTSIntegratorConstraints(self):
+        """Test the MTS integrator energy conservation on a system of constrained particles with no inner force (just constraints)"""
+
+        # Create a constrained test system
+        numParticles = 8
+        numConstraints = 5
+        system = System()
+        force = NonbondedForce()
+        for i in range(numParticles):
+            system.addParticle(5.0 if i%2==0 else 10.0)
+            force.addParticle((0.2 if i%2==0 else -0.2), 0.5, 5.0);
+        system.addConstraint(0, 1, 1.0);
+        system.addConstraint(1, 2, 1.0);
+        system.addConstraint(2, 3, 1.0);
+        system.addConstraint(4, 5, 1.0);
+        system.addConstraint(6, 7, 1.0);
+        system.addForce(force)
+
+        # Create integrator where inner timestep just evaluates constraints
+        integrator = MTSIntegrator(1*femtoseconds, [(1,1), (0,4)])
+        integrator.setConstraintTolerance(1e-5);
+
+        positions = [ (i/2., (i+1)/2., 0.) for i in range(numParticles) ]
+        velocities = [ (random.random()-0.5, random.random()-0.5, random.random()-0.5) for i in range(numParticles) ]
+
+        # Create Context
+        platform = Platform.getPlatformByName('Reference')
+        context = Context(system, integrator, platform)
+        context.setPositions(positions)
+        context.setVelocities(velocities)
+        context.applyConstraints(1e-5)
+
+        # Simulate it and see whether the constraints remain satisfied.
+        CONSTRAINT_RELATIVE_TOLERANCE = 1.e-4 # relative constraint violation tolerance
+        ENERGY_RELATIVE_TOLERANCE = 1.e-2 # relative energy violation tolerance
+        for i in range(1000):
+            state = context.getState(getPositions=True, getEnergy=True)
+            positions = state.getPositions()
+            for j in range(numConstraints):
+                [particle1, particle2, constraint_distance] = system.getConstraintParameters(j)
+                current_distance = 0.0 * nanometers**2
+                for k in range(3):
+                    current_distance += (positions[particle1][k] - positions[particle2][k])**2
+                current_distance = sqrt(current_distance)
+                # Fail test if outside of relative tolerance
+                relative_violation = (current_distance - constraint_distance) / constraint_distance
+                if (relative_violation > CONSTRAINT_RELATIVE_TOLERANCE):
+                    raise Exception('Constrained distance is violated by relative tolerance of %f (constraint %s actual %s)' % (relative_violation, str(constraint_distance), str(current_distance)))
+            # Check total energy
+            total_energy = state.getPotentialEnergy() + state.getKineticEnergy()
+            if (i == 1):
+                initial_energy = total_energy
+            elif (i > 1):
+                relative_violation = abs((total_energy - initial_energy) / initial_energy)
+                if (relative_violation > ENERGY_RELATIVE_TOLERANCE):
+                    raise Exception('Total energy is violated by relative tolerance of %f on step %d (initial %s final %s)' % (relative_violation, i, str(initial_energy), str(total_energy)))
+            # Take a step
+            integrator.step(1)
+
+if __name__ == '__main__':
+    unittest.main()