From: Pascal Merz <pascal.merz@me.com>
Date: Wed, 2 Dec 2020 07:37:43 +0000 (-0700)
Subject: Fix kinetic energy and temperatures reporting
X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=8e27c9b954709dafc870fbe594fa7a8ea5d08d61;p=alexxy%2Fgromacs.git

Fix kinetic energy and temperatures reporting

Fixes the reporting of kinetic energy and temperatures, which is
(very slightly) off in md-vv with Trotter NPT (Nose-Hoover & MTTK)
in the legacy code path.

Closes #3832
---

diff --git a/docs/release-notes/2020/2020.5.rst b/docs/release-notes/2020/2020.5.rst
index a615bbeb06..50b599f10a 100644
--- a/docs/release-notes/2020/2020.5.rst
+++ b/docs/release-notes/2020/2020.5.rst
@@ -51,6 +51,21 @@ combination of temperature / pressure coupling algorithms.
 
 :issue:`3796`
 
+Fixed kinetic energy and temperature reporting for MTTK
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+When using `pcoupl=MTTK` and `tcoupl=nose-hoover`, the reported kinetic
+energy and temperature were very slightly off. The integration of the
+temperature coupling trailed the reporting by half a time step. Note that
+these errors did not impact the dynamics, as the quantities were correctly
+integrated and only wrongly reported. Also note that the difference is so
+small that it is unlikely to have been significant for any application
+except for rigorous algorithm validation. Finally, note that this bug
+only affects this exact combination of temperature / pressure coupling
+algorithms.
+
+:issue:`3832`
+
 Fix pull error message with angles and dihedrals
 """"""""""""""""""""""""""""""""""""""""""""""""
 
diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp
index 9e519d44da..7df4a68b21 100644
--- a/src/gromacs/mdrun/md.cpp
+++ b/src/gromacs/mdrun/md.cpp
@@ -1051,7 +1051,7 @@ void gmx::LegacySimulator::do_md()
                         copy_mat(shake_vir, state->svir_prev);
                         copy_mat(force_vir, state->fvir_prev);
                     }
-                    if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
+                    if ((inputrecNptTrotter(ir) || inputrecNvtTrotter(ir)) && ir->eI == eiVV)
                     {
                         /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
                         enerd->term[F_TEMP] =