Created architecture for multiple time step integrators and completed reference implementation

openmmapi/src/PeriodicTorsionForceImpl.cpp

@@ -49,7 +49,8 @@ void PeriodicTorsionForceImpl::initialize(ContextImpl& context) {
     dynamic_cast<CalcPeriodicTorsionForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
 }
 
-double PeriodicTorsionForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy) {
+double PeriodicTorsionForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
+    if ((groups&(1<<owner.getForceGroup())) != 0)
         return dynamic_cast<CalcPeriodicTorsionForceKernel&>(kernel.getImpl()).execute(context, includeForces, includeEnergy);
 }
 

openmmapi/src/RBTorsionForceImpl.cpp

@@ -49,7 +49,8 @@ void RBTorsionForceImpl::initialize(ContextImpl& context) {
     dynamic_cast<CalcRBTorsionForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
 }
 
-double RBTorsionForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy) {
+double RBTorsionForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
+    if ((groups&(1<<owner.getForceGroup())) != 0)
         return dynamic_cast<CalcRBTorsionForceKernel&>(kernel.getImpl()).execute(context, includeForces, includeEnergy);
 }
 

platforms/cuda/src/CudaKernels.cpp

@@ -46,7 +46,7 @@ using namespace std;
 void CudaCalcForcesAndEnergyKernel::initialize(const System& system) {
 }
 
-void CudaCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, bool includeForces, bool includeEnergy) {
+void CudaCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
     _gpuContext* gpu = data.gpu;
     if (data.nonbondedMethod != NO_CUTOFF && data.computeForceCount%100 == 0)
         gpuReorderAtoms(gpu);
@@ -59,7 +59,7 @@ void CudaCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, bool
         kClearEnergy(gpu);
 }
 
-double CudaCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, bool includeForces, bool includeEnergy) {
+double CudaCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
     _gpuContext* gpu = data.gpu;
     if (gpu->bIncludeGBSA || gpu->bIncludeGBVI) {
         gpu->bRecalculateBornRadii = true;
@@ -874,7 +874,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
     data.gpu->forces.push_back(new ForceInfo(force));
 }
 
-double CudaCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+double CudaCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal) {
     return 0.0;
 }
 

platforms/cuda/src/CudaKernels.h

@@ -65,8 +65,9 @@ public:
      * @param context       the context in which to execute this kernel
      * @param includeForce  true if forces should be computed
      * @param includeEnergy true if potential energy should be computed
+     * @param groups        a set of bit flags for which force groups to include
      */
-    void beginComputation(ContextImpl& context, bool includeForce, bool includeEnergy);
+    void beginComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups);
     /**
      * This is called at the end of each force/energy computation, after calcForcesAndEnergy() has been called on
      * every ForceImpl.
@@ -74,11 +75,12 @@ public:
      * @param context       the context in which to execute this kernel
      * @param includeForce  true if forces should be computed
      * @param includeEnergy true if potential energy should be computed
+     * @param groups        a set of bit flags for which force groups to include
      * @return the potential energy of the system.  This value is added to all values returned by ForceImpls'
      * calcForcesAndEnergy() methods.  That is, each force kernel may <i>either</i> return its contribution to the
      * energy directly, <i>or</i> add it to an internal buffer so that it will be included here.
      */
-    double finishComputation(ContextImpl& context, bool includeForce, bool includeEnergy);
+    double finishComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups);
 private:
     CudaPlatform::PlatformData& data;
 };
@@ -491,9 +493,10 @@ public:
      * @param context        the context in which to execute this kernel
      * @param includeForces  true if forces should be calculated
      * @param includeEnergy  true if the energy should be calculated
+     * @param includeReciprocal  true if reciprocal space interactions should be included
      * @return the potential energy due to the force
      */
-    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal);
 private:
     class ForceInfo;
     CudaPlatform::PlatformData& data;

platforms/opencl/src/OpenCLKernels.cpp

@@ -93,7 +93,7 @@ static pair<ExpressionTreeNode, string> makeVariable(const string& name, const s
 void OpenCLCalcForcesAndEnergyKernel::initialize(const System& system) {
 }
 
-void OpenCLCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, bool includeForces, bool includeEnergy) {
+void OpenCLCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
     cl.setAtomsWereReordered(false);
     if (cl.getNonbondedUtilities().getUseCutoff() && cl.getComputeForceCount()%100 == 0) {
         cl.reorderAtoms();
@@ -104,7 +104,7 @@ void OpenCLCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, boo
     cl.getNonbondedUtilities().prepareInteractions();
 }
 
-double OpenCLCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, bool includeForces, bool includeEnergy) {
+double OpenCLCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
     cl.getBondedUtilities().computeInteractions();
     cl.getNonbondedUtilities().computeInteractions();
     cl.reduceForces();
@@ -1197,7 +1197,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
     cl.addForce(new OpenCLNonbondedForceInfo(cl.getNonbondedUtilities().getNumForceBuffers(), force));
 }
 
