OpenCL implementation of multiple time step integrators

platforms/cuda/src/CudaKernels.cpp

@@ -816,6 +816,8 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
             method = PERIODIC;
         }
         if (force.getNonbondedMethod() == NonbondedForce::Ewald || force.getNonbondedMethod() == NonbondedForce::PME) {
+            if (force.getReciprocalSpaceForceGroup() != 0)
+                throw OpenMMException("CudaPlatform does not support force groups");
             if (force.getNonbondedMethod() == NonbondedForce::Ewald) {
                 double alpha;
                 int kmaxx, kmaxy, kmaxz;

platforms/cuda/src/CudaPlatform.cpp

@@ -99,6 +99,9 @@ void CudaPlatform::contextCreated(ContextImpl& context, const map<string, string
     for (int i = 0; i < system.getNumParticles(); i++)
         if (system.isVirtualSite(i))
             throw OpenMMException("CudaPlatform does not support virtual sites");
+    for (int i = 0; i < system.getNumForces(); i++)
+        if (system.getForce(i).getForceGroup() != 0)
+            throw OpenMMException("CudaPlatform does not support force groups");
     unsigned int device = 0;
     const string& devicePropValue = (properties.find(CudaDevice()) == properties.end() ?
             getPropertyDefaultValue(CudaDevice()) : properties.find(CudaDevice())->second);

platforms/opencl/src/OpenCLBondedUtilities.cpp

@@ -43,10 +43,11 @@ OpenCLBondedUtilities::~OpenCLBondedUtilities() {
         delete bufferIndices[i];
 }
 
-void OpenCLBondedUtilities::addInteraction(const vector<vector<int> >& atoms, const string& source) {
+void OpenCLBondedUtilities::addInteraction(const vector<vector<int> >& atoms, const string& source, int group) {
     if (atoms.size() > 0) {
         forceAtoms.push_back(atoms);
         forceSource.push_back(source);
+        forceGroup.push_back(group);
         int width = 1;
         while (width < atoms[0].size())
             width *= 2;
@@ -154,7 +155,7 @@ void OpenCLBondedUtilities::initialize(const System& system) {
         const vector<int>& set = *iter;
         int setSize = set.size();
         stringstream s;
-        s<<"__kernel void computeBondedForces(__global float4* restrict forceBuffers, __global float* restrict energyBuffer, __global const float4* restrict posq";
+        s<<"__kernel void computeBondedForces(__global float4* restrict forceBuffers, __global float* restrict energyBuffer, __global const float4* restrict posq, int groups";
         for (int i = 0; i < setSize; i++) {
             int force = set[i];
             string indexType = "uint"+(indexWidth[force] == 1 ? "" : OpenCLExpressionUtilities::intToString(indexWidth[force]));
@@ -167,7 +168,7 @@ void OpenCLBondedUtilities::initialize(const System& system) {
         s<<"float energy = 0.0f;\n";
         for (int i = 0; i < setSize; i++) {
             int force = set[i];
-            s<<createForceSource(i, forceAtoms[force].size(), forceAtoms[force][0].size(), forceSource[force]);
+            s<<createForceSource(i, forceAtoms[force].size(), forceAtoms[force][0].size(), forceGroup[force], forceSource[force]);
         }
         s<<"energyBuffer[get_global_id(0)] += energy;\n";
         s<<"}\n";
@@ -180,7 +181,7 @@ void OpenCLBondedUtilities::initialize(const System& system) {
     forceSource.clear();
 }
 
-string OpenCLBondedUtilities::createForceSource(int forceIndex, int numBonds, int numAtoms, const string& computeForce) {
+string OpenCLBondedUtilities::createForceSource(int forceIndex, int numBonds, int numAtoms, int group, const string& computeForce) {
     maxBonds = max(maxBonds, numBonds);
     int width = 1;
     while (width < numAtoms)
@@ -198,6 +199,7 @@ string OpenCLBondedUtilities::createForceSource(int forceIndex, int numBonds, in
         suffix = suffix16;
     string indexType = "uint"+(width == 1 ? "" : OpenCLExpressionUtilities::intToString(width));
     stringstream s;
+    s<<"if ((groups&"<<(1<<group)<<") != 0)\n";
     s<<"for (unsigned int index = get_global_id(0); index < "<<numBonds<<"; index += get_global_size(0)) {\n";
     s<<"    "<<indexType<<" atoms = atomIndices"<<forceIndex<<"[index];\n";
     s<<"    "<<indexType<<" buffers = bufferIndices"<<forceIndex<<"[index];\n";
@@ -218,7 +220,7 @@ string OpenCLBondedUtilities::createForceSource(int forceIndex, int numBonds, in
     return s.str();
 }
 
-void OpenCLBondedUtilities::computeInteractions() {
+void OpenCLBondedUtilities::computeInteractions(int groups) {
     if (!hasInitializedKernels) {
         hasInitializedKernels = true;
         for (int i = 0; i < (int) forceSets.size(); i++) {
@@ -227,6 +229,7 @@ void OpenCLBondedUtilities::computeInteractions() {
             kernel.setArg<cl::Buffer>(index++, context.getForceBuffers().getDeviceBuffer());
             kernel.setArg<cl::Buffer>(index++, context.getEnergyBuffer().getDeviceBuffer());
             kernel.setArg<cl::Buffer>(index++, context.getPosq().getDeviceBuffer());
+            index++;
             for (int j = 0; j < (int) forceSets[i].size(); j++) {
                 kernel.setArg<cl::Buffer>(index++, atomIndices[forceSets[i][j]]->getDeviceBuffer());
                 kernel.setArg<cl::Buffer>(index++, bufferIndices[forceSets[i][j]]->getDeviceBuffer());
@@ -235,6 +238,8 @@ void OpenCLBondedUtilities::computeInteractions() {
                 kernel.setArg<cl::Memory>(index++, *arguments[j]);
         }
     }
-    for (int i = 0; i < (int) kernels.size(); i++)
+    for (int i = 0; i < (int) kernels.size(); i++) {
+        kernels[i].setArg<cl_int>(3, groups);
         context.executeKernel(kernels[i], maxBonds);
+    }
 }

platforms/opencl/src/OpenCLBondedUtilities.h

@@ -88,8 +88,9 @@ public:
      * @param atoms    this should have one entry for each bond, and that entry should contain the list
      *                 of atoms involved in the bond.  Every entry must have the same number of atoms.
      * @param source   the code to evaluate the interaction
+     * @param group    the force group in which the interaction should be calculated
      */
-    void addInteraction(const std::vector<std::vector<int> >& atoms, const std::string& source);
+    void addInteraction(const std::vector<std::vector<int> >& atoms, const std::string& source, int group);
     /**
      * Add an argument that should be passed to the interaction kernel.
      * 
@@ -111,15 +112,18 @@ public:
     }
     /**
      * Compute the bonded interactions.
+     * 
+     * @param groups        a set of bit flags for which force groups to include
      */
-    void computeInteractions();
+    void computeInteractions(int groups);
 private:
-    std::string createForceSource(int forceIndex, int numBonds, int numAtoms, const std::string& computeForce);
+    std::string createForceSource(int forceIndex, int numBonds, int numAtoms, int group, const std::string& computeForce);
     OpenCLContext& context;
     std::vector<cl::Kernel> kernels;
     std::vector<std::vector<std::vector<int> > > forceAtoms;
     std::vector<int> indexWidth;
     std::vector<std::string> forceSource;
+    std::vector<int> forceGroup;
     std::vector<std::vector<int> > forceSets;
     std::vector<cl::Memory*> arguments;
     std::vector<std::string> argTypes;

platforms/opencl/src/OpenCLKernels.cpp

@@ -94,18 +94,22 @@ void OpenCLCalcForcesAndEnergyKernel::initialize(const System& system) {
 }
 
 void OpenCLCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
+    OpenCLNonbondedUtilities& nb = cl.getNonbondedUtilities();
+    bool includeNonbonded = ((groups&(1<<nb.getForceGroup())) != 0);
     cl.setAtomsWereReordered(false);
-    if (cl.getNonbondedUtilities().getUseCutoff() && cl.getComputeForceCount()%100 == 0) {
+    if (nb.getUseCutoff() && includeNonbonded && cl.getComputeForceCount()%100 == 0) {
         cl.reorderAtoms();
-        cl.getNonbondedUtilities().updateNeighborListSize();
-    }
+        nb.updateNeighborListSize();
         cl.setComputeForceCount(cl.getComputeForceCount()+1);
+    }
     cl.clearAutoclearBuffers();
-    cl.getNonbondedUtilities().prepareInteractions();
+    if (includeNonbonded)
+        nb.prepareInteractions();
 }
 
 double OpenCLCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
-    cl.getBondedUtilities().computeInteractions();
+    cl.getBondedUtilities().computeInteractions(groups);
+    if ((groups&(1<<cl.getNonbondedUtilities().getForceGroup())) != 0)
         cl.getNonbondedUtilities().computeInteractions();
     cl.reduceForces();
     cl.getIntegrationUtilities().distributeForcesFromVirtualSites();
@@ -293,7 +297,7 @@ void OpenCLCalcHarmonicBondForceKernel::initialize(const System& system, const H
     params->upload(paramVector);
     map<string, string> replacements;
     replacements["PARAMS"] = cl.getBondedUtilities().addArgument(params->getDeviceBuffer(), "float2");
-    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::harmonicBondForce, replacements));
+    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::harmonicBondForce, replacements), force.getForceGroup());
     cl.addForce(new OpenCLBondForceInfo(0, force));
 }
 
@@ -399,7 +403,7 @@ void OpenCLCalcCustomBondForceKernel::initialize(const System& system, const Cus
     compute << OpenCLExpressionUtilities::createExpressions(expressions, variables, functions, "temp", "");
     map<string, string> replacements;
     replacements["COMPUTE_FORCE"] = compute.str();
-    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::customBondForce, replacements));
+    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::customBondForce, replacements), force.getForceGroup());
 }
 
 double OpenCLCalcCustomBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -468,7 +472,7 @@ void OpenCLCalcHarmonicAngleForceKernel::initialize(const System& system, const
     params->upload(paramVector);
     map<string, string> replacements;
     replacements["PARAMS"] = cl.getBondedUtilities().addArgument(params->getDeviceBuffer(), "float2");
-    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::harmonicAngleForce, replacements));
+    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::harmonicAngleForce, replacements), force.getForceGroup());
     cl.addForce(new OpenCLAngleForceInfo(0, force));
 }
 
@@ -575,7 +579,7 @@ void OpenCLCalcCustomAngleForceKernel::initialize(const System& system, const Cu
     compute << OpenCLExpressionUtilities::createExpressions(expressions, variables, functions, "temp", "");
     map<string, string> replacements;
     replacements["COMPUTE_FORCE"] = compute.str();
-    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::customAngleForce, replacements));
+    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::customAngleForce, replacements), force.getForceGroup());
 }
 
 double OpenCLCalcCustomAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -646,7 +650,7 @@ void OpenCLCalcPeriodicTorsionForceKernel::initialize(const System& system, cons
     params->upload(paramVector);
     map<string, string> replacements;
     replacements["PARAMS"] = cl.getBondedUtilities().addArgument(params->getDeviceBuffer(), "float4");
-    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::periodicTorsionForce, replacements));
+    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::periodicTorsionForce, replacements), force.getForceGroup());
     cl.addForce(new OpenCLPeriodicTorsionForceInfo(0, force));
 }
 
