From 81259f9110c9c92ffd594be2c7e5695ac339f059 Mon Sep 17 00:00:00 2001 From: HanatoK Date: Thu, 10 Oct 2024 11:09:53 -0500 Subject: [PATCH] DRR: correctly separate the calculate and update 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 | 31 +++++++---- src/drr/DRR.h | 8 +-- src/drr/DynamicReferenceRestraining.cpp | 72 +++++++++++++++---------- 3 files changed, 70 insertions(+), 41 deletions(-) diff --git a/src/drr/DRR.cpp b/src/drr/DRR.cpp index 88105a6c6b..bfb115dbe0 100644 --- a/src/drr/DRR.cpp +++ b/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 &pos, const vector &f, - vector &fbias) { +bool ABF::getbias(const vector &pos, const vector &f, + vector &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(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(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 &pos, + const vector &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 dR = merge(aWA.dimensions, aWB.dimensions); const string suffix = ".abf"; diff --git a/src/drr/DRR.h b/src/drr/DRR.h index 3f6fd88802..9dc63c5712 100644 --- a/src/drr/DRR.h +++ b/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 &pos, - const vector &f, - vector &fbias); + bool getbias(const vector &pos, + const vector &f, + vector &fbias) const; + void store_force(const vector &pos, + const vector &f); static ABF mergewindow(const ABF &aWA, const ABF &aWB); ~ABF() {} diff --git a/src/drr/DynamicReferenceRestraining.cpp b/src/drr/DynamicReferenceRestraining.cpp index c85ece3ece..a4a0fdddd2 100644 --- a/src/drr/DynamicReferenceRestraining.cpp +++ b/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