Add more info to NH integrator docs

docs-source/usersguide/application.rst

@@ -1068,18 +1068,54 @@ temperature control with a velocity scaling algorithm that produces more
 accurate transport properties :cite:`Basconi2013`.  This velocity scaling
 results from propagating a chain of extra variables, which slightly reduces the
 computational efficiency with respect to :code:`LangevinMiddleIntegrator`.  The
-thermostated integrator is minimally created similarly to the
+thermostated integrator is minimally created with syntax analogous to the
 :code:`LangevinMiddleIntegrator` example above::
 
-    NoseHooverIntegrator integrator(300*kelvin, 25/picosecond, 0.004*picoseconds);
+    NoseHooverIntegrator integrator(300*kelvin, 25/picosecond,
+                                    0.004*picoseconds);
 
 The first argument specifies the target temperature.  The second specifies the
 frequency of interaction with the heat bath: a lower value interacts minimally,
 yielding the microcanonical ensemble in the limit of a zero frequency, while a
 larger frequency will perturb the system greater, keeping it closer to the
 target temperature.  The third argument is the integration timestep that, like
-the other arguments, must be specified with units.  Note that for this
-integrator, the velocities are offset by half a time step from the positions.
+the other arguments, must be specified with units.
+
+More verbose invocation is possible::
+
+    NoseHooverIntegrator integrator(300*kelvin, 25/picosecond,
+                                    0.004*picoseconds, 3, 3, 7);
+
+where the extra arguments are the chain length, number of multi-timesteps and
+number of Yoshida-Suzuki points, as described in section
+:ref:`nosehoover-integrators-theory`.
+
+This integrator supports partitioning of the system into an arbitrary number of
+subsystems, where each can be coupled to a different heat bath.  As an example::
+
+    NoseHooverIntegrator integrator(0.004*picoseconds)
+    integrator.addSubsystemThermostat(thermostatedParticles, thermostatedPairs,
+                                      centerOfMassTemperature, centerOfMassCollisionFrequency,
+                                      relativeTemperature, relativeCollisionFrequency,
+                                      chainLength, numMTS, numYoshidaSuzuki)
+
+Here the first argument to :code:`addSubsystemThermostat` is a list of the
+particles to be controlled by that thermostat, which will be thermostated to
+the :code:`centerOfMassTemperature`.  The :code:`thermostatedPairs` is a list
+of pairs whose center of mass motion will be thermostated to
+:code:`centerOfMassTemperature`, while the relative motion of each pair will be
+thermostated to :code:`relativeTemperature`.  It is possible for one of
+:code:`thermostatedAtoms` or :code:`thermostatedPairs` to be empty.  The final
+three arguments have the same default values as the simple constructor
+described above.  Making multiple calls to :code:`addSubsystemThermostat` with
+different lists of atoms and pairs permits the creation of differently
+thermostated subsystems: note that because some platforms can reorder similar
+molecules, all molecules of a given type must be controled by the same
+thermostat.  This more advanced syntax allows the setup of Drude oscillator
+simulations, but a more convenient form is provided through the
+:code:`DrudeNoseHooverIntegrator` class, described in :ref:`drude-plugin`.
+Note that for this integrator, the velocities are offset by half a time step
+from the positions.
 
 Leapfrog Verlet Integrator
 --------------------------

docs-source/usersguide/library.rst

@@ -3654,4 +3654,15 @@ The equations of motion can be integrated with three different methods:
    the target temperature with one thermostat.  Another thermostat is used to keep
    relative motion of Drude pairs to a different, typically much lower,
    temperature to maintain separation of nuclear and electronic degrees of
-   freedom.
+   freedom.  The minimal specification is as follows::
+
+      DrudeNoseHooverIntegrator integrator(temperature, frequency,
+                                           temperatureDrude, frequencyDrude,
+                                           1*femtoseconds)
+
+   Where the first and third arguments specify the center-of-mass temperature and
+   relative temperature for each Drude pair, respecitvely.  The second and fourth
+   arguments describe the frequency of interaction with the center-of-mass and
+   relative heat baths, respectively, and should be specified with inverse time
+   units.  The fifth argument is the timestep.  The multi-timestep and Nosé-Hoover
+   chain length may also be specified, but sensible defaults are provided.

docs-source/usersguide/references.bib

@@ -325,6 +325,16 @@
     journal = {Europhysics Letters ({EPL})},
 }
 
+@article{Martyna1992,
+   author = {Glenn J. Martyna and Michael L. Klein and Mark Tuckerman},
+   year = 1992,
+   title = {Nos\'{e}–Hoover chains: The canonical ensemble via continuous dynamics},
+   pages = {2635-2643},
+   journal  = {Journal of Chemical Physics},
+   volume = {97},
+   issue = {4},
+}
+
 @article{Markland2008
    author = {Markland, Thomas E. and Manolopoulos, David E.},
    title = {An efficient ring polymer contraction scheme for imaginary time path integral simulations},

docs-source/usersguide/theory.rst

@@ -1391,6 +1391,52 @@ 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.
 
+.. _nosehoover-integrators-theory:
+
+NoseHooverIntegrator
+********************
+
+Like LangevinMiddleIntegerator, this 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-\Delta t/2) + \mathbf{f}_{i}(t)\Delta t/{m}_{i}
+
+
+.. math::
+   \mathbf{r}_{i}(t+\Delta t/2) = \mathbf{r}_{i}(t) + \mathbf{v}_{i}(t+\Delta t/2)\Delta t/2
+
+
+.. math::
+   \mathbf{v'}_{i}(t+\Delta t/2) = \mathrm{scale}\times\mathbf{v}_{i}(t+\Delta t/2)
+
+
+.. math::
+   \mathbf{r}_{i}(t+\Delta t) = \mathbf{r}_{i}(t+\Delta t/2) + \mathbf{v'}_{i}(t+\Delta t/2)\Delta t/2
+
+
+The universal scale factor used in the third step is determined by propagating
+auxilliary degrees of freedom alongside the regular particles.  The original
+Nosé-Hoover formulation used a single harmonic oscillator for the heat bath,
+but this is problematic in small or stiff systems, which are non-ergodic, so
+the chain formulation extends this by replacing the single oscillator
+thermostat with a chain of connected oscillators.  :cite:`Martyna1992`  For
+large systems a single oscillator (*i.e.* a chain length of one) will suffice,
+but longer chains are necessary to properly thermostat non-ergodic systems.
+The OpenMM default is to use a chain length of three to cover the latter case,
+but this can be safely reduced to increase efficiency in large systems.
+
+The heat bath propagation is performed using a multi-timestep algorithm.  Each
+propagation step is discretized into substeps using a factorization from
+Yoshida and Suzuki; the default discretization uses a :math:`\mathcal{O}(\Delta
+t^6)` approach that uses 7 points, but 1, 3 or 5 points may also be used to
+increase performace, at the expense of accuracy.  Each step is further
+subdivided into multi-timesteps with a default of 3 multi time steps per
+propagation; as with the number of Yoshida-Suziki points this value may be
+increase to increase accuracy but with additional computational expense.
+
 BrownianIntegrator
 ******************