Nosé--Hoover documentation has several typos and is not consistent with source code
Summary
The description of Nosé--Hoover temperature coupling in the reference manual has several typos and unexplained variables. In addition, the equations presented in the manual are inconsistent with the implementation in the source code
Relevant files
src/gromacs/modularsimulator/nosehooverchains.cpp
docs/reference-manual/algorithms/molecular-dynamics.rst
References
- M.P. Allen and D.J. Tildesley. Computer Simulation of Liquids, Chapter 3.8.2, 2nd Ed. New York: Oxford University Press, 2017.
Use cases
Difficult to reproduce NVT calculations from other software, as the appropriate choice of Nosé--Hoover parameters (specifically, the Q masses) cannot easily be determined
Impact
Fix requires relatively little work: simply correcting the documentation.
The fix is useful in better understanding the choice of tau-t
for the Nosé--Hoover thermostat.
Detailed description
Typos
In the "Nosé--Hoover temperature coupling" section of the reference manual, there are several missing factors of the Boltzmann constant k
, as well as confusions between the reference temperature T_0
and instantaneous temperature T
. These include:
- Eq. (5.43): missing factor of
k
on the right-hand side - Eq. (5.45): missing factor of
k
in the numerator of the right-hand side - equations above (5.46):
- missing
k
multiplying(T - T_0)
in the second equation - subscript of
\xi_i
should go from2 \ldots M-1
in the third equation - subscript
\xi_N
should be\xi_M
in the third equation - both the third and fourth equations should have
T_0
rather thanT
- missing
- missing factor of
k
in the equation below (5.47)
Additional typos:
- Eqs. (5.44) and (5.46) should have
\boldsymbol{p}_i \boldsymbol{\cdot} \boldsymbol{p}_i
on the right-hand side, for the kinetic energy - the equation for
\xi
is not given in the single-chain case, and the equation for\xi_i
is not given in the chained case - the quantity
Q'_k
is not defined in Eq. (5.46) - the quantities
Q_k
are not related to the period\tau_T
in the chained case, in an analogous way to Eq. (5.45)
Major inconsistency:
Given the definition of p_{\xi_i}
and Q_i
in the manual, the driving force in the last two equations above (5.46) should be k T_0 / N_f
, not k T_0
. This results from the manual's choice to scale both p_{\xi_i}
and Q_i
by a factor of N_f
relative to Allen and Tildesley [1]. This can easily be seen from the governing equation for p_{\xi_1}
, and comparing it to Eq. (3.68) in Ref. [1]. According to this scaling, the definition of Q'_k
highlighted above would be given by Q'_k := Q_k / N_f
.
Unfortunately, however, the current implementation of the Nosé--Hoover thermostat in the source code does not use this additional scaling by N_f
. This is easily seen in the nosehooverchains.cpp
file mentioned above, in which invXiMass
(i.e. 1/Q_k
) tells us that
Q_k = \frac{\tau_T^2 k T_0}{4 \pi^2} \cdot \Big( N_f, 1, 1, \ldots, 1 \Big)
Thus, the source code in fact uses the definition of the "bath particle" mass and momenta in Allen and Tildesley. Changing the documentation to accurately reflect the source code would simply require adding factors of N_f
in several location. The choice of the Q_k
would then be straightforward to compare with other numerical implementations of the Nosé--Hoover thermostat.