@@ -706,7 +710,7 @@ void OpenCLCalcRBTorsionForceKernel::initialize(const System& system, const RBTo
     params->upload(paramVector);
     map<string, string> replacements;
     replacements["PARAMS"] = cl.getBondedUtilities().addArgument(params->getDeviceBuffer(), "float8");
-    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::rbTorsionForce, replacements));
+    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::rbTorsionForce, replacements), force.getForceGroup());
     cl.addForce(new OpenCLRBTorsionForceInfo(0, force));
 }
 
@@ -793,7 +797,7 @@ void OpenCLCalcCMAPTorsionForceKernel::initialize(const System& system, const CM
     replacements["COEFF"] = cl.getBondedUtilities().addArgument(coefficients->getDeviceBuffer(), "float4");
     replacements["MAP_POS"] = cl.getBondedUtilities().addArgument(mapPositions->getDeviceBuffer(), "int2");
     replacements["MAPS"] = cl.getBondedUtilities().addArgument(torsionMaps->getDeviceBuffer(), "int");
-    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::cmapTorsionForce, replacements));
+    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::cmapTorsionForce, replacements), force.getForceGroup());
     cl.addForce(new OpenCLCMAPTorsionForceInfo(0, force));
 }
 
@@ -902,7 +906,7 @@ void OpenCLCalcCustomTorsionForceKernel::initialize(const System& system, const
     map<string, string> replacements;
     replacements["COMPUTE_FORCE"] = compute.str();
     replacements["M_PI"] = doubleToString(M_PI);
-    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::customTorsionForce, replacements));
+    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::customTorsionForce, replacements), force.getForceGroup());
 }
 
 double OpenCLCalcCustomTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -1170,7 +1174,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
     // Add the interaction to the default nonbonded kernel.
     
     string source = cl.replaceStrings(OpenCLKernelSources::coulombLennardJones, defines);
