Fix MTTK constant energy bugs

Fixes #3796

This also fixes an unclear comment relating to the barostat
Nose-Hoover chain variables in t_state.

diff --git a/docs/release-notes/2020/2020.5.rst b/docs/release-notes/2020/2020.5.rst
index a20bef7710eebdc67a13aa3c0e9832d680f4f6fc..76caa756f75f80fad4ee147a2d1fa70a03fa20cd 100644 (file)
--- a/docs/release-notes/2020/2020.5.rst
+++ b/docs/release-notes/2020/2020.5.rst
@@ -26,6 +26,21 @@ results.
 
 :issue:`3750`
 
+Fixed conserved energy for MTTK
+"""""""""""""""""""""""""""""""
+
+When using `pcoupl=MTTK` and `tcoupl=nose-hoover`, the calculated conserved
+energy was incorrect due to two errors dating back to GROMACS 4.6 and 2018,
+respectively. As a result, all reported conserved energies using this
+combination of temperature and pressure coupling algorithms in any GROMACS
+version since GROMACS 4.6 are likely to be wrong. Note that these errors did
+not impact the dynamics, as the conserved energy is only reported, but never
+used in calculations. Also note that this bug only affects this exact
+combination of temperature / pressure coupling algorithms.
+
+:issue:`3796`
+
+
 Fixes for ``gmx`` tools
 ^^^^^^^^^^^^^^^^^^^^^^^
 

diff --git a/src/gromacs/mdlib/coupling.cpp b/src/gromacs/mdlib/coupling.cpp
index 909564f2af96958b0b15f0b413864324ba70aed9..bdf47ab3e84324b059e8a586b94e5f2ae9bf3ac4 100644 (file)
--- a/src/gromacs/mdlib/coupling.cpp
+++ b/src/gromacs/mdlib/coupling.cpp
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,2018,2019,2020, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -1312,7 +1312,7 @@ static real energyNoseHoover(const t_inputrec* ir, const t_state* state, const t
 
         if (nd > 0.0)
         {
-            if (inputrecNvtTrotter(ir))
+            if (inputrecNvtTrotter(ir) || inputrecNptTrotter(ir))
             {
                 /* contribution from the thermal momenta of the NH chain */
                 for (int j = 0; j < nh; j++)
@@ -1363,7 +1363,7 @@ static real energyPressureMTTK(const t_inputrec* ir, const t_state* state, const
             double iQinv = MassQ->QPinv[i * nh + j];
             if (iQinv > 0)
             {
-                energy += 0.5 * gmx::square(state->nhpres_vxi[i * nh + j] / iQinv);
+                energy += 0.5 * gmx::square(state->nhpres_vxi[i * nh + j]) / iQinv;
                 /* contribution from the thermal variable of the NH chain */
                 energy += state->nhpres_xi[i * nh + j] * kT;
             }

diff --git a/src/gromacs/mdtypes/state.h b/src/gromacs/mdtypes/state.h
index e3a3ce5c6bd0e1923f899ef41db8d891d65109bf..a54bff29bbed91811711587af795d3195ea34fcc 100644 (file)
--- a/src/gromacs/mdtypes/state.h
+++ b/src/gromacs/mdtypes/state.h
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,2018,2019,2020, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -222,8 +222,8 @@ public:
     // All things public
     int natoms; //!< Number of atoms, local + non-local; this is the size of \p x, \p v and \p cg_p, when used
     int ngtc;          //!< The number of temperature coupling groups
-    int nnhpres;       //!< The NH-chain length for the MTTK barostat
-    int nhchainlength; //!< The NH-chain length for temperature coupling
+    int nnhpres;       //!< The number of NH-chains for the MTTK barostat (always 1 or 0)
+    int nhchainlength; //!< The NH-chain length for temperature coupling and MTTK barostat
     int flags; //!< Set of bit-flags telling which entries are present, see enum at the top of the file
     int                      fep_state;      //!< indicates which of the alchemical states we are in
     std::array<real, efptNR> lambda;         //!< Free-energy lambda vector