Trotter decomposition NPT conserved energy is incorrect
The MTTK (NPT Trotter decomposition) conserved energy is incorrect. My current understanding is that they were initially implemented correctly, but bugs have crept in with subsequent changes. Both bugs only affect the conserved energy and do not influence the dynamics.
Bug 1: MTTK lacks system NHC integral
MTTK uses separate Nose-Hoover chains to couple the different temperature groups (system chains), plus one to couple the barostat. The integral of the system chains is only calculated in the Trotter NVT case, but not in the NPT case. The system chains integral is only calculated if inputrecNvtTrotter(ir)
. This function used to be a macro called IR_NVT_TROTTER(ir)
. The bug was introduced in c7a82654 , when the macros IR_NPT_TROTTER(ir)
and IR_NVT_TROTTER(ir)
in legacyheaders/types/inputrec.h
were changed from meaning "has Trotter pressure coupling" and "has Trotter temperature coupling", respectively, to "is Trotter NPT" and "is Trotter NVT". In other words, IR_NVT_TROTTER(ir)
was true in the NPT case before the change, but not anymore after.
Versions affected: 4.6 and later
Bug 2: MTTK NHC integral formula
The contribution to the conserved energy from the momenta of Nose-Hoover chain degrees of freedom is \sum_{k}\frac{p^{2}_{\eta k}}{2 Q_{k}} = \sum_{k}\frac{1}{2}Q_{k}v^{2}_{\eta k}
(see, e.g., eq. (5.4) in Tuckerman et al 2006, http://dx.doi.org/10.1088/0305-4470/39/19/S18). In GROMACS, the single thermostat contributions are written as 0.5*gmx::square(ivxi[j])/iQinv[j]
. In f2854f89, this was refactored to 0.5*gmx::square(state->nhpres_vxi[i*nh + j]/iQinv)
for the barostat chain only, erroneously squaring the masses.
Versions affected: 2018 and later