DRR: correctly separate the calculate and update (#1132)

The previous implementation updates the internal states partially in
calculate(), which should be done in update(). This commit fixes the
issue.

src/drr/DRR.cpp

@@ -414,35 +414,48 @@ void DRRForceGrid::writeDivergence(const string &filename, const string &fmt) co
   fclose(pDiv);
 }
 
-bool ABF::store_getbias(const vector<double> &pos, const vector<double> &f,
-                        vector<double> &fbias) {
+bool ABF::getbias(const vector<double> &pos, const vector<double> &f,
+                  vector<double> &fbias) const {
   if (!isInBoundary(pos)) {
     std::fill(begin(fbias), end(fbias), 0.0);
     return false;
   }
   const size_t baseaddr = sampleAddress(pos);
-  unsigned long int &count = samples[baseaddr];
-  ++count;
+  const unsigned long int count = samples[baseaddr] + 1;
   double factor = 2 * (static_cast<double>(count)) / mFullSamples - 1;
   auto it_fa = begin(forces) + baseaddr * ndims;
   auto it_fb = begin(fbias);
   auto it_f = begin(f);
   auto it_maxFactor = begin(mMaxFactors);
-  do {
+  while (it_f != end(f)) {
     // Clamp to [0,maxFactors]
     factor = factor < 0 ? 0 : factor > (*it_maxFactor) ? (*it_maxFactor) : factor;
-    (*it_fa) += (*it_f); // Accumulate instantaneous force
-    (*it_fb) = factor * (*it_fa) * (-1.0) /
+    (*it_fb) = factor * ((*it_fa) + (*it_f)) * (-1.0) /
                static_cast<double>(count); // Calculate bias force
     ++it_fa;
     ++it_fb;
     ++it_f;
     ++it_maxFactor;
-  } while (it_f != end(f));
-
+  };
   return true;
 }
 
+void ABF::store_force(const vector<double> &pos,
+                      const vector<double> &f) {
+  if (!isInBoundary(pos)) {
+    return;
+  }
+  const size_t baseaddr = sampleAddress(pos);
+  samples[baseaddr] += 1;
+  auto it_fa = begin(forces) + baseaddr * ndims;
+  auto it_f = begin(f);
+  while (it_f != end(f)) {
+    *it_fa += *it_f;
+    ++it_fa;
+    ++it_f;
+  };
+}
+
 ABF ABF::mergewindow(const ABF &aWA, const ABF &aWB) {
   const vector<DRRAxis> dR = merge(aWA.dimensions, aWB.dimensions);
   const string suffix = ".abf";

src/drr/DRR.h

@@ -310,9 +310,11 @@ class ABF : public DRRForceGrid {
     mMaxFactors = maxFactors;
   }
   // Store the "instantaneous" spring force of a point and get ABF bias forces.
-  bool store_getbias(const vector<double> &pos,
-                     const vector<double> &f,
-                     vector<double> &fbias);
+  bool getbias(const vector<double> &pos,
+               const vector<double> &f,
+               vector<double> &fbias) const;
+  void store_force(const vector<double> &pos,
+                   const vector<double> &f);
   static ABF mergewindow(const ABF &aWA, const ABF &aWB);
   ~ABF() {}
 

src/drr/DynamicReferenceRestraining.cpp

@@ -701,7 +701,6 @@ DynamicReferenceRestraining::DynamicReferenceRestraining(
 }
 
 void DynamicReferenceRestraining::calculate() {
-  long long int step_now = getStep();
   if (firsttime) {
     for (size_t i = 0; i < ndims; ++i) {
       fict[i] = getArgument(i);
@@ -711,6 +710,46 @@ void DynamicReferenceRestraining::calculate() {
     }
     firsttime = false;
   }
+  if (withExternalForce == false) {
+    double ene = 0.0;
+    for (size_t i = 0; i < ndims; ++i) {
+      real[i] = getArgument(i);
+      springlength[i] = difference(i, fict[i], real[i]);
+      fictNoPBC[i] = real[i] - springlength[i];
+      double f = -kappa[i] * springlength[i];
+      ffict_measured[i] = -f;
+      ene += 0.5 * kappa[i] * springlength[i] * springlength[i];
+      setOutputForce(i, f);
+      ffict[i] = -f;
+      fict[i] = fictValue[i]->bringBackInPbc(fict[i]);
+      fictValue[i]->set(fict[i]);
+      vfictValue[i]->set(vfict_laststep[i]);
+      springforceValue[i]->set(ffict_measured[i]);
+      fictNoPBCValue[i]->set(fictNoPBC[i]);
+    }
+    setBias(ene);
+    ABFGrid.getbias(fict, ffict_measured, fbias);
+  } else {
+    for (size_t i = 0; i < ndims; ++i) {
+      real[i] = getArgument(i);
+      ffict_measured[i] = externalForceValue[i]->get();
+      if (withExternalFict) {
+        fictNoPBC[i] = externalFictValue[i]->get();
+      }
+      springforceValue[i]->set(ffict_measured[i]);
+      fictNoPBCValue[i]->set(fictNoPBC[i]);
+    }
+    ABFGrid.getbias(real, ffict_measured, fbias);
+    if (!nobias) {
+      for (size_t i = 0; i < ndims; ++i) {
+        setOutputForce(i, fbias[i]);
+      }
+    }
+  }
+}
+
+void DynamicReferenceRestraining::update() {
+  const long long int step_now = getStep();
   if (step_now != 0) {
     if ((step_now % int(outputfreq)) == 0) {
       save(outputname, step_now);
@@ -748,40 +787,18 @@ void DynamicReferenceRestraining::calculate() {
     }
   }
   if (withExternalForce == false) {
-    double ene = 0.0;
     for (size_t i = 0; i < ndims; ++i) {
       real[i] = getArgument(i);
       springlength[i] = difference(i, fict[i], real[i]);
-      fictNoPBC[i] = real[i] - springlength[i];
-      double f = -kappa[i] * springlength[i];
-      ffict_measured[i] = -f;
-      ene += 0.5 * kappa[i] * springlength[i] * springlength[i];
-      setOutputForce(i, f);
-      ffict[i] = -f;
-      fict[i] = fictValue[i]->bringBackInPbc(fict[i]);
-      fictValue[i]->set(fict[i]);
-      vfictValue[i]->set(vfict_laststep[i]);
-      springforceValue[i]->set(ffict_measured[i]);
-      fictNoPBCValue[i]->set(fictNoPBC[i]);
+      ffict_measured[i] = kappa[i] * springlength[i];
     }
-    setBias(ene);
-    ABFGrid.store_getbias(fict, ffict_measured, fbias);
+    ABFGrid.store_force(fict, ffict_measured);
   } else {
     for (size_t i = 0; i < ndims; ++i) {
       real[i] = getArgument(i);
       ffict_measured[i] = externalForceValue[i]->get();
-      if (withExternalFict) {
-        fictNoPBC[i] = externalFictValue[i]->get();
-      }
-      springforceValue[i]->set(ffict_measured[i]);
-      fictNoPBCValue[i]->set(fictNoPBC[i]);
-    }
-    ABFGrid.store_getbias(real, ffict_measured, fbias);
-    if (!nobias) {
-      for (size_t i = 0; i < ndims; ++i) {
-        setOutputForce(i, fbias[i]);
-      }
     }
+    ABFGrid.store_force(real, ffict_measured);
   }
   if (useCZARestimator) {
     CZARestimator.store(real, ffict_measured);
@@ -790,9 +807,6 @@ void DynamicReferenceRestraining::calculate() {
     eabf_UI.update_output_filename(outputprefix);
     eabf_UI.update(int(step_now), real, fictNoPBC);
   }
-}
-
-void DynamicReferenceRestraining::update() {
   if (withExternalForce == false) {
     for (size_t i = 0; i < ndims; ++i) {
       // consider additional forces on the fictitious particle