Zero pressure and virial tensors at beginning of step

This sets pressure and virial tensors to zero when they are
accessed for the first time in a step. Not doing this does not
yield wrong behavior currently, but carries the risk of
hard-to-spot bugs moving forward.

Change-Id: I2553bd52fa56a497a61c31738ec4436f7f587ece

diff --git a/src/gromacs/modularsimulator/energyelement.cpp b/src/gromacs/modularsimulator/energyelement.cpp
index a927ac9b3225fee8475e985edc83a18c4c299e16..70d885462b1471b17a318cd33cfe6b3a95c8f222 100644 (file)
--- a/src/gromacs/modularsimulator/energyelement.cpp
+++ b/src/gromacs/modularsimulator/energyelement.cpp
@@ -93,10 +93,8 @@ EnergyElement::EnergyElement(
     logWritingStep_(-1),
     forceVirialStep_(-1),
     shakeVirialStep_(-1),
-#ifndef NDEBUG
     totalVirialStep_(-1),
     pressureStep_(-1),
-#endif
     needToSumEkinhOld_(false),
     startingBehavior_(startingBehavior),
     dummyLegacyState_(),
@@ -319,6 +317,11 @@ void EnergyElement::addToConstraintVirial(const tensor virial, Step step)
 
 rvec* EnergyElement::forceVirial(Step gmx_unused step)
 {
+    if (step > forceVirialStep_)
+    {
+        forceVirialStep_ = step;
+        clear_mat(forceVirial_);
+    }
     GMX_ASSERT(step >= forceVirialStep_ || forceVirialStep_ == -1,
                "Asked for force virial of previous step.");
     return forceVirial_;
@@ -326,6 +329,11 @@ rvec* EnergyElement::forceVirial(Step gmx_unused step)
 
 rvec* EnergyElement::constraintVirial(Step gmx_unused step)
 {
+    if (step > shakeVirialStep_)
+    {
+        shakeVirialStep_ = step;
+        clear_mat(shakeVirial_);
+    }
     GMX_ASSERT(step >= shakeVirialStep_ || shakeVirialStep_ == -1,
                "Asked for constraint virial of previous step.");
     return shakeVirial_;
@@ -333,6 +341,11 @@ rvec* EnergyElement::constraintVirial(Step gmx_unused step)
 
 rvec* EnergyElement::totalVirial(Step gmx_unused step)
 {
+    if (step > totalVirialStep_)
+    {
+        totalVirialStep_ = step;
+        clear_mat(totalVirial_);
+    }
     GMX_ASSERT(step >= totalVirialStep_ || totalVirialStep_ == -1,
                "Asked for total virial of previous step.");
     return totalVirial_;
@@ -340,6 +353,11 @@ rvec* EnergyElement::totalVirial(Step gmx_unused step)
 
 rvec* EnergyElement::pressure(Step gmx_unused step)
 {
+    if (step > pressureStep_)
+    {
+        pressureStep_ = step;
+        clear_mat(pressure_);
+    }
     GMX_ASSERT(step >= pressureStep_ || pressureStep_ == -1,
                "Asked for pressure of previous step.");
     return pressure_;

diff --git a/src/gromacs/modularsimulator/energyelement.h b/src/gromacs/modularsimulator/energyelement.h
index 63cb0e5eb3a694d5540136df4a838d9f0892ec32..8420513b98f8655dfcc7fcbe2e062bb32413eb1b 100644 (file)
--- a/src/gromacs/modularsimulator/energyelement.h
+++ b/src/gromacs/modularsimulator/energyelement.h
@@ -299,12 +299,10 @@ class EnergyElement final :
         Step forceVirialStep_;
         //! The step number of the current constraint virial tensor
         Step shakeVirialStep_;
-#ifndef NDEBUG
         //! The step number of the current total virial tensor
         Step totalVirialStep_;
         //! The step number of the current pressure tensor
         Step pressureStep_;
-#endif
 
         //! Whether ekinh_old needs to be summed up (set by compute globals)
         bool needToSumEkinhOld_;