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 a615bbeb060c267aa7e9632acff767b1dba52581..50b599f10ac026e2ce8c6249574472bb64896472 100644 (file)
--- 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 9e519d44daffc17591d92f82307a15ab93159128..7df4a68b2186a08b74a08f54434ed547547b1785 100644 (file)
--- 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] =