-    cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source);
+    cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source, force.getForceGroup());
     if (hasLJ)
         cl.getNonbondedUtilities().addParameter(OpenCLNonbondedUtilities::ParameterInfo("sigmaEpsilon", "float", 2, sizeof(cl_float2), sigmaEpsilon->getDeviceBuffer()));
 
@@ -1192,7 +1196,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
         exceptionParams->upload(exceptionParamsVector);
         map<string, string> replacements;
         replacements["PARAMS"] = cl.getBondedUtilities().addArgument(exceptionParams->getDeviceBuffer(), "float4");
-        cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::nonbondedExceptions, replacements));
+        cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::nonbondedExceptions, replacements), force.getForceGroup());
     }
     cl.addForce(new OpenCLNonbondedForceInfo(cl.getNonbondedUtilities().getNumForceBuffers(), force));
 }
@@ -1252,7 +1256,7 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ
             }
        }
     }
-    if (cosSinSums != NULL && cl.getContextIndex() == 0) {
+    if (cosSinSums != NULL && cl.getContextIndex() == 0 && includeReciprocal) {
         mm_float4 boxSize = cl.getPeriodicBoxSize();
         mm_float4 recipBoxSize = mm_float4((float) (2*M_PI/boxSize.x), (float) (2*M_PI/boxSize.y), (float) (2*M_PI/boxSize.z), 0);
         float recipCoefficient = (float) (ONE_4PI_EPS0*4*M_PI/(boxSize.x*boxSize.y*boxSize.z));
@@ -1263,7 +1267,7 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ
         ewaldForcesKernel.setArg<cl_float>(4, recipCoefficient);
         cl.executeKernel(ewaldForcesKernel, cl.getNumAtoms());
     }
-    if (pmeGrid != NULL && cl.getContextIndex() == 0) {
+    if (pmeGrid != NULL && cl.getContextIndex() == 0 && includeReciprocal) {
         mm_float4 boxSize = cl.getPeriodicBoxSize();
         mm_float4 invBoxSize = cl.getInvPeriodicBoxSize();
         pmeUpdateBsplinesKernel.setArg<mm_float4>(4, boxSize);
@@ -1443,7 +1447,7 @@ void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, cons
     map<string, string> replacements;
     replacements["COMPUTE_FORCE"] = compute.str();
     string source = cl.replaceStrings(OpenCLKernelSources::customNonbonded, replacements);
-    cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source);
+    cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source, force.getForceGroup());
     for (int i = 0; i < (int) params->getBuffers().size(); i++) {
         const OpenCLNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
         cl.getNonbondedUtilities().addParameter(OpenCLNonbondedUtilities::ParameterInfo(prefix+"params"+intToString(i+1), buffer.getComponentType(), buffer.getNumComponents(), buffer.getSize(), buffer.getMemory()));
@@ -1538,7 +1542,7 @@ void OpenCLCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOB
     bool useCutoff = (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff);
     bool usePeriodic = (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff && force.getNonbondedMethod() != GBSAOBCForce::CutoffNonPeriodic);
     string source = OpenCLKernelSources::gbsaObc2;
-    nb.addInteraction(useCutoff, usePeriodic, false, force.getCutoffDistance(), vector<vector<int> >(), source);
+    nb.addInteraction(useCutoff, usePeriodic, false, force.getCutoffDistance(), vector<vector<int> >(), source, force.getForceGroup());
     nb.addParameter(OpenCLNonbondedUtilities::ParameterInfo("obcParams", "float", 2, sizeof(cl_float2), params->getDeviceBuffer()));;
     nb.addParameter(OpenCLNonbondedUtilities::ParameterInfo("bornForce", "float", 1, sizeof(cl_float), bornForce->getDeviceBuffer()));;
     cl.addForce(new OpenCLGBSAOBCForceInfo(nb.getNumForceBuffers(), force));
@@ -2351,7 +2355,7 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
             globals->upload(globalParamValues);
             arguments.push_back(OpenCLNonbondedUtilities::ParameterInfo(prefix+"globals", "float", 1, sizeof(cl_float), globals->getDeviceBuffer()));
         }
-        cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, force.getNumExclusions() > 0, force.getCutoffDistance(), exclusionList, source);
+        cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, force.getNumExclusions() > 0, force.getCutoffDistance(), exclusionList, source, force.getForceGroup());
         for (int i = 0; i < (int) parameters.size(); i++)
             cl.getNonbondedUtilities().addParameter(parameters[i]);
         for (int i = 0; i < (int) arguments.size(); i++)
