Prevent extra force evaluations with NoseHooverIntegrator (#4753)

olla/include/openmm/kernels.h

@@ -1124,10 +1124,8 @@ public:
      * 
      * @param context    the context in which to execute this kernel
      * @param integrator the NoseHooverIntegrator this kernel is being used for
-     * @param forcesAreValid a reference to the parent integrator's boolean for keeping
-     *                       track of the validity of the current forces.
      */
-    virtual void execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid) = 0;
+    virtual void execute(ContextImpl& context, const NoseHooverIntegrator& integrator) = 0;
     /**
      * Compute the kinetic energy.
      * 

openmmapi/src/NoseHooverIntegrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2019-2023 Stanford University and the Authors.      *
+ * Portions copyright (c) 2019-2024 Stanford University and the Authors.      *
  * Authors: Andreas Krämer and Andrew C. Simmonett                            *
  * Contributors: Peter Eastman                                                *
  *                                                                            *
@@ -47,19 +47,13 @@ using namespace OpenMM;
 using std::string;
 using std::vector;
 
-NoseHooverIntegrator::NoseHooverIntegrator(double stepSize):
-    forcesAreValid(false),
-    hasSubsystemThermostats_(true)
-{
+NoseHooverIntegrator::NoseHooverIntegrator(double stepSize) : hasSubsystemThermostats_(true) {
     setStepSize(stepSize);
     setConstraintTolerance(1e-5);
     setMaximumPairDistance(0.0);
 }
 NoseHooverIntegrator::NoseHooverIntegrator(double temperature, double collisionFrequency, double stepSize,
-                                           int chainLength, int numMTS, int numYoshidaSuzuki) : 
-    forcesAreValid(false),
-    hasSubsystemThermostats_(false) {
-
+                                           int chainLength, int numMTS, int numYoshidaSuzuki) : hasSubsystemThermostats_(false) {
     setStepSize(stepSize);
     setConstraintTolerance(1e-5);
     setMaximumPairDistance(0.0);
@@ -294,7 +288,6 @@ void NoseHooverIntegrator::setRelativeCollisionFrequency(double frequency, int c
 }
 
 double NoseHooverIntegrator::computeKineticEnergy() {
-    forcesAreValid = false;
     double kE = 0.0;
     if(noseHooverChains.size() > 0) {
         for (const auto &nhc: noseHooverChains){
@@ -329,7 +322,6 @@ void NoseHooverIntegrator::initialize(ContextImpl& contextRef) {
     owner = &contextRef.getOwner();
     kernel = context->getPlatform().createKernel(IntegrateNoseHooverStepKernel::Name(), contextRef);
     kernel.getAs<IntegrateNoseHooverStepKernel>().initialize(contextRef.getSystem(), *this);
-    forcesAreValid = false;
 
     // check for drude particles and build the Nose-Hoover Chains
     for (auto& thermostat: noseHooverChains){
@@ -367,10 +359,8 @@ void NoseHooverIntegrator::step(int steps) {
     if (context == NULL)
         throw OpenMMException("This Integrator is not bound to a context!");
     for (int i = 0; i < steps; ++i) {
-        if(context->updateContextState())
-            forcesAreValid = false;
         context->calcForcesAndEnergy(true, false, getIntegrationForceGroups());
-        kernel.getAs<IntegrateNoseHooverStepKernel>().execute(*context, *this, forcesAreValid);
+        kernel.getAs<IntegrateNoseHooverStepKernel>().execute(*context, *this);
     }
 }
 

platforms/common/include/openmm/common/CommonKernels.h

@@ -1190,10 +1190,8 @@ public:
      *
      * @param context    the context in which to execute this kernel
      * @param integrator the VerletIntegrator this kernel is being used for
-     * @param forcesAreValid a reference to the parent integrator's boolean for keeping
-     *                       track of the validity of the current forces.
      */
