## 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