Merge pull request #2621 from andysim/nhdocs

Nosé Hoover Integrator docs

docs-source/usersguide/application.rst

@@ -1059,6 +1059,34 @@ sampling, and therefore is preferred for most applications.  Also note that
 :code:`LangevinIntegrator`\ , like :code:`LangevinMiddleIntegrator`\ , is a leapfrog
 integrator, so the velocities are offset by half a time step from the positions.
 
+Nosé-Hoover Integrator
+----------------------
+
+The :code:`NoseHooverIntegrator` uses the same "middle" leapfrog propagation
+algorithm as :code:`LangevinMiddleIntegrator`, but replaces the stochastic
+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 with syntax analogous to the
+:code:`LangevinMiddleIntegrator` example above::
+
+    NoseHooverIntegrator integrator(300*kelvin, 1/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.  For initial equilibration
+to the target temperature, a larger interaction frequency is recommended,
+*e.g.* 25 ps\ :sup:`-1`.
+
+This integrator supports lots of other options, including the ability to couple
+different parts of the system to thermostats at different temperatures. See the
+API documentation for details.
+
 Leapfrog Verlet Integrator
 --------------------------
 

docs-source/usersguide/library.rst

@@ -3634,7 +3634,7 @@ the Drude particle and the spring constant *k* by
 A damped interaction\ :cite:`Thole1981` is used between dipoles that are
 bonded to each other.
 
-The equations of motion can be integrated with two different methods:
+The equations of motion can be integrated with three different methods:
 
 #. In the Self Consistent Field (SCF) method, the ordinary particles are first
    updated as usual.  A local energy minimization is then performed to select new
@@ -3649,3 +3649,20 @@ The equations of motion can be integrated with two different methods:
    temperature, while using a much lower temperature for their relative internal
    motion.  In practice, this produces dipole moments very close to those from the
    SCF solution while being much faster to compute.
+#. The Nosé-Hoover dual thermostat method.  In this approach the motion of
+   non-Drude sites and center of mass motion of Drude pairs are thermostated to
+   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.  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

@@ -41,6 +41,17 @@
   doi = {10.1103/PhysRevLett.100.020603},
 }
 
+@article{Basconi2013,
+  title = {Effects of Temperature Control Algorithms on Transport Properties and Kinetics in Molecular Dynamics Simulations},
+  author = {Joseph E. Bascon and Michael R. Shirts},
+  journal = {Journal of Chemical Theory and Computation},
+  volume = {9},
+  issue = {7},
+  pages = {2887-2899},
+  year = {2013},
+  doi= {10.1021/ct400109a}
+}
+
 @article{Berendsen1987
    author = {Berendsen, H. J. C. and Grigera, J. R. and Straatsma, T. P.},
    title = {The missing term in effective pair potentials},
@@ -314,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
 ******************