Fix handling of previous virials

These quantities get written to checkpoint files only for the Trotter
pressure-coupling integrators that need them, but they were being
copied in do_md for all Trotter integrators. This meant that an
appending restart of md-vv plus nose-hoover plus no pressure coupling
truncated off a correct edr frame and wrote one with zero virial and
wrong pressure. And in the same case, a no-append restart writes a
duplicate frame that does not agree with the one written before
termination.

Fixed by copying them only when they are useful to copy, so that in
the problem case, the correctly computed post-restart virial is not
wiped with zeroes from state fields that were never read from the
checkpoint. Cases that use the previous-virial values are not
affected.

Refs #1793

Change-Id: I84908122aefbe8658f423eaf4e5bd4ae25a93d24

diff --git a/src/programs/mdrun/md.cpp b/src/programs/mdrun/md.cpp
index afa3d62a75a4e225841c05283444184851e8fa20..340d692c2abceebe83ad9c9c67ccc5814fba1c31 100644 (file)
--- a/src/programs/mdrun/md.cpp
+++ b/src/programs/mdrun/md.cpp
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2011,2012,2013,2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2011,2012,2013,2014,2015,2016,2017, 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.
@@ -1215,8 +1215,15 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             }
             if (bTrotter && !bInitStep)
             {
-                copy_mat(shake_vir, state->svir_prev);
-                copy_mat(force_vir, state->fvir_prev);
+                /* TODO This is only needed when we're about to write
+                 * a checkpoint, because we use it after the restart
+                 * (in a kludge?). But what should we be doing if
+                 * startingFromCheckpoint or bInitStep are true? */
+                if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
+                {
+                    copy_mat(shake_vir, state->svir_prev);
+                    copy_mat(force_vir, state->fvir_prev);
+                }
                 if (IR_NVT_TROTTER(ir) && ir->eI == eiVV)
                 {
                     /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
@@ -1280,7 +1287,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, state->x, ir, t, wcycle);
 
         /* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
-        if (bStartingFromCpt && bTrotter)
+        if (bStartingFromCpt && (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir)))
         {
             copy_mat(state->svir_prev, shake_vir);
             copy_mat(state->fvir_prev, force_vir);