-double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal) {
     bool deviceIsCpu = (cl.getDevice().getInfo<CL_DEVICE_TYPE>() == CL_DEVICE_TYPE_CPU);
     if (!hasInitializedKernel) {
         hasInitializedKernel = true;

platforms/opencl/src/OpenCLKernels.h

@@ -60,8 +60,9 @@ public:
      * @param context       the context in which to execute this kernel
      * @param includeForce  true if forces should be computed
      * @param includeEnergy true if potential energy should be computed
+     * @param groups        a set of bit flags for which force groups to include
      */
-    void beginComputation(ContextImpl& context, bool includeForce, bool includeEnergy);
+    void beginComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups);
     /**
      * This is called at the end of each force/energy computation, after calcForcesAndEnergy() has been called on
      * every ForceImpl.
@@ -69,11 +70,12 @@ public:
      * @param context       the context in which to execute this kernel
      * @param includeForce  true if forces should be computed
      * @param includeEnergy true if potential energy should be computed
+     * @param groups        a set of bit flags for which force groups to include
      * @return the potential energy of the system.  This value is added to all values returned by ForceImpls'
      * calcForcesAndEnergy() methods.  That is, each force kernel may <i>either</i> return its contribution to the
      * energy directly, <i>or</i> add it to an internal buffer so that it will be included here.
      */
-    double finishComputation(ContextImpl& context, bool includeForce, bool includeEnergy);
+    double finishComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups);
 private:
    OpenCLContext& cl;
 };
@@ -503,9 +505,11 @@ public:
      * @param context        the context in which to execute this kernel
      * @param includeForces  true if forces should be calculated
      * @param includeEnergy  true if the energy should be calculated
+     * @param includeDirect  true if direct space interactions should be included
+     * @param includeReciprocal  true if reciprocal space interactions should be included
      * @return the potential energy due to the force
      */
