Merge pull request #2561 from peastman/middle

Converted BAOABLangevinIntegrator to LangevinMiddleIntegrator

docs-source/usersguide/application.rst

@@ -120,7 +120,7 @@ steps.
         forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
         system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME,
                 nonbondedCutoff=1*nanometer, constraints=HBonds)
-        integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
+        integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
         simulation = Simulation(pdb.topology, system, integrator)
         simulation.context.setPositions(pdb.positions)
         simulation.minimizeEnergy()
@@ -210,14 +210,17 @@ convenient and less error-prone.  We could have equivalently specified
 The units system will be described in more detail later, in Section :ref:`units-and-dimensional-analysis`.
 ::
 
-    integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
+    integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
 
 This line creates the integrator to use for advancing the equations of motion.
-It specifies a :class:`BAOABLangevinIntegrator`, which performs Langevin dynamics,
+It specifies a :class:`LangevinMiddleIntegrator`, which performs Langevin dynamics,
 and assigns it to a variable called :code:`integrator`\ .  It also specifies
 the values of three parameters that are specific to Langevin dynamics: the
 simulation temperature (300 K), the friction coefficient (1 ps\ :sup:`-1`\ ), and
-the step size (0.004 ps).
+the step size (0.004 ps).  Lots of other integration methods are also available.
+For example, if you wanted to simulate the system at constant energy rather than
+constant temperature you would use a :code:`VerletIntegrator`\ .  The available
+integration methods are listed in Section :ref:`integrators`.
 ::
 
     simulation = Simulation(pdb.topology, system, integrator)
@@ -295,7 +298,7 @@ found in OpenMM’s :file:`examples` folder with the name :file:`simulateAmber.p
         inpcrd = AmberInpcrdFile('input.inpcrd')
         system = prmtop.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
                 constraints=HBonds)
-        integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
+        integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
         simulation = Simulation(prmtop.topology, system, integrator)
         simulation.context.setPositions(inpcrd.positions)
         if inpcrd.boxVectors is not None:
@@ -389,7 +392,7 @@ with the name :file:`simulateGromacs.py`.
                 includeDir='/usr/local/gromacs/share/gromacs/top')
         system = top.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
                 constraints=HBonds)
-        integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
+        integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
         simulation = Simulation(top.topology, system, integrator)
         simulation.context.setPositions(gro.positions)
         simulation.minimizeEnergy()
@@ -453,7 +456,7 @@ on the :class:`CharmmPsfFile`.
         params = CharmmParameterSet('charmm22.rtf', 'charmm22.prm')
         system = psf.createSystem(params, nonbondedMethod=NoCutoff,
                 nonbondedCutoff=1*nanometer, constraints=HBonds)
-        integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
+        integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
         simulation = Simulation(psf.topology, system, integrator)
         simulation.context.setPositions(pdb.positions)
         simulation.minimizeEnergy()
@@ -1018,36 +1021,42 @@ mass constant while slowing down the fast motions of hydrogens.  When combined
 with constraints (typically :code:`constraints=AllBonds`\ ), this often allows a
 further increase in integration step size.
 
+.. _integrators:
+
 Integrators
 ===========
 
 
 OpenMM offers a choice of several different integration methods.  You select
 which one to use by creating an integrator object of the appropriate type.
+Detailed descriptions of all these integrators can be found in Chapter
+:ref:`integrators-theory`.  In addition to these built in integrators, lots of
+others are available as part of the `OpenMMTools <https://openmmtools.readthedocs.io>`_ package.
 
-BAOAB Langevin Integrator
--------------------------
+Langevin Middle Integrator
+--------------------------
 
 In the examples of the previous sections, we used Langevin integration:
 ::
 
-    integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
+    integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
 
 The three parameter values in this line are the simulation temperature (300 K),
 the friction coefficient (1 ps\ :sup:`-1`\ ), and the step size (0.004 ps).  You
 are free to change these to whatever values you want.  Be sure to specify units
 on all values.  For example, the step size could be written either as
 :code:`0.004*picoseconds` or :code:`4*femtoseconds`\ .  They are exactly
-equivalent.
+equivalent.  Note that :code:`LangevinMiddleIntegrator` is a leapfrog
+integrator, so the velocities are offset by half a time step from the positions.
 
 Langevin Integrator
 -------------------
 
-:code:`LangevinIntegrator` is very similar to :code:`BAOABLangevinIntegrator`,
+:code:`LangevinIntegrator` is very similar to :code:`LangevinMiddleIntegrator`,
 but it uses a different discretization of the Langevin equation.
-:code:`BAOABLangevinIntegrator` tends to produce more accurate configurational
+:code:`LangevinMiddleIntegrator` tends to produce more accurate configurational
 sampling, and therefore is preferred for most applications.  Also note that
-:code:`LangevinIntegrator` (unlike :code:`BAOABLangevinIntegrator`) is a leapfrog
+:code:`LangevinIntegrator`\ , like :code:`LangevinMiddleIntegrator`\ , is a leapfrog
 integrator, so the velocities are offset by half a time step from the positions.
 
 Leapfrog Verlet Integrator
@@ -1158,7 +1167,7 @@ previous section:
     system = prmtop.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
             constraints=HBonds)
     system.addForce(MonteCarloBarostat(1*bar, 300*kelvin))
-    integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
+    integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
     ...
 
 The parameters of the Monte Carlo barostat are the pressure (1 bar) and
@@ -1718,7 +1727,7 @@ executing 1000 time steps at each temperature:
         :autonumber:`Example,simulated annealing`
 
 This code needs very little explanation.  The loop is executed 100 times.  Each
-time through, it adjusts the temperature of the :class:`BAOABLangevinIntegrator` and then
+time through, it adjusts the temperature of the :class:`LangevinMiddleIntegrator` and then
 calls :code:`step(1000)` to take 1000 time steps.
 
 Applying an External Force to Particles: a Spherical Container
@@ -1750,7 +1759,7 @@ coordinates.  Here is the code to do it:
         system.addForce(force)
         for i in range(system.getNumParticles()):
             force.addParticle(i, [])
