Updated documentation and examples to include BAOABLangevinIntegrator

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 = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
+        integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
         simulation = Simulation(pdb.topology, system, integrator)
         simulation.context.setPositions(pdb.positions)
         simulation.minimizeEnergy()
@@ -210,10 +210,10 @@ 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 = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
+    integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
 
 This line creates the integrator to use for advancing the equations of motion.
-It specifies a :class:`LangevinIntegrator`, which performs Langevin dynamics,
+It specifies a :class:`BAOABLangevinIntegrator`, 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
@@ -295,7 +295,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 = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
+        integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
         simulation = Simulation(prmtop.topology, system, integrator)
         simulation.context.setPositions(inpcrd.positions)
         if inpcrd.boxVectors is not None:
@@ -389,7 +389,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 = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
+        integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
         simulation = Simulation(top.topology, system, integrator)
         simulation.context.setPositions(gro.positions)
         simulation.minimizeEnergy()
@@ -453,7 +453,7 @@ on the :class:`CharmmPsfFile`.
         params = CharmmParameterSet('charmm22.rtf', 'charmm22.prm')
         system = psf.createSystem(params, nonbondedMethod=NoCutoff,
                 nonbondedCutoff=1*nanometer, constraints=HBonds)
-        integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
+        integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
         simulation = Simulation(psf.topology, system, integrator)
         simulation.context.setPositions(pdb.positions)
         simulation.minimizeEnergy()
@@ -1022,13 +1022,13 @@ 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.
 
-Langevin Integrator
--------------------
+BAOAB Langevin Integrator
+-------------------------
 
 In the examples of the previous sections, we used Langevin integration:
 ::
 
-    integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
+    integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*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.002 ps).  You
@@ -1037,6 +1037,16 @@ on all values.  For example, the step size could be written either as
 :code:`0.002*picoseconds` or :code:`2*femtoseconds`\ .  They are exactly
 equivalent.
 
+Langevin Integrator
+-------------------
+
+:code:`LangevinIntegrator` is very similar to :code:`BAOABLangevinIntegrator`,
+but it uses a different discretization of the Langevin equation.
+:code:`BAOABLangevinIntegrator` 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
+integrator, so the velocities are offset by half a time step from the positions.
+
 Leapfrog Verlet Integrator
 --------------------------
 
@@ -1145,7 +1155,7 @@ previous section:
     system = prmtop.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
             constraints=HBonds)
     system.addForce(MonteCarloBarostat(1*bar, 300*kelvin))
-    integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
+    integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
     ...
 
 The parameters of the Monte Carlo barostat are the pressure (1 bar) and
@@ -1705,7 +1715,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:`LangevinIntegrator` and then
+time through, it adjusts the temperature of the :class:`BAOABLangevinIntegrator` and then
 calls :code:`step(1000)` to take 1000 time steps.
 
 Applying an External Force to Particles: a Spherical Container
@@ -1737,7 +1747,7 @@ coordinates.  Here is the code to do it:
         system.addForce(force)
         for i in range(system.getNumParticles()):
             force.addParticle(i, [])
-        integrator = LangevinIntegrator(300*kelvin, 91/picosecond, 0.002*picoseconds)
+        integrator = BAOABLangevinIntegrator(300*kelvin, 91/picosecond, 0.002*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
-    LangevinIntegrator integrator(temperature, friction, stepSize);
+    BAOABLangevinIntegrator 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 LangevinIntegrator, initialize a
+System and the Forces.  We then create a BAOABLangevinIntegrator, initialize a
 Context in which to run a simulation, and instruct the Integrator to advance the
 simulation for 10,000 time steps.
 

docs-source/usersguide/references.bib

@@ -240,6 +240,19 @@
    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},

docs-source/usersguide/theory.rst

@@ -1338,6 +1338,21 @@ The integration is done using a leap-frog method similar to VerletIntegrator.
 :cite:`Izaguirre2010` The same comments about the offset between positions and
 velocities apply to this integrator as to that one.
 
+BAOABLangevinIntegrator
+***********************
+
+This integrator is similar to LangevinIntegerator, but it instead uses the BAOAB
+discretization. :cite:`Leimkuhler2013` 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
+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.
+
+Unlike LangevinIntegrator, this does not use a leap-frog algorithm.  The
+positions and velocities all correspond to the same point in time.
+
 BrownianIntegrator
 ******************
 

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 = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
+integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
 simulation = Simulation(prmtop.topology, system, integrator)
 simulation.context.setPositions(inpcrd.positions)
 if inpcrd.boxVectors is not None:

examples/simulateCharmm.py

@@ -22,7 +22,7 @@ params = CharmmParameterSet('charmm22.rtf', 'charmm22.par')
 # Instantiate the system
 system = psf.createSystem(params, nonbondedMethod=NoCutoff,
                           nonbondedCutoff=None)
-integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
+integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*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 = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
+integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*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 = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
+integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
 simulation = Simulation(pdb.topology, system, integrator)
 simulation.context.setPositions(pdb.positions)
 simulation.minimizeEnergy()