-    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal);
 private:
     struct SortTrait {
         typedef mm_int2 DataType;

platforms/opencl/src/OpenCLParallelKernels.cpp

@@ -54,35 +54,36 @@ using namespace std;
 class OpenCLParallelCalcForcesAndEnergyKernel::BeginComputationTask : public OpenCLContext::WorkTask {
 public:
     BeginComputationTask(ContextImpl& context, OpenCLContext& cl, OpenCLCalcForcesAndEnergyKernel& kernel,
-            bool includeForce, bool includeEnergy, mm_float4* pinnedMemory) : context(context), cl(cl), kernel(kernel),
-            includeForce(includeForce), includeEnergy(includeEnergy), pinnedMemory(pinnedMemory) {
+            bool includeForce, bool includeEnergy, int groups, mm_float4* pinnedMemory) : context(context), cl(cl), kernel(kernel),
+            includeForce(includeForce), includeEnergy(includeEnergy), groups(groups), pinnedMemory(pinnedMemory) {
     }
     void execute() {
         // Copy coordinates over to this device and execute the kernel.
 
         if (cl.getContextIndex() > 0)
             cl.getQueue().enqueueWriteBuffer(cl.getPosq().getDeviceBuffer(), CL_FALSE, 0, cl.getPaddedNumAtoms()*sizeof(mm_float4), pinnedMemory);
-        kernel.beginComputation(context, includeForce, includeEnergy);
+        kernel.beginComputation(context, includeForce, includeEnergy, groups);
     }
 private:
     ContextImpl& context;
     OpenCLContext& cl;
     OpenCLCalcForcesAndEnergyKernel& kernel;
     bool includeForce, includeEnergy;
+    int groups;
     mm_float4* pinnedMemory;
 };
 
 class OpenCLParallelCalcForcesAndEnergyKernel::FinishComputationTask : public OpenCLContext::WorkTask {
 public:
     FinishComputationTask(ContextImpl& context, OpenCLContext& cl, OpenCLCalcForcesAndEnergyKernel& kernel,
-            bool includeForce, bool includeEnergy, double& energy, long long& completionTime, mm_float4* pinnedMemory) :
-            context(context), cl(cl), kernel(kernel), includeForce(includeForce), includeEnergy(includeEnergy), energy(energy),
+            bool includeForce, bool includeEnergy, int groups, double& energy, long long& completionTime, mm_float4* pinnedMemory) :
+            context(context), cl(cl), kernel(kernel), includeForce(includeForce), includeEnergy(includeEnergy), groups(groups), energy(energy),
             completionTime(completionTime), pinnedMemory(pinnedMemory) {
     }
     void execute() {
         // Execute the kernel, then download forces.
         
-        energy += kernel.finishComputation(context, includeForce, includeEnergy);
+        energy += kernel.finishComputation(context, includeForce, includeEnergy, groups);
         if (includeForce) {
             if (cl.getContextIndex() > 0) {
                 int numAtoms = cl.getPaddedNumAtoms();
@@ -99,6 +100,7 @@ private:
     OpenCLContext& cl;
     OpenCLCalcForcesAndEnergyKernel& kernel;
     bool includeForce, includeEnergy;
+    int groups;
     double& energy;
     long long& completionTime;
     mm_float4* pinnedMemory;
@@ -125,7 +127,7 @@ void OpenCLParallelCalcForcesAndEnergyKernel::initialize(const System& system) {
         getKernel(i).initialize(system);
 }
 
-void OpenCLParallelCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, bool includeForce, bool includeEnergy) {
+void OpenCLParallelCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups) {
     OpenCLContext& cl0 = *data.contexts[0];
     if (contextForces == NULL) {
         contextForces = new OpenCLArray<mm_float4>(cl0, &cl0.getForceBuffers().getDeviceBuffer(),
@@ -144,15 +146,15 @@ void OpenCLParallelCalcForcesAndEnergyKernel::beginComputation(ContextImpl& cont
         data.contextEnergy[i] = 0.0;
         OpenCLContext& cl = *data.contexts[i];
         OpenCLContext::WorkThread& thread = cl.getWorkThread();
-        thread.addTask(new BeginComputationTask(context, cl, getKernel(i), includeForce, includeEnergy, pinnedPositionMemory));
+        thread.addTask(new BeginComputationTask(context, cl, getKernel(i), includeForce, includeEnergy, groups, pinnedPositionMemory));
     }
 }
 
-double OpenCLParallelCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, bool includeForce, bool includeEnergy) {
+double OpenCLParallelCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups) {
     for (int i = 0; i < (int) data.contexts.size(); i++) {
         OpenCLContext& cl = *data.contexts[i];
         OpenCLContext::WorkThread& thread = cl.getWorkThread();
-        thread.addTask(new FinishComputationTask(context, cl, getKernel(i), includeForce, includeEnergy, data.contextEnergy[i], completionTimes[i], pinnedForceMemory));
+        thread.addTask(new FinishComputationTask(context, cl, getKernel(i), includeForce, includeEnergy, groups, data.contextEnergy[i], completionTimes[i], pinnedForceMemory));
     }
     data.syncContexts();
     double energy = 0.0;
@@ -487,16 +489,16 @@ double OpenCLParallelCalcCustomTorsionForceKernel::execute(ContextImpl& context,
 class OpenCLParallelCalcNonbondedForceKernel::Task : public OpenCLContext::WorkTask {
 public:
     Task(ContextImpl& context, OpenCLCalcNonbondedForceKernel& kernel, bool includeForce,
-            bool includeEnergy, double& energy) : context(context), kernel(kernel),
-            includeForce(includeForce), includeEnergy(includeEnergy), energy(energy) {
+            bool includeEnergy, bool includeDirect, bool includeReciprocal, double& energy) : context(context), kernel(kernel),
+            includeForce(includeForce), includeEnergy(includeEnergy), includeDirect(includeDirect), includeReciprocal(includeReciprocal), energy(energy) {
     }
     void execute() {
-        energy += kernel.execute(context, includeForce, includeEnergy);
+        energy += kernel.execute(context, includeForce, includeEnergy, includeDirect, includeReciprocal);
     }
 private:
     ContextImpl& context;
     OpenCLCalcNonbondedForceKernel& kernel;
-    bool includeForce, includeEnergy;
+    bool includeForce, includeEnergy, includeDirect, includeReciprocal;
     double& energy;
 };
 
@@ -511,11 +513,11 @@ void OpenCLParallelCalcNonbondedForceKernel::initialize(const System& system, co
         getKernel(i).initialize(system, force);
 }
 
-double OpenCLParallelCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+double OpenCLParallelCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal) {
     for (int i = 0; i < (int) data.contexts.size(); i++) {
         OpenCLContext& cl = *data.contexts[i];
         OpenCLContext::WorkThread& thread = cl.getWorkThread();
-        thread.addTask(new Task(context, getKernel(i), includeForces, includeEnergy, data.contextEnergy[i]));
+        thread.addTask(new Task(context, getKernel(i), includeForces, includeEnergy, includeDirect, includeReciprocal, data.contextEnergy[i]));
     }
     return 0.0;
 }

platforms/opencl/src/OpenCLParallelKernels.h

@@ -58,8 +58,9 @@ public:
      * @param context       the context in which to execute this kernel
      * @param includeForce  true if forces should be computed
      * @param includeEnergy true if potential energy should be computed
+     * @param groups        a set of bit flags for which force groups to include
      */
-    void beginComputation(ContextImpl& context, bool includeForce, bool includeEnergy);
+    void beginComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups);
     /**
      * This is called at the end of each force/energy computation, after calcForcesAndEnergy() has been called on
      * every ForceImpl.
@@ -67,11 +68,12 @@ public:
      * @param context       the context in which to execute this kernel
      * @param includeForce  true if forces should be computed
      * @param includeEnergy true if potential energy should be computed
+     * @param groups        a set of bit flags for which force groups to include
      * @return the potential energy of the system.  This value is added to all values returned by ForceImpls'
      * calcForcesAndEnergy() methods.  That is, each force kernel may <i>either</i> return its contribution to the
      * energy directly, <i>or</i> add it to an internal buffer so that it will be included here.
      */
-    double finishComputation(ContextImpl& context, bool includeForce, bool includeEnergy);
+    double finishComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups);
 private:
     class BeginComputationTask;
     class FinishComputationTask;
@@ -356,9 +358,11 @@ public:
      * @param context        the context in which to execute this kernel
      * @param includeForces  true if forces should be calculated
      * @param includeEnergy  true if the energy should be calculated
+     * @param includeReciprocal  true if reciprocal space interactions should be included
+     * @param includeReciprocal  true if reciprocal space interactions should be included
      * @return the potential energy due to the force
      */
-    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal);
 private:
     class Task;
     OpenCLPlatform::PlatformData& data;

platforms/reference/src/ReferenceKernels.cpp

@@ -145,7 +145,7 @@ static void findAnglesForCCMA(const System& system, vector<ReferenceCCMAAlgorith
 void ReferenceCalcForcesAndEnergyKernel::initialize(const System& system) {
 }
 
-void ReferenceCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, bool includeForces, bool includeEnergy) {
+void ReferenceCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
     vector<RealVec>& forceData = extractForces(context);
     if (includeForces) {
         int numParticles = context.getSystem().getNumParticles();
@@ -159,7 +159,7 @@ void ReferenceCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context,
         savedForces = forceData;
 }
 
-double ReferenceCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, bool includeForces, bool includeEnergy) {
+double ReferenceCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
     if (!includeForces)
         extractForces(context) = savedForces; // Restore the forces so computing the energy doesn't overwrite the forces with incorrect values.
     else
@@ -676,7 +676,7 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
         dispersionCoefficient = 0.0;
 }
 
-double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal) {
     vector<RealVec>& posData = extractPositions(context);
     vector<RealVec>& forceData = extractForces(context);
     RealOpenMM energy = 0;
@@ -694,10 +694,12 @@ double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool inc
         clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
     if (pme)
         clj.setUsePME(ewaldAlpha, gridSize);
-    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, includeEnergy ? &energy : NULL);
+    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, includeEnergy ? &energy : NULL, includeDirect, includeReciprocal);
+    if (includeDirect) {
         ReferenceBondForce refBondForce;
         ReferenceLJCoulomb14 nonbonded14;
         refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, includeEnergy ? &energy : NULL, nonbonded14);