@@ -2660,7 +2664,7 @@ void OpenCLCalcCustomExternalForceKernel::initialize(const System& system, const
     compute << OpenCLExpressionUtilities::createExpressions(expressions, variables, functions, "temp", "");
     map<string, string> replacements;
     replacements["COMPUTE_FORCE"] = compute.str();
-    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::customExternalForce, replacements));
+    cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::customExternalForce, replacements), force.getForceGroup());
 }
 
 double OpenCLCalcCustomExternalForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -3536,6 +3540,8 @@ public:
     void execute() {
         // Reorder the per-DOF variables to reflect the new atom order.
 
+        if (perDofValues.getNumParameters() == 0)
+            return;
         int numAtoms = cl.getNumAtoms();
         if (deviceValuesAreCurrent)
             perDofValues.getParameterValues(localPerDofValues);
@@ -3623,7 +3629,7 @@ string OpenCLIntegrateCustomStepKernel::createGlobalComputation(const string& va
     return OpenCLExpressionUtilities::createExpressions(expressions, variables, functions, "temp", "");
 }
 
-string OpenCLIntegrateCustomStepKernel::createPerDofComputation(const string& variable, const Lepton::ParsedExpression& expr, int component, CustomIntegrator& integrator) {
+string OpenCLIntegrateCustomStepKernel::createPerDofComputation(const string& variable, const Lepton::ParsedExpression& expr, int component, CustomIntegrator& integrator, const string& forceName) {
     const string suffixes[] = {".x", ".y", ".z"};
     string suffix = suffixes[component];
     map<string, Lepton::ParsedExpression> expressions;
@@ -3643,7 +3649,7 @@ string OpenCLIntegrateCustomStepKernel::createPerDofComputation(const string& va
     map<string, string> variables;
     variables["x"] = "position"+suffix;
     variables["v"] = "velocity"+suffix;
-    variables["f"] = "f"+suffix;
+    variables[forceName] = "f"+suffix;
     variables["gaussian"] = "gaussian"+suffix;
     variables["uniform"] = "uniform"+suffix;
     variables["m"] = "mass";
@@ -3681,6 +3687,7 @@ void OpenCLIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegr
         requiredUniform.resize(integrator.getNumComputations(), 0);
         needsForces.resize(numSteps, false);
         needsEnergy.resize(numSteps, false);
+        forceGroup.resize(numSteps, -2);
         invalidatesForces.resize(numSteps, false);
         merged.resize(numSteps, false);
         modifiesParameters = false;
@@ -3722,15 +3729,39 @@ void OpenCLIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegr
         stepType.resize(numSteps);
         vector<string> variable(numSteps);
         vector<Lepton::ParsedExpression> expression(numSteps);
+        vector<string> forceGroupName;
+        for (int i = 0; i < 32; i++) {
+            stringstream str;
+            str << "f" << i;
+            forceGroupName.push_back(str.str());
+        }
+        vector<string> forceName(numSteps, "f");
         for (int step = 0; step < numSteps; step++) {
             string expr;
             integrator.getComputationStep(step, stepType[step], variable[step], expr);
             if (expr.size() > 0) {
                 expression[step] = Lepton::Parser::parse(expr).optimize();
-                needsForces[step] = usesVariable(expression[step], "f");
-                needsEnergy[step] = usesVariable(expression[step], "energy");
+                if (usesVariable(expression[step], "f")) {
+                    needsForces[step] = true;
+                    forceGroup[step] = -1;
+                }
+                if (usesVariable(expression[step], "energy")) {
+                    needsEnergy[step] = true;
+                    forceGroup[step] = -1;
+                }
+                for (int i = 0; i < 32; i++) {
+                    if (usesVariable(expression[step], forceGroupName[i])) {
+                        if (forceGroup[step] != -2)
+                            throw OpenMMException("A single computation step cannot depend on multiple force groups");
+                        needsForces[step] = true;
+                        forceGroup[step] = 1<<i;
+                        forceName[step] = forceGroupName[i];
+                    }
+                }
             }
             invalidatesForces[step] = (stepType[step] == CustomIntegrator::ConstrainPositions || affectsForce.find(variable[step]) != affectsForce.end());
+            if (forceGroup[step] == -2 && step > 0)
+                forceGroup[step] = forceGroup[step-1];
         }
         
         // Determine how each step will represent the position (as just a value, or a value plus a delta).
@@ -3756,7 +3787,7 @@ void OpenCLIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegr
         // Identify steps that can be merged into a single kernel.
         
         for (int step = 1; step < numSteps; step++) {
-            if ((needsForces[step] || needsEnergy[step]) && invalidatesForces[step-1])
+            if ((needsForces[step] || needsEnergy[step]) && (invalidatesForces[step-1] || forceGroup[step] != forceGroup[step-1]))
                 continue;
             if (stepType[step-1] == CustomIntegrator::ComputeGlobal && stepType[step] == CustomIntegrator::ComputeGlobal)
                 merged[step] = true;
@@ -3784,7 +3815,7 @@ void OpenCLIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegr
                 for (int j = step; j < numSteps && (j == step || merged[j]); j++) {
                     compute << "{\n";
                     for (int i = 0; i < 3; i++)
-                        compute << createPerDofComputation(stepType[j] == CustomIntegrator::ComputePerDof ? variable[j] : "", expression[j], i, integrator);
+                        compute << createPerDofComputation(stepType[j] == CustomIntegrator::ComputePerDof ? variable[j] : "", expression[j], i, integrator, forceName[j]);
                     if (variable[j] == "x") {
                         if (storePosAsDelta[j]) {
                             if (cl.getSupportsDoublePrecision())
@@ -3934,7 +3965,7 @@ void OpenCLIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegr
     // Loop over computation steps in the integrator and execute them.
 
     for (int i = 0; i < numSteps; i++) {
-        if ((needsForces[i] || needsEnergy[i]) && !forcesAreValid) {
+        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.
             
@@ -3952,7 +3983,7 @@ void OpenCLIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegr
                     break;
             }
             recordChangedParameters(context);
-            context.calcForcesAndEnergy(computeForce, false);
+            context.calcForcesAndEnergy(computeForce, false, forceGroup[i]);
             if (computeEnergy)
                 cl.executeKernel(sumEnergyKernel, OpenCLContext::ThreadBlockSize, OpenCLContext::ThreadBlockSize);
             forcesAreValid = true;

platforms/opencl/src/OpenCLKernels.h

@@ -978,7 +978,7 @@ public:
 private:
     class ReorderListener;
     std::string createGlobalComputation(const std::string& variable, const Lepton::ParsedExpression& expr, CustomIntegrator& integrator);
-    std::string createPerDofComputation(const std::string& variable, const Lepton::ParsedExpression& expr, int component, CustomIntegrator& integrator);
+    std::string createPerDofComputation(const std::string& variable, const Lepton::ParsedExpression& expr, int component, CustomIntegrator& integrator, const std::string& forceName);
     void recordChangedParameters(ContextImpl& context);
     OpenCLContext& cl;
     double prevStepSize;
@@ -1000,6 +1000,7 @@ private:
     std::vector<bool> needsEnergy;
     std::vector<bool> invalidatesForces;
     std::vector<bool> merged;
+    std::vector<int> forceGroup;
     std::vector<int> requiredGaussian;
     std::vector<int> requiredUniform;
     std::vector<std::string> parameterNames;

platforms/opencl/src/OpenCLNonbondedUtilities.cpp

@@ -37,7 +37,7 @@ using namespace std;
 
 OpenCLNonbondedUtilities::OpenCLNonbondedUtilities(OpenCLContext& context) : context(context), cutoff(-1.0), useCutoff(false), anyExclusions(false),
         numForceBuffers(0), exclusionIndices(NULL), exclusionRowIndices(NULL), exclusions(NULL), interactingTiles(NULL), interactionFlags(NULL),
-        interactionCount(NULL), blockCenter(NULL), blockBoundingBox(NULL) {
+        interactionCount(NULL), blockCenter(NULL), blockBoundingBox(NULL), nonbondedForceGroup(0) {
     // Decide how many thread blocks and force buffers to use.
 
     deviceIsCpu = (context.getDevice().getInfo<CL_DEVICE_TYPE>() == CL_DEVICE_TYPE_CPU);
@@ -91,7 +91,7 @@ OpenCLNonbondedUtilities::~OpenCLNonbondedUtilities() {
         delete blockBoundingBox;
 }
 
-void OpenCLNonbondedUtilities::addInteraction(bool usesCutoff, bool usesPeriodic, bool usesExclusions, double cutoffDistance, const vector<vector<int> >& exclusionList, const string& kernel) {
+void OpenCLNonbondedUtilities::addInteraction(bool usesCutoff, bool usesPeriodic, bool usesExclusions, double cutoffDistance, const vector<vector<int> >& exclusionList, const string& kernel, int forceGroup) {
     if (cutoff != -1.0) {
         if (usesCutoff != useCutoff)
             throw OpenMMException("All Forces must agree on whether to use a cutoff");
@@ -99,6 +99,8 @@ void OpenCLNonbondedUtilities::addInteraction(bool usesCutoff, bool usesPeriodic
             throw OpenMMException("All Forces must agree on whether to use periodic boundary conditions");
         if (cutoffDistance != cutoff)
             throw OpenMMException("All Forces must use the same cutoff distance");
+        if (forceGroup != nonbondedForceGroup)
+            throw OpenMMException("All nonbonded forces must be in the same force group");
     }
     if (usesExclusions)
         requestExclusions(exclusionList);
@@ -106,6 +108,7 @@ void OpenCLNonbondedUtilities::addInteraction(bool usesCutoff, bool usesPeriodic
     usePeriodic = usesPeriodic;
     cutoff = cutoffDistance;
     kernelSource += kernel+"\n";
+    nonbondedForceGroup = forceGroup;
 }
 
 void OpenCLNonbondedUtilities::addParameter(const ParameterInfo& parameter) {

platforms/opencl/src/OpenCLNonbondedUtilities.h

@@ -74,8 +74,9 @@ public:
      * @param cutoffDistance the cutoff distance for this interaction (ignored if usesCutoff is false)
      * @param exclusionList  for each atom, specifies the list of other atoms whose interactions should be excluded
      * @param kernel         the code to evaluate the interaction
+     * @param forceGroup     the force group in which the interaction should be calculated
      */
-    void addInteraction(bool usesCutoff, bool usesPeriodic, bool usesExclusions, double cutoffDistance, const std::vector<std::vector<int> >& exclusionList, const std::string& kernel);
+    void addInteraction(bool usesCutoff, bool usesPeriodic, bool usesExclusions, double cutoffDistance, const std::vector<std::vector<int> >& exclusionList, const std::string& kernel, int forceGroup);
     /**
      * Add a per-atom parameter that the default interaction kernel may depend on.
      */
@@ -148,6 +149,12 @@ public:
     bool getHasInteractions() {
         return cutoff != -1.0;
     }
+    /**
+     * Get the force group in which nonbonded interactions should be computed.
+     */
+    int getForceGroup() {
+        return nonbondedForceGroup;
+    }
     /**
      * Prepare to compute interactions.  This updates the neighbor list.
      */
@@ -258,7 +265,7 @@ private:
     std::map<std::string, std::string> kernelDefines;
     double cutoff;
     bool useCutoff, usePeriodic, forceBufferPerAtomBlock, deviceIsCpu, anyExclusions;
-    int numForceBuffers, startTileIndex, numTiles, numForceThreadBlocks, forceThreadBlockSize;
+    int numForceBuffers, startTileIndex, numTiles, numForceThreadBlocks, forceThreadBlockSize, nonbondedForceGroup;
 };
 
 /**

platforms/opencl/tests/TestOpenCLCustomIntegrator.cpp

@@ -525,6 +525,107 @@ void testPerDofVariables() {
     }
 }
 
+/**
+ * Test evaluating force groups separately.
+ */
+void testForceGroups() {
+    OpenCLPlatform 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;
+    OpenCLPlatform 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();
@@ -536,6 +637,8 @@ int main() {
         testParameter();
         testRandomDistributions();
         testPerDofVariables();
+        testForceGroups();
+        testRespa();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;