-        integrator = BAOABLangevinIntegrator(300*kelvin, 91/picosecond, 0.004*picoseconds)
+        integrator = LangevinMiddleIntegrator(300*kelvin, 91/picosecond, 0.004*picoseconds)
         ...
 
     .. caption::

docs-source/usersguide/library.rst

@@ -243,14 +243,14 @@ simulation might look like:
         angles->addAngle(angle[i].particle1, angle[i].particle2,
             angle[i].particle3, angle[i].angle, angle[i].k);
     // ...create and initialize other force field terms in the same way
-    BAOABLangevinIntegrator integrator(temperature, friction, stepSize);
+    LangevinMiddleIntegrator integrator(temperature, friction, stepSize);
     Context context(system, integrator);
     context.setPositions(initialPositions);
     context.setVelocities(initialVelocities);
     integrator.step(10000);
 
 We create a System, add various Forces to it, and set parameters on both the
-System and the Forces.  We then create a BAOABLangevinIntegrator, initialize a
+System and the Forces.  We then create a LangevinMiddleIntegrator, initialize a
 Context in which to run a simulation, and instruct the Integrator to advance the
 simulation for 10,000 time steps.
 
@@ -1631,7 +1631,7 @@ along with the handle :code:`omm`\ , back to the calling function.
     // best available Platform. Initialize the configuration from the default
     // positions we collected above. Initial velocities will be zero but could
     // have been set here.
-    omm->integrator = new OpenMM::BAOABLangevinIntegrator(temperature,
+    omm->integrator = new OpenMM::LangevinMiddleIntegrator(temperature,
     frictionInPs,
     stepSizeInFs * OpenMM::PsPerFs);
     omm->context    = new OpenMM::Context(*omm->system, *omm->integrator);

docs-source/usersguide/references.bib

@@ -240,19 +240,6 @@
    type = {Journal Article}
 }
 