+    }
     if (periodic || ewald || pme) {
         RealVec& boxSize = extractBoxSize(context);
         energy += dispersionCoefficient/(boxSize[0]*boxSize[1]*boxSize[2]);

platforms/reference/src/ReferenceKernels.h

@@ -75,8 +75,9 @@ public:
      * @param context       the context in which to execute this kernel
      * @param includeForce  true if forces should be computed
      * @param includeEnergy true if potential energy should be computed
+     * @param groups        a set of bit flags for which force groups to include
      */
-    void beginComputation(ContextImpl& context, bool includeForce, bool includeEnergy);
+    void beginComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups);
     /**
      * This is called at the end of each force/energy computation, after calcForcesAndEnergy() has been called on
      * every ForceImpl.
@@ -84,11 +85,12 @@ public:
      * @param context       the context in which to execute this kernel
      * @param includeForce  true if forces should be computed
      * @param includeEnergy true if potential energy should be computed
+     * @param groups        a set of bit flags for which force groups to include
      * @return the potential energy of the system.  This value is added to all values returned by ForceImpls'
      * calcForcesAndEnergy() methods.  That is, each force kernel may <i>either</i> return its contribution to the
      * energy directly, <i>or</i> add it to an internal buffer so that it will be included here.
      */
-    double finishComputation(ContextImpl& context, bool includeForce, bool includeEnergy);
+    double finishComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups);
 private:
     std::vector<RealVec> savedForces;
 };
@@ -488,9 +490,10 @@ public:
      * @param context        the context in which to execute this kernel
      * @param includeForces  true if forces should be calculated
      * @param includeEnergy  true if the energy should be calculated
+     * @param includeReciprocal  true if reciprocal space interactions should be included
      * @return the potential energy due to the force
      */