-    void execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid);
+    void execute(ContextImpl& context, const NoseHooverIntegrator& integrator);
     /**
      * Compute the kinetic energy.
      *

platforms/common/src/CommonKernels.cpp

@@ -5895,7 +5895,7 @@ void CommonIntegrateNoseHooverStepKernel::initialize(const System& system, const
         energyBuffer.initialize<mm_float2>(cc, energyBufferSize, "energyBuffer");
 }
 
-void CommonIntegrateNoseHooverStepKernel::execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid) {
+void CommonIntegrateNoseHooverStepKernel::execute(ContextImpl& context, const NoseHooverIntegrator& integrator) {
     ContextSelector selector(cc);
     IntegrationUtilities& integration = cc.getIntegrationUtilities();
     int paddedNumAtoms = cc.getPaddedNumAtoms();
@@ -5904,7 +5904,8 @@ void CommonIntegrateNoseHooverStepKernel::execute(ContextImpl& context, const No
 
     // If the atom reordering has occured, the forces from the previous step are permuted and thus invalid.
     // They need to be either sorted or recomputed; here we choose the latter.
-    if (!forcesAreValid || cc.getAtomsWereReordered()) context.calcForcesAndEnergy(true, false, integrator.getIntegrationForceGroups());
+    if (cc.getAtomsWereReordered())
+        context.calcForcesAndEnergy(true, false, integrator.getIntegrationForceGroups());
 
     const auto& atomList = integrator.getAllThermostatedIndividualParticles();
     const auto& pairList = integrator.getAllThermostatedPairs();

platforms/reference/include/ReferenceKernels.h

@@ -1204,10 +1204,8 @@ public:
      * 
      * @param context    the context in which to execute this kernel
      * @param integrator the VerletIntegrator this kernel is being used for
-     * @param forcesAreValid a reference to the parent integrator's boolean for keeping
-     *                       track of the validity of the current forces.
      */
-    void execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid);
+    void execute(ContextImpl& context, const NoseHooverIntegrator& integrator);
     /**
      * Compute the kinetic energy.
      * 

platforms/reference/include/ReferenceNoseHooverDynamics.h

 
-/* Portions copyright (c) 2006-2020 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2024 Stanford University and Simbios.
  * Contributors: Andy Simmonett, Peter Eastman, Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -73,7 +73,6 @@ class ReferenceNoseHooverDynamics : public ReferenceDynamics {
          @param forces              forces
          @param masses              atom masses
          @param tolerance           the constraint tolerance
-         @param forcesAreValid      whether the forces are valid (duh!)
          @param allAtoms            a list of all atoms not involved in a Drude-like pair
          @param allPairs            a list of all Drude-like pairs, and their KT values, in the system
          @param maxPairDistance     the maximum separation allowed for a Drude-like pair
@@ -81,7 +80,7 @@ class ReferenceNoseHooverDynamics : public ReferenceDynamics {
          --------------------------------------------------------------------------------------- */
      
       void step1(OpenMM::ContextImpl &context, const OpenMM::System& system, std::vector<OpenMM::Vec3>& atomCoordinates,
-                 std::vector<OpenMM::Vec3>& velocities, std::vector<OpenMM::Vec3>& forces, std::vector<double>& masses, double tolerance, bool &forcesAreValid,
+                 std::vector<OpenMM::Vec3>& velocities, std::vector<OpenMM::Vec3>& forces, std::vector<double>& masses, double tolerance,
                  const std::vector<int> & allAtoms, const std::vector<std::tuple<int, int, double>> & allPairs, double maxPairDistance);
       /**---------------------------------------------------------------------------------------
       
@@ -93,14 +92,13 @@ class ReferenceNoseHooverDynamics : public ReferenceDynamics {
          @param forces              forces
          @param masses              atom masses
          @param tolerance           the constraint tolerance
-         @param forcesAreValid      whether the forces are valid (duh!)
          @param allAtoms            a list of all atoms not involved in a Drude-like pair
          @param allPairs            a list of all Drude-like pairs, and their KT values, in the system
          @param maxPairDistance     the maximum separation allowed for a Drude-like pair
       
          --------------------------------------------------------------------------------------- */
       void step2(OpenMM::ContextImpl &context, const OpenMM::System& system, std::vector<OpenMM::Vec3>& atomCoordinates,
-                 std::vector<OpenMM::Vec3>& velocities, std::vector<OpenMM::Vec3>& forces, std::vector<double>& masses, double tolerance, bool &forcesAreValid,
+                 std::vector<OpenMM::Vec3>& velocities, std::vector<OpenMM::Vec3>& forces, std::vector<double>& masses, double tolerance,
                  const std::vector<int> & allAtoms, const std::vector<std::tuple<int, int, double>> & allPairs, double maxPairDistance);
       
 };