-@article{Leimkuhler2013,
-    author = {Leimkuhler, Benedict and Matthews, Charles},
-    title = {Rational Construction of Stochastic Numerical Methods for Molecular Sampling},
-    journal = {Applied Mathematics Research eXpress},
-    volume = {2013},
-    number = {1},
-    pages = {34-56},
-    year = {2012},
-    month = {06},
-    doi = {10.1093/amrx/abs010},
-    type = {Journal Article}
-}
-
 @article{Li2010
    author = {Li, D.W. and Br{\"u}schweiler, R.},
    title = {{NMR}-based protein potentials},
@@ -629,3 +616,15 @@
    year = {2015},
    type = {Journal Article}
 }
+
+@article{Zhang2019,
+    author = {Zhang, Zhijun and Liu, Xinzijian and Yan, Kangyu and Tuckerman, Mark E. and Liu, Jian},
+    title = {Unified Efficient Thermostat Scheme for the Canonical Ensemble with Holonomic or Isokinetic Constraints via Molecular Dynamics},
+    journal = {The Journal of Physical Chemistry A},
+    volume = {123},
+    number = {28},
+    pages = {6056-6079},
+    year = {2019},
+    doi = {10.1021/acs.jpca.9b02771},
+    type = {Journal Article}
+}

docs-source/usersguide/theory.rst

@@ -1278,6 +1278,8 @@ algorithm.  This can be used to implement algorithms such as lambda-dynamics,
 where a global parameter is integrated as a dynamic variable.
 
 
+.. _integrators-theory:
+
 Integrators
 ###########
 
@@ -1353,16 +1355,16 @@ random number, and :math:`\alpha=\exp(-\gamma\Delta t)`.
 The same comments about the offset between positions and velocities apply to
 this integrator as to VerletIntegrator.
 
-BAOABLangevinIntegrator
-***********************
+LangevinMiddleIntegrator
+************************
 
-This integrator is similar to LangevinIntegerator, but it instead uses the BAOAB
-discretization. :cite:`Leimkuhler2013` In each step, the positions and velocities
+This integrator is similar to LangevinIntegerator, but it instead uses the LFMiddle
+discretization. :cite:`Zhang2019` In each step, the positions and velocities
 are updated as follows:
 
 
 .. math::
-   \mathbf{v}_{i}(t+\Delta t/2) = \mathbf{v}_{i}(t) + \mathbf{f}_{i}(t)\Delta t/2{m}_{i}
+   \mathbf{v}_{i}(t+\Delta t/2) = \mathbf{v}_{i}(t-\Delta t/2) + \mathbf{f}_{i}(t)\Delta t/{m}_{i}
 
 
 .. math::
@@ -1377,26 +1379,18 @@ are updated as follows:
    \mathbf{r}_{i}(t+\Delta t) = \mathbf{r}_{i}(t+\Delta t/2) + \mathbf{v'}_{i}(t+\Delta t/2)\Delta t/2
 
 
-.. math::
-   \mathbf{v}_{i}(t+\Delta t) = \mathbf{v'}_{i}(t+\Delta t/2) + \mathbf{f}_{i}(t+\Delta t)\Delta t/2{m}_{i}
-
-
-This tends to produce more accurate
-sampling of configurational properties (such as free energies), but less
-accurate sampling of kinetic properties (such as mean kinetic energy).  Because
+This tends to produce more accurate sampling of configurational properties (such
+as free energies), but less accurate sampling of kinetic properties.  Because
 configurational properties are much more important than kinetic ones in most
 simulations, this integrator is generally preferred over LangevinIntegrator.  It
 often allows one to use a larger time step while still maintaining similar or
 better accuracy.
 
 One disadvantage of this integrator is that it requires applying constraints
-three times per time step, compared to only once for LangevinIntegrator.  This
+twice per time step, compared to only once for LangevinIntegrator.  This
 can make it slightly slower for systems that involve constraints.  However, this
 usually is more than compensated by allowing you to use a larger time step.
 
-Unlike LangevinIntegrator, this does not use a leap-frog algorithm.  The
-positions and velocities all correspond to the same point in time.
-
 BrownianIntegrator
 ******************
 

examples/HelloSodiumChloride.cpp

@@ -249,7 +249,7 @@ myInitializeOpenMM( const MyAtomInfo    atoms[],
     // best available Platform. Initialize the configuration from the default
     // positions we collected above. Initial velocities will be zero but could
     // have been set here.
-    omm->integrator = new OpenMM::BAOABLangevinIntegrator(temperature, frictionInPs, 
+    omm->integrator = new OpenMM::LangevinMiddleIntegrator(temperature, frictionInPs, 
                                                      stepSizeInFs * OpenMM::PsPerFs);
     omm->context    = new OpenMM::Context(*omm->system, *omm->integrator);
     omm->context->setPositions(initialPosInNm);

examples/HelloSodiumChlorideInC.c

@@ -252,7 +252,7 @@ myInitializeOpenMM( const MyAtomInfo    atoms[],
      * best available Platform. Initialize the configuration from the default
      * positions we collected above. Initial velocities will be zero but could
      * have been set here. */
-    omm->integrator = (OpenMM_Integrator*)OpenMM_BAOABLangevinIntegrator_create(
+    omm->integrator = (OpenMM_Integrator*)OpenMM_LangevinMiddleIntegrator_create(
                                             temperature, frictionInPerPs, 
                                             stepSizeInFs * OpenMM_PsPerFs);
     omm->context    = OpenMM_Context_create(omm->system, omm->integrator);

examples/HelloSodiumChlorideInFortran.f90

@@ -171,7 +171,7 @@ SUBROUTINE myInitializeOpenMM(ommHandle, platformName)
     ! These are the objects we'll create here thare are stored in the
     ! Context for later access. Don't forget to delete them at the end.
     type (OpenMM_System)             system
-    type (OpenMM_BAOABLangevinIntegrator) langevin
+    type (OpenMM_LangevinMiddleIntegrator) langevin
     type (OpenMM_Context)            context
 
     ! These are temporary OpenMM objects used and discarded here.
@@ -236,11 +236,11 @@ SUBROUTINE myInitializeOpenMM(ommHandle, platformName)
     ! best available Platform. Initialize the configuration from the default
     ! positions we collected above. Initial velocities will be zero but could
     ! have been set here.
-    call OpenMM_BAOABLangevinIntegrator_create(langevin,                     &
+    call OpenMM_LangevinMiddleIntegrator_create(langevin,                     &
                                           Temperature, FrictionInPerPs, &
                                           StepSizeInFs * OpenMM_PsPerFs)
 
-    ! Convert BAOABLangevinIntegrator to generic Integrator type for this call.
+    ! Convert LangevinMiddleIntegrator to generic Integrator type for this call.
     call OpenMM_Context_create(context, system,                         &
                                transfer(langevin, OpenMM_Integrator(0)))
     call OpenMM_Context_setPositions(context, initialPosInNm)

examples/simulateAmber.py

@@ -6,7 +6,7 @@ from sys import stdout
 prmtop = AmberPrmtopFile('input.prmtop')
 inpcrd = AmberInpcrdFile('input.inpcrd')
 system = prmtop.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer, constraints=HBonds)
-integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
+integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
 simulation = Simulation(prmtop.topology, system, integrator)
 simulation.context.setPositions(inpcrd.positions)
 if inpcrd.boxVectors is not None:

examples/simulateCharmm.py

@@ -21,7 +21,7 @@ params = CharmmParameterSet('charmm22.rtf', 'charmm22.par')
 
 # Instantiate the system
 system = psf.createSystem(params, nonbondedMethod=NoCutoff, constraints=HBonds)
-integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
+integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
 simulation = Simulation(psf.topology, system, integrator)
 simulation.context.setPositions(pdb.getPositions())
 simulation.minimizeEnergy()

examples/simulateGromacs.py

@@ -6,7 +6,7 @@ from sys import stdout
 gro = GromacsGroFile('input.gro')
 top = GromacsTopFile('input.top', periodicBoxVectors=gro.getPeriodicBoxVectors())
 system = top.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer, constraints=HBonds)
-integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
+integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
 simulation = Simulation(top.topology, system, integrator)
 simulation.context.setPositions(gro.positions)
 simulation.minimizeEnergy()

examples/simulatePdb.py

@@ -6,7 +6,7 @@ from sys import stdout
 pdb = PDBFile('input.pdb')
 forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
 system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCutoff=1*nanometer, constraints=HBonds)
-integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
+integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
 simulation = Simulation(pdb.topology, system, integrator)
 simulation.context.setPositions(pdb.positions)
 simulation.minimizeEnergy()

olla/include/openmm/kernels.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -33,7 +33,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "openmm/AndersenThermostat.h"
-#include "openmm/BAOABLangevinIntegrator.h"
+#include "openmm/LangevinMiddleIntegrator.h"
 #include "openmm/BrownianIntegrator.h"
 #include "openmm/CMAPTorsionForce.h"
 #include "openmm/CMMotionRemover.h"
@@ -1129,40 +1129,36 @@ public:
 };
 
 /**
- * This kernel is invoked by BAOABLangevinIntegrator to take one time step.
+ * This kernel is invoked by LangevinMiddleIntegrator to take one time step.
  */
-class IntegrateBAOABStepKernel : public KernelImpl {
+class IntegrateLangevinMiddleStepKernel : public KernelImpl {
 public:
     static std::string Name() {
-        return "IntegrateBAOABStep";
+        return "IntegrateLangevinMiddleStep";
     }
-    IntegrateBAOABStepKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    IntegrateLangevinMiddleStepKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
     }
     /**
      * Initialize the kernel.
      * 
      * @param system     the System this kernel will be applied to
-     * @param integrator the BAOABLangevinIntegrator this kernel will be used for
+     * @param integrator the LangevinMiddleIntegrator this kernel will be used for
      */
-    virtual void initialize(const System& system, const BAOABLangevinIntegrator& integrator) = 0;
+    virtual void initialize(const System& system, const LangevinMiddleIntegrator& integrator) = 0;
     /**
      * Execute the kernel.
      * 
      * @param context    the context in which to execute this kernel
-     * @param integrator the BAOABLangevinIntegrator this kernel is being used for
-     * @param forcesAreValid if the context has been modified since the last time step, this will be
-     *                       false to show that cached forces are invalid and must be recalculated.
-     *                       On exit, this should specify whether the cached forces are valid at the
-     *                       end of the step.
+     * @param integrator the LangevinMiddleIntegrator this kernel is being used for
      */
-    virtual void execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator, bool& forcesAreValid) = 0;
+    virtual void execute(ContextImpl& context, const LangevinMiddleIntegrator& integrator) = 0;
     /**
      * Compute the kinetic energy.
      * 
      * @param context    the context in which to execute this kernel
-     * @param integrator the BAOABLangevinIntegrator this kernel is being used for
+     * @param integrator the LangevinMiddleIntegrator this kernel is being used for
      */
-    virtual double computeKineticEnergy(ContextImpl& context, const BAOABLangevinIntegrator& integrator) = 0;
+    virtual double computeKineticEnergy(ContextImpl& context, const LangevinMiddleIntegrator& integrator) = 0;
 };
 
 /**

openmmapi/include/OpenMM.h

@@ -33,7 +33,6 @@
  * -------------------------------------------------------------------------- */
 
 #include "openmm/AndersenThermostat.h"
-#include "openmm/BAOABLangevinIntegrator.h"
 #include "openmm/BrownianIntegrator.h"
 #include "openmm/CMAPTorsionForce.h"
 #include "openmm/CMMotionRemover.h"
@@ -57,6 +56,7 @@
 #include "openmm/HarmonicBondForce.h"
 #include "openmm/Integrator.h"
 #include "openmm/LangevinIntegrator.h"
+#include "openmm/LangevinMiddleIntegrator.h"
 #include "openmm/LocalEnergyMinimizer.h"
 #include "openmm/MonteCarloAnisotropicBarostat.h"
 #include "openmm/MonteCarloBarostat.h"

openmmapi/include/openmm/CustomIntegrator.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2011-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -236,7 +236,7 @@ namespace OpenMM {
  * </pre></tt>
  * 
  * The second one implements the algorithm used by the standard
- * BAOABLangevinIntegrator class.  kB is Boltzmann's constant.
+ * LangevinMiddleIntegrator class.  kB is Boltzmann's constant.
  * 
  * <tt><pre>
  * CustomIntegrator integrator(dt);
@@ -245,18 +245,14 @@ namespace OpenMM {
  * integrator.addGlobalVariable("kT", kB*temperature);
  * integrator.addPerDofVariable("x1", 0);
  * integrator.addUpdateContextState();
- * integrator.addComputePerDof("v", "v + 0.5*dt*f/m");
+ * integrator.addComputePerDof("v", "v + dt*f/m");
+ * integrator.addConstrainVelocities();
  * integrator.addComputePerDof("x", "x + 0.5*dt*v");
- * integrator.addComputePerDof("x1", "x");
- * integrator.addConstrainPositions();
- * integrator.addComputePerDof("v", "v + 2*(x-x1)/dt");
  * integrator.addComputePerDof("v", "a*v + b*sqrt(kT/m)*gaussian");
  * integrator.addComputePerDof("x", "x + 0.5*dt*v");
  * integrator.addComputePerDof("x1", "x");
  * integrator.addConstrainPositions();
- * integrator.addComputePerDof("v", "v + 2*(x-x1)/dt");
- * integrator.addComputePerDof("v", "v + 0.5*dt*f/m");
- * integrator.addConstrainVelocities();
+ * integrator.addComputePerDof("v", "v + (x-x1)/dt");
  * </pre></tt>
  * 
  * Another feature of CustomIntegrator is that it can use derivatives of the

openmmapi/include/openmm/BAOABLangevinIntegrator.h → openmmapi/include/openmm/LangevinMiddleIntegrator.h

-#ifndef OPENMM_BAOABLANGEVININTEGRATOR_H_
-#define OPENMM_BAOABLANGEVININTEGRATOR_H_
+#ifndef OPENMM_LANGEVINMIDDLEINTEGRATOR_H_
+#define OPENMM_LANGEVINMIDDLEINTEGRATOR_H_
 
 /* -------------------------------------------------------------------------- *
  *                                   OpenMM                                   *
@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -40,21 +40,21 @@ namespace OpenMM {
 
 /**
  * This is an Integrator which simulates a System using Langevin dynamics, with
- * the BAOAB discretization of Leimkuhler and Matthews (http://dx.doi.org/10.1093/amrx/abs010).
+ * the LFMiddle discretization (http://dx.doi.org/10.1021/acs.jpca.9b02771).
  * This method tend to produce more accurate configurational sampling than other
  * discretizations, such as the one used in LangevinIntegrator.
  */
 
-class OPENMM_EXPORT BAOABLangevinIntegrator : public Integrator {
+class OPENMM_EXPORT LangevinMiddleIntegrator : public Integrator {
 public:
     /**
-     * Create a BAOABLangevinIntegrator.
+     * Create a LangevinMiddleIntegrator.
      * 
      * @param temperature    the temperature of the heat bath (in Kelvin)
      * @param frictionCoeff  the friction coefficient which couples the system to the heat bath (in inverse picoseconds)
      * @param stepSize       the step size with which to integrate the system (in picoseconds)
      */
-    BAOABLangevinIntegrator(double temperature, double frictionCoeff, double stepSize);
+    LangevinMiddleIntegrator(double temperature, double frictionCoeff, double stepSize);
     /**
      * Get the temperature of the heat bath (in Kelvin).
      *
@@ -128,10 +128,6 @@ protected:
      * cleanup.  It will also get called again if the application calls reinitialize() on the Context.
      */
     void cleanup();
-    /**
-     * When the user modifies the state, we need to mark that the forces need to be recalculated.
-     */
-    void stateChanged(State::DataType changed);
     /**
      * Get the names of all Kernels used by this Integrator.
      */
@@ -147,10 +143,9 @@ protected:
 private:
     double temperature, friction;
     int randomNumberSeed;
-    bool forcesAreValid;
     Kernel kernel;
 };
 
 } // namespace OpenMM
 
-#endif /*OPENMM_BAOABLANGEVININTEGRATOR_H_*/
+#endif /*OPENMM_LANGEVINMIDDLEINTEGRATOR_H_*/

openmmapi/src/BAOABLangevinIntegrator.cpp → openmmapi/src/LangevinMiddleIntegrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -29,7 +29,7 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
-#include "openmm/BAOABLangevinIntegrator.h"
+#include "openmm/LangevinMiddleIntegrator.h"
 #include "openmm/Context.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/internal/ContextImpl.h"
@@ -40,53 +40,47 @@ using namespace OpenMM;
 using std::string;
 using std::vector;
 
-BAOABLangevinIntegrator::BAOABLangevinIntegrator(double temperature, double frictionCoeff, double stepSize) {
+LangevinMiddleIntegrator::LangevinMiddleIntegrator(double temperature, double frictionCoeff, double stepSize) {
     setTemperature(temperature);
     setFriction(frictionCoeff);
     setStepSize(stepSize);
     setConstraintTolerance(1e-5);
     setRandomNumberSeed(0);
-    forcesAreValid = false;
 }
 
-void BAOABLangevinIntegrator::initialize(ContextImpl& contextRef) {
+void LangevinMiddleIntegrator::initialize(ContextImpl& contextRef) {
     if (owner != NULL && &contextRef.getOwner() != owner)
         throw OpenMMException("This Integrator is already bound to a context");
     context = &contextRef;
     owner = &contextRef.getOwner();
-    kernel = context->getPlatform().createKernel(IntegrateBAOABStepKernel::Name(), contextRef);
-    kernel.getAs<IntegrateBAOABStepKernel>().initialize(contextRef.getSystem(), *this);
+    kernel = context->getPlatform().createKernel(IntegrateLangevinMiddleStepKernel::Name(), contextRef);
+    kernel.getAs<IntegrateLangevinMiddleStepKernel>().initialize(contextRef.getSystem(), *this);
 }
 
-void BAOABLangevinIntegrator::cleanup() {
+void LangevinMiddleIntegrator::cleanup() {
     kernel = Kernel();
 }
 
-void BAOABLangevinIntegrator::stateChanged(State::DataType changed) {
-    forcesAreValid = false;
-}
-
-vector<string> BAOABLangevinIntegrator::getKernelNames() {
+vector<string> LangevinMiddleIntegrator::getKernelNames() {
     std::vector<std::string> names;
-    names.push_back(IntegrateBAOABStepKernel::Name());
+    names.push_back(IntegrateLangevinMiddleStepKernel::Name());
     return names;
 }
 
-double BAOABLangevinIntegrator::computeKineticEnergy() {
-    forcesAreValid = false;
-    return kernel.getAs<IntegrateBAOABStepKernel>().computeKineticEnergy(*context, *this);
+double LangevinMiddleIntegrator::computeKineticEnergy() {
+    return kernel.getAs<IntegrateLangevinMiddleStepKernel>().computeKineticEnergy(*context, *this);
 }
 
-bool BAOABLangevinIntegrator::kineticEnergyRequiresForce() const {
+bool LangevinMiddleIntegrator::kineticEnergyRequiresForce() const {
     return false;
 }
 
-void BAOABLangevinIntegrator::step(int steps) {
+void LangevinMiddleIntegrator::step(int steps) {
     if (context == NULL)
         throw OpenMMException("This Integrator is not bound to a context!");  
     for (int i = 0; i < steps; ++i) {
-        if (context->updateContextState())
-            forcesAreValid = false;
-        kernel.getAs<IntegrateBAOABStepKernel>().execute(*context, *this, forcesAreValid);
+        context->updateContextState();
+        context->calcForcesAndEnergy(true, false);
+        kernel.getAs<IntegrateLangevinMiddleStepKernel>().execute(*context, *this);
     }
 }

platforms/common/include/openmm/common/CommonKernels.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -905,44 +905,40 @@ private:
 };
 
 /**
- * This kernel is invoked by BAOABLangevinIntegrator to take one time step.
+ * This kernel is invoked by LangevinMiddleIntegrator to take one time step.
  */
-class CommonIntegrateBAOABStepKernel : public IntegrateBAOABStepKernel {
+class CommonIntegrateLangevinMiddleStepKernel : public IntegrateLangevinMiddleStepKernel {
 public:
-    CommonIntegrateBAOABStepKernel(std::string name, const Platform& platform, ComputeContext& cc) : IntegrateBAOABStepKernel(name, platform), cc(cc),
+    CommonIntegrateLangevinMiddleStepKernel(std::string name, const Platform& platform, ComputeContext& cc) : IntegrateLangevinMiddleStepKernel(name, platform), cc(cc),
             hasInitializedKernels(false) {
     }
     /**
      * Initialize the kernel, setting up the particle masses.
      * 
      * @param system     the System this kernel will be applied to
-     * @param integrator the BAOABLangevinIntegrator this kernel will be used for
+     * @param integrator the LangevinMiddleIntegrator this kernel will be used for
      */
-    void initialize(const System& system, const BAOABLangevinIntegrator& integrator);
+    void initialize(const System& system, const LangevinMiddleIntegrator& integrator);
     /**
      * Execute the kernel.
      * 
      * @param context    the context in which to execute this kernel
-     * @param integrator the BAOABLangevinIntegrator this kernel is being used for
-     * @param forcesAreValid if the context has been modified since the last time step, this will be
-     *                       false to show that cached forces are invalid and must be recalculated.
-     *                       On exit, this should specify whether the cached forces are valid at the
-     *                       end of the step.
+     * @param integrator the LangevinMiddleIntegrator this kernel is being used for
      */
-    void execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator, bool& forcesAreValid);
+    void execute(ContextImpl& context, const LangevinMiddleIntegrator& integrator);
     /**
      * Compute the kinetic energy.
      * 
      * @param context    the context in which to execute this kernel
-     * @param integrator the BAOABLangevinIntegrator this kernel is being used for
+     * @param integrator the LangevinMiddleIntegrator this kernel is being used for
      */
-    double computeKineticEnergy(ContextImpl& context, const BAOABLangevinIntegrator& integrator);
+    double computeKineticEnergy(ContextImpl& context, const LangevinMiddleIntegrator& integrator);
 private:
     ComputeContext& cc;
     double prevTemp, prevFriction, prevStepSize;
     bool hasInitializedKernels;
     ComputeArray params, oldDelta;
-    ComputeKernel kernel1, kernel2, kernel3, kernel4;
+    ComputeKernel kernel1, kernel2, kernel3;
 };
 
 /**

platforms/common/src/CommonKernels.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -5547,27 +5547,26 @@ double CommonIntegrateLangevinStepKernel::computeKineticEnergy(ContextImpl& cont
     return cc.getIntegrationUtilities().computeKineticEnergy(0.5*integrator.getStepSize());
 }
 
-void CommonIntegrateBAOABStepKernel::initialize(const System& system, const BAOABLangevinIntegrator& integrator) {
+void CommonIntegrateLangevinMiddleStepKernel::initialize(const System& system, const LangevinMiddleIntegrator& integrator) {
     cc.initializeContexts();
     cc.setAsCurrent();
     cc.getIntegrationUtilities().initRandomNumberGenerator(integrator.getRandomNumberSeed());
-    ComputeProgram program = cc.compileProgram(CommonKernelSources::baoab);
-    kernel1 = program->createKernel("integrateBAOABPart1");
-    kernel2 = program->createKernel("integrateBAOABPart2");
-    kernel3 = program->createKernel("integrateBAOABPart3");
-    kernel4 = program->createKernel("integrateBAOABPart4");
+    ComputeProgram program = cc.compileProgram(CommonKernelSources::langevinMiddle);
+    kernel1 = program->createKernel("integrateLangevinMiddlePart1");
+    kernel2 = program->createKernel("integrateLangevinMiddlePart2");
+    kernel3 = program->createKernel("integrateLangevinMiddlePart3");
     if (cc.getUseDoublePrecision() || cc.getUseMixedPrecision()) {
-        params.initialize<double>(cc, 2, "baoabParams");
+        params.initialize<double>(cc, 2, "langevinMiddleParams");
         oldDelta.initialize<mm_double4>(cc, cc.getPaddedNumAtoms(), "oldDelta");
     }
     else {
-        params.initialize<float>(cc, 2, "baoabParams");
+        params.initialize<float>(cc, 2, "langevinMiddleParams");
         oldDelta.initialize<mm_float4>(cc, cc.getPaddedNumAtoms(), "oldDelta");
     }
     prevStepSize = -1.0;
 }
 
-void CommonIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator, bool& forcesAreValid) {
+void CommonIntegrateLangevinMiddleStepKernel::execute(ContextImpl& context, const LangevinMiddleIntegrator& integrator) {
     IntegrationUtilities& integration = cc.getIntegrationUtilities();
     int numAtoms = cc.getNumAtoms();
     int paddedNumAtoms = cc.getPaddedNumAtoms();
@@ -5577,11 +5576,8 @@ void CommonIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLa
         kernel1->addArg(paddedNumAtoms);
         kernel1->addArg(cc.getVelm());
         kernel1->addArg(cc.getLongForceBuffer());
-        kernel1->addArg(integration.getPosDelta());
-        kernel1->addArg(oldDelta);
         kernel1->addArg(integration.getStepSize());
         kernel2->addArg(numAtoms);
-        kernel2->addArg(cc.getPosq());
         kernel2->addArg(cc.getVelm());
         kernel2->addArg(integration.getPosDelta());
         kernel2->addArg(oldDelta);
@@ -5589,8 +5585,6 @@ void CommonIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLa
         kernel2->addArg(integration.getStepSize());
         kernel2->addArg(integration.getRandom());
         kernel2->addArg(); // Random index will be set just before it is executed.
-        if (cc.getUseMixedPrecision())
-            kernel2->addArg(cc.getPosqCorrection());
         kernel3->addArg(numAtoms);
         kernel3->addArg(cc.getPosq());
         kernel3->addArg(cc.getVelm());
@@ -5599,15 +5593,6 @@ void CommonIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLa
         kernel3->addArg(integration.getStepSize());
         if (cc.getUseMixedPrecision())
             kernel3->addArg(cc.getPosqCorrection());
-        kernel4->addArg(numAtoms);
-        kernel4->addArg(paddedNumAtoms);
-        kernel4->addArg(cc.getVelm());
-        kernel4->addArg(cc.getLongForceBuffer());
-        kernel4->addArg(integration.getStepSize());
-    }
-    if (!forcesAreValid) {
-        context.calcForcesAndEnergy(true, false);
-        forcesAreValid = true;
     }
     double temperature = integrator.getTemperature();
     double friction = integrator.getFriction();
@@ -5630,15 +5615,12 @@ void CommonIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLa
 
     // Perform the integration.
 
-    kernel2->setArg(8, integration.prepareRandomNumbers(cc.getPaddedNumAtoms()));
+    kernel2->setArg(7, integration.prepareRandomNumbers(cc.getPaddedNumAtoms()));
     kernel1->execute(numAtoms);
-    integration.applyConstraints(integrator.getConstraintTolerance());
+    integration.applyVelocityConstraints(integrator.getConstraintTolerance());
     kernel2->execute(numAtoms);
     integration.applyConstraints(integrator.getConstraintTolerance());
     kernel3->execute(numAtoms);
-    context.calcForcesAndEnergy(true, false);
-    kernel4->execute(numAtoms);
-    integration.applyVelocityConstraints(integrator.getConstraintTolerance());
     integration.computeVirtualSites();
 
     // Update the time and step count.
@@ -5646,8 +5628,6 @@ void CommonIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLa
     cc.setTime(cc.getTime()+stepSize);
     cc.setStepCount(cc.getStepCount()+1);
     cc.reorderAtoms();
-    if (cc.getAtomsWereReordered())
-        forcesAreValid = false;
     
     // Reduce UI lag.
     
@@ -5656,7 +5636,7 @@ void CommonIntegrateBAOABStepKernel::execute(ContextImpl& context, const BAOABLa
 #endif
 }
 
-double CommonIntegrateBAOABStepKernel::computeKineticEnergy(ContextImpl& context, const BAOABLangevinIntegrator& integrator) {
+double CommonIntegrateLangevinMiddleStepKernel::computeKineticEnergy(ContextImpl& context, const LangevinMiddleIntegrator& integrator) {
     return cc.getIntegrationUtilities().computeKineticEnergy(0.0);
 }
 

platforms/common/src/kernels/baoab.cc → platforms/common/src/kernels/langevinMiddle.cc

 enum {VelScale, NoiseScale};
 
 /**
- * Perform the first part of BAOAB integration: velocity half step, then position half step.
+ * Perform the first part of integration: velocity step.
  */
 
-KERNEL void integrateBAOABPart1(int numAtoms, int paddedNumAtoms, GLOBAL mixed4* RESTRICT velm, GLOBAL const mm_long* RESTRICT force, GLOBAL mixed4* RESTRICT posDelta,
-        GLOBAL mixed4* RESTRICT oldDelta, GLOBAL const mixed2* RESTRICT dt) {
-    mixed halfdt = 0.5*dt[0].y;
-    mixed fscale = halfdt/(mixed) 0x100000000;
+KERNEL void integrateLangevinMiddlePart1(int numAtoms, int paddedNumAtoms, GLOBAL mixed4* RESTRICT velm, GLOBAL const mm_long* RESTRICT force,
+        GLOBAL const mixed2* RESTRICT dt) {
+    mixed fscale = dt[0].y/(mixed) 0x100000000;
     for (int index = GLOBAL_ID; index < numAtoms; index += GLOBAL_SIZE) {
         mixed4 velocity = velm[index];
         if (velocity.w != 0.0) {
@@ -15,59 +14,33 @@ KERNEL void integrateBAOABPart1(int numAtoms, int paddedNumAtoms, GLOBAL mixed4*
             velocity.y += fscale*velocity.w*force[index+paddedNumAtoms];
             velocity.z += fscale*velocity.w*force[index+paddedNumAtoms*2];
             velm[index] = velocity;
-            mixed4 delta = make_mixed4(halfdt*velocity.x, halfdt*velocity.y, halfdt*velocity.z, 0);
-            posDelta[index] = delta;
-            oldDelta[index] = delta;
         }
     }
 }
 
 /**
- * Perform the second part of BAOAB integration: apply constraint forces to velocities, then interact with heat bath,
- * then position half step.
+ * Perform the second part of integration: position half step, then interact with heat bath,
+ * then another position half step.
  */
 
-KERNEL void integrateBAOABPart2(int numAtoms, GLOBAL real4* RESTRICT posq, GLOBAL mixed4* RESTRICT velm, GLOBAL mixed4* RESTRICT posDelta,
+KERNEL void integrateLangevinMiddlePart2(int numAtoms, GLOBAL mixed4* RESTRICT velm, GLOBAL mixed4* RESTRICT posDelta,
         GLOBAL mixed4* RESTRICT oldDelta, GLOBAL const mixed* RESTRICT paramBuffer, GLOBAL const mixed2* RESTRICT dt, GLOBAL const float4* RESTRICT random, unsigned int randomIndex
-#ifdef USE_MIXED_PRECISION
-        , GLOBAL real4* RESTRICT posqCorrection
-#endif
         ) {
     mixed vscale = paramBuffer[VelScale];
     mixed noisescale = paramBuffer[NoiseScale];
     mixed halfdt = 0.5*dt[0].y;
-    mixed invHalfdt = 1/halfdt;
     int index = GLOBAL_ID;
     randomIndex += index;
     while (index < numAtoms) {
         mixed4 velocity = velm[index];
         if (velocity.w != 0.0) {
-            mixed4 delta = posDelta[index];
+            mixed4 delta = make_mixed4(halfdt*velocity.x, halfdt*velocity.y, halfdt*velocity.z, 0);
             mixed sqrtInvMass = SQRT(velocity.w);
-            velocity.x += (delta.x-oldDelta[index].x)*invHalfdt;
-            velocity.y += (delta.y-oldDelta[index].y)*invHalfdt;
-            velocity.z += (delta.z-oldDelta[index].z)*invHalfdt;
             velocity.x = vscale*velocity.x + noisescale*sqrtInvMass*random[randomIndex].x;
             velocity.y = vscale*velocity.y + noisescale*sqrtInvMass*random[randomIndex].y;
             velocity.z = vscale*velocity.z + noisescale*sqrtInvMass*random[randomIndex].z;
             velm[index] = velocity;
-#ifdef USE_MIXED_PRECISION
-            real4 pos1 = posq[index];
-            real4 pos2 = posqCorrection[index];
-            mixed4 pos = make_mixed4(pos1.x+(mixed)pos2.x, pos1.y+(mixed)pos2.y, pos1.z+(mixed)pos2.z, pos1.w);
-#else
-            real4 pos = posq[index];
-#endif
-            pos.x += delta.x;
-            pos.y += delta.y;
-            pos.z += delta.z;
-#ifdef USE_MIXED_PRECISION
-            posq[index] = make_real4((real) pos.x, (real) pos.y, (real) pos.z, (real) pos.w);
-            posqCorrection[index] = make_real4(pos.x-(real) pos.x, pos.y-(real) pos.y, pos.z-(real) pos.z, 0);
-#else
-            posq[index] = pos;
-#endif
-            delta = make_mixed4(halfdt*velocity.x, halfdt*velocity.y, halfdt*velocity.z, 0);
+            delta += make_mixed4(halfdt*velocity.x, halfdt*velocity.y, halfdt*velocity.z, 0);
             posDelta[index] = delta;
             oldDelta[index] = delta;
         }
@@ -77,25 +50,24 @@ KERNEL void integrateBAOABPart2(int numAtoms, GLOBAL real4* RESTRICT posq, GLOBA
 }
 
 /**
- * Perform the third part of BAOAB integration: apply constraint forces to velocities, then record
- * the constrained positions in preparation for computing forces.
+ * Perform the third part of integration: apply constraint forces to velocities, then record
+ * the constrained positions.
  */
 
-KERNEL void integrateBAOABPart3(int numAtoms, GLOBAL real4* RESTRICT posq, GLOBAL mixed4* RESTRICT velm,
+KERNEL void integrateLangevinMiddlePart3(int numAtoms, GLOBAL real4* RESTRICT posq, GLOBAL mixed4* RESTRICT velm,
          GLOBAL mixed4* RESTRICT posDelta, GLOBAL mixed4* RESTRICT oldDelta, GLOBAL const mixed2* RESTRICT dt
 #ifdef USE_MIXED_PRECISION
         , GLOBAL real4* RESTRICT posqCorrection
 #endif
         ) {
-    mixed halfdt = 0.5*dt[0].y;
-    mixed invHalfdt = 1/halfdt;
+    mixed invDt = 1/dt[0].y;
     for (int index = GLOBAL_ID; index < numAtoms; index += GLOBAL_SIZE) {
         mixed4 velocity = velm[index];
         if (velocity.w != 0.0) {
             mixed4 delta = posDelta[index];
-            velocity.x += (delta.x-oldDelta[index].x)*invHalfdt;
-            velocity.y += (delta.y-oldDelta[index].y)*invHalfdt;
-            velocity.z += (delta.z-oldDelta[index].z)*invHalfdt;
+            velocity.x += (delta.x-oldDelta[index].x)*invDt;
+            velocity.y += (delta.y-oldDelta[index].y)*invDt;
+            velocity.z += (delta.z-oldDelta[index].z)*invDt;
             velm[index] = velocity;
 #ifdef USE_MIXED_PRECISION
             real4 pos1 = posq[index];
@@ -116,22 +88,3 @@ KERNEL void integrateBAOABPart3(int numAtoms, GLOBAL real4* RESTRICT posq, GLOBA
         }
     }
 }
-
-/**
- * Perform the fourth part of BAOAB integration: velocity half step.
- */
-
-KERNEL void integrateBAOABPart4(int numAtoms, int paddedNumAtoms, GLOBAL mixed4* RESTRICT velm,
-        GLOBAL const mm_long* RESTRICT force, GLOBAL const mixed2* RESTRICT dt) {
-    mixed halfdt = 0.5*dt[0].y;
-    mixed fscale = halfdt/(mixed) 0x100000000;
-    for (int index = GLOBAL_ID; index < numAtoms; index += GLOBAL_SIZE) {
-        mixed4 velocity = velm[index];
-        if (velocity.w != 0.0) {
-            velocity.x += fscale*velocity.w*force[index];
-            velocity.y += fscale*velocity.w*force[index+paddedNumAtoms];
-            velocity.z += fscale*velocity.w*force[index+paddedNumAtoms*2];
-            velm[index] = velocity;
-        }
-    }
-}

platforms/cpu/include/CpuKernels.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2013-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2013-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -32,7 +32,6 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
-#include "CpuBAOABDynamics.h"
 #include "CpuBondForce.h"
 #include "CpuCustomGBForce.h"
 #include "CpuCustomManyParticleForce.h"
@@ -40,6 +39,7 @@
 #include "CpuGayBerneForce.h"
 #include "CpuGBSAOBCForce.h"
 #include "CpuLangevinDynamics.h"
+#include "CpuLangevinMiddleDynamics.h"
 #include "CpuNeighborList.h"
 #include "CpuNonbondedForce.h"
 #include "CpuPlatform.h"
@@ -538,42 +538,38 @@ private:
 };
 
 /**
- * This kernel is invoked by BAOABLangevinIntegrator to take one time step.
+ * This kernel is invoked by LangevinMiddleIntegrator to take one time step.
  */
-class CpuIntegrateBAOABStepKernel : public IntegrateBAOABStepKernel {
+class CpuIntegrateLangevinMiddleStepKernel : public IntegrateLangevinMiddleStepKernel {
 public:
-    CpuIntegrateBAOABStepKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data) : IntegrateBAOABStepKernel(name, platform),
+    CpuIntegrateLangevinMiddleStepKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data) : IntegrateLangevinMiddleStepKernel(name, platform),
             data(data), dynamics(0) {
     }
-    ~CpuIntegrateBAOABStepKernel();
+    ~CpuIntegrateLangevinMiddleStepKernel();
     /**
      * Initialize the kernel, setting up the particle masses.
      * 
      * @param system     the System this kernel will be applied to
-     * @param integrator the BAOABLangevinIntegrator this kernel will be used for
+     * @param integrator the LangevinMiddleIntegrator this kernel will be used for
      */
-    void initialize(const System& system, const BAOABLangevinIntegrator& integrator);
+    void initialize(const System& system, const LangevinMiddleIntegrator& integrator);
     /**
      * Execute the kernel.
      * 
      * @param context    the context in which to execute this kernel
-     * @param integrator the BAOABLangevinIntegrator this kernel is being used for
-     * @param forcesAreValid if the context has been modified since the last time step, this will be
-     *                       false to show that cached forces are invalid and must be recalculated.
-     *                       On exit, this should specify whether the cached forces are valid at the
-     *                       end of the step.
+     * @param integrator the LangevinMiddleIntegrator this kernel is being used for
      */
-    void execute(ContextImpl& context, const BAOABLangevinIntegrator& integrator, bool& forcesAreValid);
+    void execute(ContextImpl& context, const LangevinMiddleIntegrator& integrator);
     /**
      * Compute the kinetic energy.
      * 
      * @param context    the context in which to execute this kernel
-     * @param integrator the BAOABLangevinIntegrator this kernel is being used for
+     * @param integrator the LangevinMiddleIntegrator this kernel is being used for
      */
-    double computeKineticEnergy(ContextImpl& context, const BAOABLangevinIntegrator& integrator);
+    double computeKineticEnergy(ContextImpl& context, const LangevinMiddleIntegrator& integrator);
 private:
     CpuPlatform::PlatformData& data;
-    CpuBAOABDynamics* dynamics;
+    CpuLangevinMiddleDynamics* dynamics;
     std::vector<double> masses;
     double prevTemp, prevFriction, prevStepSize;
 };