-    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal);
 private:
     int numParticles, num14;
     int **exclusionArray, **bonded14IndexArray;

platforms/reference/src/SimTKReference/ReferenceCustomDynamics.cpp

@@ -34,6 +34,7 @@
 #include "lepton/ParsedExpression.h"
 #include "lepton/Parser.h"
 #include <set>
+#include <sstream>
 
 using namespace std;
 using namespace OpenMM;
@@ -90,6 +91,7 @@ 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
@@ -98,6 +100,8 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
         invalidatesForces.resize(numSteps, false);
         needsForces.resize(numSteps, false);
         needsEnergy.resize(numSteps, false);
+        forceGroup.resize(numSteps, -2);
+        forceName.resize(numSteps, "f");
         set<string> affectsForce;
         affectsForce.insert("x");
         for (vector<ForceImpl*>::const_iterator iter = context.getForceImpls().begin(); iter != context.getForceImpls().end(); ++iter) {
@@ -110,15 +114,40 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
         
         // Make a list of which steps require valid forces or energy to be known.
         
+        vector<string> forceGroupName;
+        for (int i = 0; i < 32; i++) {
+            stringstream str;
+            str << "f" << i;
+            forceGroupName.push_back(str.str());
+        }
         for (int i = 0; i < numSteps; i++) {
             if (stepType[i] == CustomIntegrator::ComputeGlobal || stepType[i] == CustomIntegrator::ComputePerDof || stepType[i] == CustomIntegrator::ComputeSum) {
                 for (int j = 0; j < stepExpression[i].getNumOperations(); j++) {
                     const Lepton::Operation& op = stepExpression[i].getOperation(j);
                     if (op.getId() == Lepton::Operation::VARIABLE) {
-                        if (op.getName() == "f")
-                            needsForces[i] = true;
-                        else if (op.getName() == "energy")
+                        if (op.getName() == "energy") {
+                            if (forceGroup[i] != -2)
+                                throw OpenMMException("A single computation step cannot depend on multiple force groups");
                             needsEnergy[i] = true;
+                            forceGroup[i] = -1;
+                        }
+                        else if (op.getName() == "f") {
+                            if (forceGroup[i] != -2)
+                                throw OpenMMException("A single computation step cannot depend on multiple force groups");
+                            needsForces[i] = true;
+                            forceGroup[i] = -1;
+                        }
+                        else if (op.getName()[0] == 'f') {
+                            for (int k = 0; k < (int) forceGroupName.size(); k++)
+                                if (op.getName() == forceGroupName[k]) {
+                                    if (forceGroup[i] != -2)
+                                        throw OpenMMException("A single computation step cannot depend on multiple force groups");
+                                    needsForces[i] = true;
+                                    forceGroup[i] = 1<<k;
+                                    forceName[i] = forceGroupName[k];
+                                    break;
+                                }
+                        }
                     }
                 }
             }
@@ -138,7 +167,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) {
+        if ((needsForces[i] || needsEnergy[i]) && (!forcesAreValid || currentForceGroup != forceGroup[i])) {
             // Recompute forces and or energy.  Figure out what is actually needed
             // between now and the next time they get invalidated again.
             
@@ -156,10 +185,11 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
                     break;
             }
             recordChangedParameters(context, globals);
-            RealOpenMM e = context.calcForcesAndEnergy(computeForce, computeEnergy);
+            RealOpenMM e = context.calcForcesAndEnergy(computeForce, computeEnergy, forceGroup[i]);
             if (computeEnergy)
                 energy = e;
             forcesAreValid = true;
+            currentForceGroup = forceGroup[i];
         }
         globals["energy"] = energy;
         
@@ -186,11 +216,11 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
                 }
                 if (results == NULL)
                     throw OpenMMException("Illegal per-DOF output variable: "+stepVariable[i]);