platforms/reference/src/ReferenceKernels.cpp

@@ -2348,7 +2348,7 @@ void ReferenceIntegrateNoseHooverStepKernel::initialize(const System& system, co
     this->chainPropagator = new ReferenceNoseHooverChain();
 }
 
-void ReferenceIntegrateNoseHooverStepKernel::execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid) {
+void ReferenceIntegrateNoseHooverStepKernel::execute(ContextImpl& context, const NoseHooverIntegrator& integrator) {
     double stepSize = integrator.getStepSize();
     vector<Vec3>& posData = extractPositions(context);
     vector<Vec3>& velData = extractVelocities(context);
@@ -2362,7 +2362,7 @@ void ReferenceIntegrateNoseHooverStepKernel::execute(ContextImpl& context, const
         dynamics->setVirtualSites(extractVirtualSites(context));
         prevStepSize = stepSize;
     }
-    dynamics->step1(context, context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance(), forcesAreValid,
+    dynamics->step1(context, context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance(),
                      integrator.getAllThermostatedIndividualParticles(), integrator.getAllThermostatedPairs(), integrator.getMaximumPairDistance());
     int numChains = integrator.getNumThermostats();
     for(int chain = 0; chain < numChains; ++chain) {
@@ -2371,7 +2371,7 @@ void ReferenceIntegrateNoseHooverStepKernel::execute(ContextImpl& context, const
         std::pair<double, double> scaleFactors = propagateChain(context, thermostatChain, KEs, stepSize);
         scaleVelocities(context, thermostatChain, scaleFactors);
     }
-    dynamics->step2(context, context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance(), forcesAreValid,
+    dynamics->step2(context, context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance(),
                      integrator.getAllThermostatedIndividualParticles(), integrator.getAllThermostatedPairs(), integrator.getMaximumPairDistance());
     data.time += stepSize;
     data.stepCount++;

platforms/reference/src/SimTKReference/ReferenceNoseHooverDynamics.cpp

 
-/* Portions copyright (c) 2006-2020 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2024 Stanford University and Simbios.
  * Contributors: Andy Simmonett, Peter Eastman, Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -51,13 +51,11 @@ ReferenceNoseHooverDynamics::~ReferenceNoseHooverDynamics() {
 
 void ReferenceNoseHooverDynamics::step1(OpenMM::ContextImpl &context, const OpenMM::System& system, vector<Vec3>& atomCoordinates,
                                           vector<Vec3>& velocities,
-                                          vector<Vec3>& forces, vector<double>& masses, double tolerance, bool &forcesAreValid,
+                                          vector<Vec3>& forces, vector<double>& masses, double tolerance,
                                           const std::vector<int> & atomList, const std::vector<std::tuple<int, int, double>> &pairList,
                                           double maxPairDistance) {
 
     // first-time-through initialization
-    if (!forcesAreValid) context.calcForcesAndEnergy(true, false, context.getIntegrator().getIntegrationForceGroups());
-
     if (getTimeStep() == 0) {
        // invert masses
        for (int ii = 0; ii < numberOfAtoms; ii++) {
@@ -115,7 +113,7 @@ void ReferenceNoseHooverDynamics::step1(OpenMM::ContextImpl &context, const Open
 
 void ReferenceNoseHooverDynamics::step2(OpenMM::ContextImpl &context, const OpenMM::System& system, vector<Vec3>& atomCoordinates,
                                           vector<Vec3>& velocities,
-                                          vector<Vec3>& forces, vector<double>& masses, double tolerance, bool &forcesAreValid,
+                                          vector<Vec3>& forces, vector<double>& masses, double tolerance,
                                           const std::vector<int> & atomList, const std::vector<std::tuple<int, int, double>> &pairList,
                                           double maxPairDistance) {
     const double halfdt = 0.5*getDeltaT();