Merge pull request #2695 from peastman/mtslangevin

Created MTSLangevinIntegrator

docs-source/usersguide/application.rst

@@ -1123,6 +1123,13 @@ algorithm\ :cite:`Tuckerman1992`.  This allows some forces in the system to be e
 frequently than others.  For details on how to use it, consult the API
 documentation.
 
+Multiple Time Step Langevin Integrator
+--------------------------------------
+
+:class:`MTSLangevinIntegrator` is similar to :class:`MTSIntegrator`, but it uses
+the Langevin method to perform constant temperature dynamics.  For details on
+how to use it, consult the API documentation.
+
 Compound Integrator
 -------------------
 

wrappers/python/simtk/openmm/__init__.py

@@ -18,7 +18,7 @@ if sys.platform == 'win32':
 
 from simtk.openmm.openmm import *
 from simtk.openmm.vec3 import Vec3
-from simtk.openmm.mtsintegrator import MTSIntegrator
+from simtk.openmm.mtsintegrator import MTSIntegrator, MTSLangevinIntegrator
 from simtk.openmm.amd import AMDIntegrator, AMDForceGroupIntegrator, DualAMDIntegrator
 
 if os.getenv('OPENMM_PLUGIN_DIR') is None and os.path.isdir(version.openmm_library_path):

wrappers/python/simtk/openmm/mtsintegrator.py

@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
 Biological Structures at Stanford, funded under the NIH Roadmap for
 Medical Research, grant U54 GM072970. See https://simtk.org.
 
-Portions copyright (c) 2013-2015 Stanford University and the Authors.
+Portions copyright (c) 2013-2020 Stanford University and the Authors.
 Authors: Peter Eastman
 Contributors:
 
@@ -107,3 +107,93 @@ class MTSIntegrator(CustomIntegrator):
             else:
                 self._createSubsteps(substeps, groups[1:])
             self.addComputePerDof("v", "v+0.5*(dt/"+str(substeps)+")*f"+str(group)+"/m")
+
+
+class MTSLangevinIntegrator(CustomIntegrator):
+    """MTSLangevinIntegrator implements the BAOAB-RESPA multiple time step algorithm for
+    constant temperature dynamics.
+
+    This integrator allows different forces to be evaluated at different frequencies,
+    for example to evaluate the expensive, slowly changing forces less frequently than
+    the inexpensive, quickly changing forces.
+
+    To use it, you must first divide your forces into two or more groups (by calling
+    setForceGroup() on them) that should be evaluated at different frequencies.  When
+    you create the integrator, you provide a tuple for each group specifying the index
+    of the force group and the frequency (as a fraction of the outermost time step) at
+    which to evaluate it.  For example::
+
+        integrator = MTSLangevinIntegrator(300*kelvin, 1/picosecond, 4*femtoseconds, [(0,1), (1,2), (2,8)])
+
+    This specifies that the outermost time step is 4 fs, so each step of the integrator
+    will advance time by that much.  It also says that force group 0 should be evaluated
+    once per time step, force group 1 should be evaluated twice per time step (every 2 fs),
+    and force group 2 should be evaluated eight times per time step (every 0.5 fs).
+
+    A common use of this algorithm is to evaluate reciprocal space nonbonded interactions
+    less often than the bonded and direct space nonbonded interactions.  The following
+    example looks up the NonbondedForce, sets the reciprocal space interactions to their
+    own force group, and then creates an integrator that evaluates them once every 4 fs,
+    but all other interactions every 2 fs::
+
+        nonbonded = [f for f in system.getForces() if isinstance(f, NonbondedForce)][0]
+        nonbonded.setReciprocalSpaceForceGroup(1)
+        integrator = MTSLangevinIntegrator(300*kelvin, 1/picosecond, 4*femtoseconds, [(1,1), (0,2)])
+
+    For details, see Tuckerman et al., J. Chem. Phys. 97(3) pp. 1990-2001 (1992) and
+    Lagardere et al., J. Phys. Chem. Lett. 10(10) pp. 2593-2599 (2019).
+    """
+
+    def __init__(self, temperature, friction, dt, groups):
+        """Create an MTSLangevinIntegrator.
+
+        Parameters
+        ----------
+        temperature : temperature
+            the temperature of the heat bath
+        friction : 1/temperature
+            the friction coefficient which couples the system to the heat bath
+        dt : time
+            The largest (outermost) integration time step to use
+        groups : list
+            A list of tuples defining the force groups.  The first element of
+            each tuple is the force group index, and the second element is the
+            number of times that force group should be evaluated in one time step.
+        """
+        if len(groups) == 0:
+            raise ValueError("No force groups specified")
+        groups = sorted(groups, key=lambda x: x[1])
+        CustomIntegrator.__init__(self, dt)
+        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)))
+        from simtk.unit import MOLAR_GAS_CONSTANT_R
+        self.addGlobalVariable('kT', MOLAR_GAS_CONSTANT_R*temperature)
+        self.addPerDofVariable("x1", 0)
+        self.addUpdateContextState();
+        self._createSubsteps(1, groups)
+        self.addConstrainVelocities();
+
+    def _createSubsteps(self, parentSubsteps, groups):
+        group, substeps = groups[0]
+        stepsPerParentStep = substeps / parentSubsteps
+        if stepsPerParentStep < 1 or stepsPerParentStep != int(stepsPerParentStep):
+            raise ValueError("The number for substeps for each group must be a multiple of the number for the previous group")
+        stepsPerParentStep = int(stepsPerParentStep)
+        if group < 0 or group > 31:
+            raise ValueError("Force group must be between 0 and 31")
+        for i in range(stepsPerParentStep):
+            self.addComputePerDof("v", "v+0.5*(dt/"+str(substeps)+")*f"+str(group)+"/m")
+            if len(groups) == 1:
+                self.addComputePerDof("x", "x+(dt/"+str(2*substeps)+")*v")
+                self.addComputePerDof("v", "a*v + b*sqrt(kT/m)*gaussian")
+                self.addComputePerDof("x", "x+(dt/"+str(2*substeps)+")*v")
+                self.addComputePerDof("x1", "x")
+                self.addConstrainPositions();
+                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

@@ -121,6 +121,41 @@ class TestIntegrators(unittest.TestCase):
             context.setPositions(pdb.positions)
             integrator.step(10)
 
+    def testMTSLangevinIntegrator(self):
+        """Test the MTSLangevinIntegrator 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 = MTSLangevinIntegrator(300*kelvin, 1/picosecond, 4*femtoseconds, [(2,1), (1,2), (0,4)])
+
+        # Run some equilibration.
+        context = Context(system, integrator)
+        context.setPositions(pdb.positions)
+        context.setVelocitiesToTemperature(300*kelvin)
+        integrator.step(500)
+
+        # See if the temperature is correct.
+        totalEnergy = 0*kilojoules_per_mole
+        steps = 50
+        for i in range(steps):
+            integrator.step(10)
+            totalEnergy += context.getState(getEnergy=True).getKineticEnergy()
+        averageEnergy = totalEnergy/steps
+        dofs = 3*system.getNumParticles() - system.getNumConstraints() - 3
+        temperature = averageEnergy*2/(dofs*MOLAR_GAS_CONSTANT_R)
+        self.assertTrue(290*kelvin < temperature < 310*kelvin)
+
     def testNoseHooverIntegrator(self):
         """Test partial thermostating in the NoseHooverIntegrator (only API)"""
         pdb = PDBFile('systems/alanine-dipeptide-explicit.pdb')