Fix #3790: Collision rate for MTS BAOAB Langevin integrator now correctly...

Fix #3790: Collision rate for MTS BAOAB Langevin integrator now correctly accounts for number of substeps (#3791)

* Fix #3790: Collision rate for MTS BAOAB Langevin integrator now correctly accounts for number of substeps. Fix suggested by Charlie Matthews (@c-matthews)

* Fix computation of total number of substeps

* Add test for MTS friction

* Fix typo

* Fix yet another typo

* Fix typo and check against analytical result

* Fix typo

* Fix typos

* Fix more typos

* Integrate MTSLangevinIntegrator for longer to allow thermalization

* Revert number of integrator steps

* Update TestIntegrators.py

Fixed a failing test case
Co-authored-by: Peter Eastman <peter.eastman@gmail.com>

wrappers/python/openmm/mtsintegrator.py

@@ -167,8 +167,9 @@ class MTSLangevinIntegrator(CustomIntegrator):
         self.temperature = temperature
         self.friction = friction
         import math
-        self.addGlobalVariable("a", math.exp(-friction*dt))
-        self.addGlobalVariable("b", math.sqrt(1-math.exp(-2*friction*dt)))
+        total_substeps = groups[-1][1]
+        self.addGlobalVariable("a", math.exp(-friction*dt / total_substeps))
+        self.addGlobalVariable("b", math.sqrt(1-math.exp(-2*friction*dt / total_substeps)))
         from openmm.unit import MOLAR_GAS_CONSTANT_R
         self.addGlobalVariable('kT', MOLAR_GAS_CONSTANT_R*temperature)
         self.addPerDofVariable("x1", 0)

wrappers/python/tests/TestIntegrators.py

@@ -137,7 +137,7 @@ class TestIntegrators(unittest.TestCase):
                 force.setForceGroup(0)
 
         # Create an integrator
-        integrator = MTSLangevinIntegrator(300*kelvin, 1/picosecond, 4*femtoseconds, [(2,1), (1,2), (0,4)])
+        integrator = MTSLangevinIntegrator(300*kelvin, 5/picosecond, 4*femtoseconds, [(2,1), (1,2), (0,4)])
 
         # Run some equilibration.
         context = Context(system, integrator)
@@ -156,6 +156,36 @@ class TestIntegrators(unittest.TestCase):
         temperature = averageEnergy*2/(dofs*MOLAR_GAS_CONSTANT_R)
         self.assertTrue(290*kelvin < temperature < 310*kelvin)
 
+    def testMTSLangevinIntegratorFriction(self):
+        """Test the MTSLangevinIntegrator on a force-free particle to ensure friction is properly accounted for (issue #3790)"""
+        # Create a System with a single particle and no forces
+        system = System()
+        system.addParticle(12.0*amu)
+        platform = Platform.getPlatformByName('Reference')
+        initial_positions = [Vec3(0,0,0)]
+        initial_velocities = [Vec3(1,0,0)]
+        nsteps = 125 # number of steps to take
+        collision_rate = 1/picosecond
+        timestep = 4*femtoseconds
+
+        def get_final_velocities(nsubsteps):
+            """Get the final velocity vector after a fixed number of steps for the specified number of substeps"""
+            integrator = MTSLangevinIntegrator(0*kelvin, collision_rate, timestep, [(0,nsubsteps)])
+            context = Context(system, integrator, platform)
+            context.setPositions(initial_positions)
+            context.setVelocities(initial_velocities)
+            integrator.step(nsteps)
+            final_velocities = context.getState(getVelocities=True).getVelocities()
+            del context, integrator
+            return final_velocities
+
+        # Compare sub-stepped MTS with single-step MTS
+        for nsubsteps in range(2,6):
+            mts_velocities = get_final_velocities(nsubsteps)
+            self.assertAlmostEqual(math.exp(-timestep*nsteps*collision_rate), mts_velocities[0].x)
+            self.assertAlmostEqual(0, mts_velocities[0].y)
+            self.assertAlmostEqual(0, mts_velocities[0].z)
+
     def testNoseHooverIntegrator(self):
         """Test partial thermostating in the NoseHooverIntegrator (only API)"""
         pdb = PDBFile('systems/alanine-dipeptide-explicit.pdb')