Reduced unnecessary force evaluations with multiple time step integrators

openmmapi/include/openmm/internal/ContextImpl.h

@@ -186,6 +186,10 @@ public:
      * @return the potential energy of the system, or 0 if includeEnergy is false
      */
     double calcForcesAndEnergy(bool includeForces, bool includeEnergy, int groups=0xFFFFFFFF);
+    /**
+     * Get the set of force group flags that were passed to the most recent call to calcForcesAndEnergy().
+     */
+    int getLastForceGroups() const;
     /**
      * Calculate the kinetic energy of the system (in kJ/mol).
      */
@@ -226,6 +230,7 @@ private:
     std::map<std::string, double> parameters;
     mutable std::vector<std::vector<int> > molecules;
     bool hasInitializedForces;
+    int lastForceGroups;
     Platform* platform;
     Kernel initializeForcesKernel, kineticEnergyKernel, updateStateDataKernel, applyConstraintsKernel, virtualSitesKernel;
     void* platformData;

openmmapi/src/ContextImpl.cpp

@@ -49,7 +49,7 @@ using std::vector;
 using std::string;
 
 ContextImpl::ContextImpl(Context& owner, System& system, Integrator& integrator, Platform* platform, const map<string, string>& properties) :
-         owner(owner), system(system), integrator(integrator), hasInitializedForces(false), platform(platform), platformData(NULL) {
+         owner(owner), system(system), integrator(integrator), hasInitializedForces(false), lastForceGroups(-1), platform(platform), platformData(NULL) {
     if (system.getNumParticles() == 0)
         throw OpenMMException("Cannot create a Context for a System with no particles");
     
@@ -192,6 +192,7 @@ void ContextImpl::computeVirtualSites() {
 }
 
 double ContextImpl::calcForcesAndEnergy(bool includeForces, bool includeEnergy, int groups) {
+    lastForceGroups = groups;
     CalcForcesAndEnergyKernel& kernel = dynamic_cast<CalcForcesAndEnergyKernel&>(initializeForcesKernel.getImpl());
     double energy = 0.0;
     kernel.beginComputation(*this, includeForces, includeEnergy, groups);
@@ -201,6 +202,10 @@ double ContextImpl::calcForcesAndEnergy(bool includeForces, bool includeEnergy,
     return energy;
 }
 
+int ContextImpl::getLastForceGroups() const {
+    return lastForceGroups;
+}
+
 double ContextImpl::calcKineticEnergy() {
     return dynamic_cast<CalcKineticEnergyKernel&>(kineticEnergyKernel.getImpl()).execute(*this);
 }

platforms/reference/src/SimTKReference/ReferenceCustomDynamics.cpp

@@ -91,7 +91,6 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
                                      vector<RealVec>& velocities, vector<RealVec>& forces, vector<RealOpenMM>& masses,
                                      map<string, RealOpenMM>& globals, vector<vector<RealVec> >& perDof, bool& forcesAreValid){
     int numSteps = stepType.size();
-    int currentForceGroup = -1;
     globals.insert(context.getParameters().begin(), context.getParameters().end());
     if (invalidatesForces.size() == 0) {
         // The first time this is called, work out when to recompute forces and energy.  First build a
@@ -167,7 +166,7 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
     // Loop over steps and execute them.
     
     for (int i = 0; i < numSteps; i++) {
-        if ((needsForces[i] || needsEnergy[i]) && (!forcesAreValid || currentForceGroup != forceGroup[i])) {
+        if ((needsForces[i] || needsEnergy[i]) && (!forcesAreValid || context.getLastForceGroups() != forceGroup[i])) {
             // Recompute forces and or energy.  Figure out what is actually needed
             // between now and the next time they get invalidated again.
             
@@ -189,7 +188,6 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
             if (computeEnergy)
                 energy = e;
             forcesAreValid = true;
-            currentForceGroup = forceGroup[i];
         }
         globals["energy"] = energy;
         
@@ -248,8 +246,6 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
     ReferenceVirtualSites::computePositions(context.getSystem(), atomCoordinates);
     incrementTimeStep();
     recordChangedParameters(context, globals);
-    if (currentForceGroup != -1)
-        forcesAreValid = false;
 }
 
 void ReferenceCustomDynamics::computePerDof(int numberOfAtoms, vector<RealVec>& results, const vector<RealVec>& atomCoordinates,