-                computePerDof(numberOfAtoms, *results, atomCoordinates, velocities, forces, masses, globals, perDof, stepExpression[i]);
+                computePerDof(numberOfAtoms, *results, atomCoordinates, velocities, forces, masses, globals, perDof, stepExpression[i], forceName[i]);
                 break;
             }
             case CustomIntegrator::ComputeSum: {
-                computePerDof(numberOfAtoms, sumBuffer, atomCoordinates, velocities, forces, masses, globals, perDof, stepExpression[i]);
+                computePerDof(numberOfAtoms, sumBuffer, atomCoordinates, velocities, forces, masses, globals, perDof, stepExpression[i], forceName[i]);
                 RealOpenMM sum = 0.0;
                 for (int j = 0; j < numberOfAtoms; j++)
                     if (masses[j] != 0.0)
@@ -218,11 +248,14 @@ 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,
               const vector<RealVec>& velocities, const vector<RealVec>& forces, const vector<RealOpenMM>& masses,
-              const map<string, RealOpenMM>& globals, const vector<vector<RealVec> >& perDof, const Lepton::ExpressionProgram& expression) {
+              const map<string, RealOpenMM>& globals, const vector<vector<RealVec> >& perDof,
+              const Lepton::ExpressionProgram& expression, const std::string& forceName) {
     // Loop over all degrees of freedom.
     
     map<string, RealOpenMM> variables = globals;
@@ -234,7 +267,7 @@ void ReferenceCustomDynamics::computePerDof(int numberOfAtoms, vector<RealVec>&
 
                 variables["x"] = atomCoordinates[i][j];
                 variables["v"] = velocities[i][j];
-                variables["f"] = forces[i][j];
+                variables[forceName] = forces[i][j];
                 variables["uniform"] = SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber();
                 variables["gaussian"] = SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
                 for (int k = 0; k < (int) perDof.size(); k++)

platforms/reference/src/SimTKReference/ReferenceCustomDynamics.h

@@ -43,14 +43,16 @@ private:
     std::vector<RealOpenMM> inverseMasses;
     std::vector<OpenMM::RealVec> sumBuffer;
     std::vector<OpenMM::CustomIntegrator::ComputationType> stepType;
-    std::vector<std::string> stepVariable;
+    std::vector<std::string> stepVariable, forceName;
     std::vector<Lepton::ExpressionProgram> stepExpression;
     std::vector<bool> invalidatesForces, needsForces, needsEnergy;
+    std::vector<int> forceGroup;
     RealOpenMM energy;
     
     void computePerDof(int numberOfAtoms, std::vector<OpenMM::RealVec>& results, const std::vector<OpenMM::RealVec>& atomCoordinates,
                   const std::vector<OpenMM::RealVec>& velocities, const std::vector<OpenMM::RealVec>& forces, const std::vector<RealOpenMM>& masses,
-                  const std::map<std::string, RealOpenMM>& globals, const std::vector<std::vector<OpenMM::RealVec> >& perDof, const Lepton::ExpressionProgram& expression);
+                  const std::map<std::string, RealOpenMM>& globals, const std::vector<std::vector<OpenMM::RealVec> >& perDof,
+                  const Lepton::ExpressionProgram& expression, const std::string& forceName);
     
     void recordChangedParameters(OpenMM::ContextImpl& context, std::map<std::string, RealOpenMM>& globals);
       

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp

@@ -164,13 +164,15 @@ ReferenceLJCoulombIxn::~ReferenceLJCoulombIxn( ){
    @param forces           force array (forces added)
    @param energyByAtom     atom energy
    @param totalEnergy      total energy
+   @param includeDirect      true if direct space interactions should be included
+   @param includeReciprocal  true if reciprocal space interactions should be included
 
    --------------------------------------------------------------------------------------- */
 
-void ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, vector<RealVec>& atomCoordinates,
+void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>& atomCoordinates,
                                              RealOpenMM** atomParameters, int** exclusions,
                                              RealOpenMM* fixedParameters, vector<RealVec>& forces,
-                                             RealOpenMM* energyByAtom, RealOpenMM* totalEnergy) const {
+                                             RealOpenMM* energyByAtom, RealOpenMM* totalEnergy, bool includeDirect, bool includeReciprocal) const {
     typedef std::complex<RealOpenMM> d_complex;
 
     static const RealOpenMM epsilon     =  1.0;
@@ -212,7 +214,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, vector<RealVec
 // **************************************************************************************
     // PME
 
-  if (pme) {
+  if (pme && includeReciprocal) {
 	pme_t          pmedata; /* abstract handle for PME data */
 	RealOpenMM virial[3][3];
 
@@ -232,7 +234,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, vector<RealVec
 
     // Ewald method
 
-  else if (ewald) {
+  else if (ewald && includeReciprocal) {
 
     // setup reciprocal box
 
@@ -336,17 +338,12 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, vector<RealVec
     }
   }
 
-  else {
-      std::stringstream message;
-      message << " Wrong method for Ewald summation, Aborting" << std::endl;
-      SimTKOpenMMLog::printError( message );
-  }
-
-
 // **************************************************************************************
 // SHORT-RANGE ENERGY AND FORCES
 // **************************************************************************************
 
+    if (!includeDirect)
+        return;
     RealOpenMM totalVdwEnergy            = 0.0f;
     RealOpenMM totalRealSpaceEwaldEnergy = 0.0f;
 
@@ -454,16 +451,23 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, vector<RealVec
    @param forces           force array (forces added)
    @param energyByAtom     atom energy
    @param totalEnergy      total energy
+   @param includeDirect      true if direct space interactions should be included
+   @param includeReciprocal  true if reciprocal space interactions should be included
 
    --------------------------------------------------------------------------------------- */
 
-void ReferenceLJCoulombIxn::calculatePairIxn( int numberOfAtoms, vector<RealVec>& atomCoordinates,
+void ReferenceLJCoulombIxn::calculatePairIxn(int numberOfAtoms, vector<RealVec>& atomCoordinates,
                                              RealOpenMM** atomParameters, int** exclusions,
                                              RealOpenMM* fixedParameters, vector<RealVec>& forces,
-                                             RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
+                                             RealOpenMM* energyByAtom, RealOpenMM* totalEnergy, bool includeDirect, bool includeReciprocal) const {
 
-   if (ewald || pme)
-        return calculateEwaldIxn(numberOfAtoms, atomCoordinates, atomParameters, exclusions, fixedParameters, forces, energyByAtom, totalEnergy);
+   if (ewald || pme) {
+       calculateEwaldIxn(numberOfAtoms, atomCoordinates, atomParameters, exclusions, fixedParameters, forces, energyByAtom,
+               totalEnergy, includeDirect, includeReciprocal);
+       return;
+   }
+   if (!includeDirect)
+       return;
    if (cutoff) {
        for (int i = 0; i < (int) neighborList->size(); i++) {
            OpenMM::AtomPair pair = (*neighborList)[i];

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.h

@@ -30,7 +30,7 @@
 
 // ---------------------------------------------------------------------------------------
 
-class ReferenceLJCoulombIxn : public ReferencePairIxn {
+class ReferenceLJCoulombIxn {
 
    private:
        
@@ -153,13 +153,15 @@ class ReferenceLJCoulombIxn : public ReferencePairIxn {
          @param forces           force array (forces added)
          @param energyByAtom     atom energy
          @param totalEnergy      total energy
+         @param includeDirect      true if direct space interactions should be included
+         @param includeReciprocal  true if reciprocal space interactions should be included
       
          --------------------------------------------------------------------------------------- */
           
-      void calculatePairIxn( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
+      void calculatePairIxn(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
                             RealOpenMM** atomParameters, int** exclusions,
                             RealOpenMM* fixedParameters, std::vector<OpenMM::RealVec>& forces,
-                            RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
+                            RealOpenMM* energyByAtom, RealOpenMM* totalEnergy, bool includeDirect, bool includeReciprocal) const;
 
 private:
       /**---------------------------------------------------------------------------------------
@@ -177,13 +179,15 @@ private:
          @param forces           force array (forces added)
          @param energyByAtom     atom energy
          @param totalEnergy      total energy
+         @param includeDirect      true if direct space interactions should be included
+         @param includeReciprocal  true if reciprocal space interactions should be included
             
          --------------------------------------------------------------------------------------- */
           
-      void calculateEwaldIxn( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
+      void calculateEwaldIxn(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
                             RealOpenMM** atomParameters, int** exclusions,
                             RealOpenMM* fixedParameters, std::vector<OpenMM::RealVec>& forces,
-                            RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
+                            RealOpenMM* energyByAtom, RealOpenMM* totalEnergy, bool includeDirect, bool includeReciprocal) const;
 };
 
 // ---------------------------------------------------------------------------------------

platforms/reference/tests/TestReferenceCustomIntegrator.cpp

@@ -510,6 +510,107 @@ void testPerDofVariables() {
     }
 }
 
+/**
+ * Test evaluating force groups separately.
+ */
+void testForceGroups() {
+    ReferencePlatform platform;
+    System system;
+    system.addParticle(2.0);
+    system.addParticle(2.0);
+    CustomIntegrator integrator(0.01);
+    integrator.addPerDofVariable("outf", 0);
+    integrator.addPerDofVariable("outf1", 0);
+    integrator.addPerDofVariable("outf2", 0);
+    integrator.addComputePerDof("outf", "f");
+    integrator.addComputePerDof("outf1", "f1");
+    integrator.addComputePerDof("outf2", "f2");
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    bonds->addBond(0, 1, 1.5, 1.1);
+    bonds->setForceGroup(1);
+    system.addForce(bonds);
+    NonbondedForce* nb = new NonbondedForce();
+    nb->addParticle(0.2, 1, 0);
+    nb->addParticle(0.2, 1, 0);
+    nb->setForceGroup(2);
+    system.addForce(nb);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    context.setPositions(positions);
+    
+    // See if the various forces are computed correctly.
+    
+    integrator.step(1);
+    vector<Vec3> f, f1, f2;
+    integrator.getPerDofVariable(0, f);
+    integrator.getPerDofVariable(1, f1);
+    integrator.getPerDofVariable(2, f2);
+    ASSERT_EQUAL_VEC(Vec3(1.1*0.5, 0, 0), f1[0], 1e-5);
+    ASSERT_EQUAL_VEC(Vec3(-1.1*0.5, 0, 0), f1[1], 1e-5);
+    ASSERT_EQUAL_VEC(Vec3(-138.935456*0.2*0.2/4.0, 0, 0), f2[0], 1e-5);
+    ASSERT_EQUAL_VEC(Vec3(138.935456*0.2*0.2/4.0, 0, 0), f2[1], 1e-5);
+    ASSERT_EQUAL_VEC(f1[0]+f2[0], f[0], 1e-5);
+    ASSERT_EQUAL_VEC(f1[1]+f2[1], f[1], 1e-5);
+}
+
+/**
+ * Test a multiple time step r-RESPA integrator.
+ */
+void testRespa() {
+    const int numParticles = 8;
+    ReferencePlatform platform;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
+    CustomIntegrator integrator(0.002);
+    integrator.addComputePerDof("v", "v+0.5*dt*f1/m");
+    for (int i = 0; i < 2; i++) {
+        integrator.addComputePerDof("v", "v+0.5*(dt/2)*f0/m");
+        integrator.addComputePerDof("x", "x+(dt/2)*v");
+        integrator.addComputePerDof("v", "v+0.5*(dt/2)*f0/m");
+    }
+    integrator.addComputePerDof("v", "v+0.5*dt*f1/m");
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    for (int i = 0; i < numParticles-2; i++)
+        bonds->addBond(i, i+1, 1.0, 0.5);
+    system.addForce(bonds);
+    NonbondedForce* nb = new NonbondedForce();
+    nb->setCutoffDistance(2.0);
+    nb->setNonbondedMethod(NonbondedForce::Ewald);
+    for (int i = 0; i < numParticles; ++i) {
+        system.addParticle(i%2 == 0 ? 5.0 : 10.0);
+        nb->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
+    }
+    nb->setForceGroup(1);
+    nb->setReciprocalSpaceForceGroup(0);
+    system.addForce(nb);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numParticles; ++i) {
+        positions[i] = Vec3(i/2, (i+1)/2, 0);
+        velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
+    }
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    
+    // Simulate it and monitor energy conservations.
+    
+    double initialEnergy = 0.0;
+    for (int i = 0; i < 1000; ++i) {
+        State state = context.getState(State::Energy);
+        double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+        if (i == 1)
+            initialEnergy = energy;
+        else if (i > 1)
+            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.05);
+        integrator.step(2);
+    }
+}
+
 int main() {
     try {
         testSingleBond();
@@ -521,6 +622,8 @@ int main() {
         testParameter();
         testRandomDistributions();
         testPerDofVariables();
+        testForceGroups();
+        testRespa();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

plugins/amoeba/openmmapi/include/openmm/internal/AmoebaGeneralizedKirkwoodForceImpl.h

@@ -53,7 +53,7 @@ public:
     void updateContextState(ContextImpl& context) {
         // This force field doesn't update the state directly.
     }
-    double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy);
+    double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups);
     std::map<std::string, double> getDefaultParameters() {
         return std::map<std::string, double>(); // This force field doesn't define any parameters.
     }

plugins/amoeba/openmmapi/include/openmm/internal/AmoebaHarmonicAngleForceImpl.h

@@ -56,7 +56,7 @@ public:
     void updateContextState(ContextImpl& context) {
         // This force field doesn't update the state directly.
     }
-    double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy);
+    double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups);
     std::map<std::string, double> getDefaultParameters() {
         return std::map<std::string, double>(); // This force field doesn't define any parameters.
     }

plugins/amoeba/openmmapi/include/openmm/internal/AmoebaHarmonicBondForceImpl.h

@@ -57,7 +57,7 @@ public:
     void updateContextState(ContextImpl& context) {
         // This force field doesn't update the state directly.
     }
-    double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy);
+    double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups);
     std::map<std::string, double> getDefaultParameters() {
         return std::map<std::string, double>(); // This force field doesn't define any parameters.
     }

plugins/amoeba/openmmapi/include/openmm/internal/AmoebaHarmonicInPlaneAngleForceImpl.h

@@ -56,7 +56,7 @@ public:
     void updateContextState(ContextImpl& context) {
         // This force field doesn't update the state directly.
     }
-    double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy);
+    double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups);
     std::map<std::string, double> getDefaultParameters() {
         return std::map<std::string, double>(); // This force field doesn't define any parameters.
     }

plugins/amoeba/openmmapi/include/openmm/internal/AmoebaMultipoleForceImpl.h

@@ -55,7 +55,7 @@ public:
     void updateContextState(ContextImpl& context) {
         // This force field doesn't update the state directly.
     }
-    double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy);
+    double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups);
     std::map<std::string, double> getDefaultParameters() {
         return std::map<std::string, double>(); // This force field doesn't define any parameters.
     }