Merged parallel code (#4649)

* Unified lots of parallel computation code between platforms

* Unified test code between platforms

* Eliminated duplicated timing code

platforms/opencl/src/OpenCLContext.cpp

@@ -658,6 +658,10 @@ vector<ComputeContext*> OpenCLContext::getAllContexts() {
     return result;
 }
 
+double& OpenCLContext::getEnergyWorkspace() {
+    return platformData.contextEnergy[contextIndex];
+}
+
 cl::CommandQueue& OpenCLContext::getQueue() {
     return currentQueue;
 }

platforms/opencl/src/OpenCLKernelFactory.cpp

@@ -27,6 +27,7 @@
 #include "OpenCLKernelFactory.h"
 #include "OpenCLParallelKernels.h"
 #include "openmm/common/CommonKernels.h"
+#include "openmm/common/CommonParallelKernels.h"
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/OpenMMException.h"
 
@@ -34,39 +35,39 @@ using namespace OpenMM;
 
 KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const {
     OpenCLPlatform::PlatformData& data = *static_cast<OpenCLPlatform::PlatformData*>(context.getPlatformData());
+    OpenCLContext& cl = *data.contexts[0];
     if (data.contexts.size() > 1) {
         // We are running in parallel on multiple devices, so we may want to create a parallel kernel.
         
         if (name == CalcForcesAndEnergyKernel::Name())
             return new OpenCLParallelCalcForcesAndEnergyKernel(name, platform, data);
         if (name == CalcHarmonicBondForceKernel::Name())
-            return new OpenCLParallelCalcHarmonicBondForceKernel(name, platform, data, context.getSystem());
+            return new CommonParallelCalcHarmonicBondForceKernel(name, platform, cl, context.getSystem());
         if (name == CalcCustomBondForceKernel::Name())
-            return new OpenCLParallelCalcCustomBondForceKernel(name, platform, data, context.getSystem());
+            return new CommonParallelCalcCustomBondForceKernel(name, platform, cl, context.getSystem());
         if (name == CalcHarmonicAngleForceKernel::Name())
-            return new OpenCLParallelCalcHarmonicAngleForceKernel(name, platform, data, context.getSystem());
+            return new CommonParallelCalcHarmonicAngleForceKernel(name, platform, cl, context.getSystem());
         if (name == CalcCustomAngleForceKernel::Name())
-            return new OpenCLParallelCalcCustomAngleForceKernel(name, platform, data, context.getSystem());
+            return new CommonParallelCalcCustomAngleForceKernel(name, platform, cl, context.getSystem());
         if (name == CalcPeriodicTorsionForceKernel::Name())
-            return new OpenCLParallelCalcPeriodicTorsionForceKernel(name, platform, data, context.getSystem());
+            return new CommonParallelCalcPeriodicTorsionForceKernel(name, platform, cl, context.getSystem());
         if (name == CalcRBTorsionForceKernel::Name())
-            return new OpenCLParallelCalcRBTorsionForceKernel(name, platform, data, context.getSystem());
+            return new CommonParallelCalcRBTorsionForceKernel(name, platform, cl, context.getSystem());
         if (name == CalcCMAPTorsionForceKernel::Name())
-            return new OpenCLParallelCalcCMAPTorsionForceKernel(name, platform, data, context.getSystem());
+            return new CommonParallelCalcCMAPTorsionForceKernel(name, platform, cl, context.getSystem());
         if (name == CalcCustomTorsionForceKernel::Name())
-            return new OpenCLParallelCalcCustomTorsionForceKernel(name, platform, data, context.getSystem());
+            return new CommonParallelCalcCustomTorsionForceKernel(name, platform, cl, context.getSystem());
         if (name == CalcNonbondedForceKernel::Name())
             return new OpenCLParallelCalcNonbondedForceKernel(name, platform, data, context.getSystem());
         if (name == CalcCustomNonbondedForceKernel::Name())
-            return new OpenCLParallelCalcCustomNonbondedForceKernel(name, platform, data, context.getSystem());
+            return new CommonParallelCalcCustomNonbondedForceKernel(name, platform, cl, context.getSystem());
         if (name == CalcCustomExternalForceKernel::Name())
-            return new OpenCLParallelCalcCustomExternalForceKernel(name, platform, data, context.getSystem());
+            return new CommonParallelCalcCustomExternalForceKernel(name, platform, cl, context.getSystem());
         if (name == CalcCustomHbondForceKernel::Name())
-            return new OpenCLParallelCalcCustomHbondForceKernel(name, platform, data, context.getSystem());
+            return new CommonParallelCalcCustomHbondForceKernel(name, platform, cl, context.getSystem());
         if (name == CalcCustomCompoundBondForceKernel::Name())
-            return new OpenCLParallelCalcCustomCompoundBondForceKernel(name, platform, data, context.getSystem());
+            return new CommonParallelCalcCustomCompoundBondForceKernel(name, platform, cl, context.getSystem());
     }
-    OpenCLContext& cl = *data.contexts[0];
     if (name == CalcForcesAndEnergyKernel::Name())
         return new OpenCLCalcForcesAndEnergyKernel(name, platform, cl);
     if (name == UpdateStateDataKernel::Name())

platforms/opencl/src/OpenCLParallelKernels.cpp

@@ -25,32 +25,11 @@
  * -------------------------------------------------------------------------- */
 
 #include "OpenCLParallelKernels.h"
+#include "openmm/internal/timer.h"
 
 using namespace OpenMM;
 using namespace std;
 
-/**
- * Get the current clock time, measured in microseconds.
- */
-#ifdef _MSC_VER
-    #include <Windows.h>
-    static long long getTime() {
-        FILETIME ft;
-        GetSystemTimeAsFileTime(&ft); // 100-nanoseconds since 1-1-1601
-        ULARGE_INTEGER result;
-        result.LowPart = ft.dwLowDateTime;
-        result.HighPart = ft.dwHighDateTime;
-        return result.QuadPart/10;
-    }
-#else
-    #include <sys/time.h> 
-    static long long getTime() {
-        struct timeval tod;
-        gettimeofday(&tod, 0);
-        return 1000000*tod.tv_sec+tod.tv_usec;
-    }
-#endif
-
 class OpenCLParallelCalcForcesAndEnergyKernel::BeginComputationTask : public OpenCLContext::WorkTask {
 public:
     BeginComputationTask(ContextImpl& context, OpenCLContext& cl, OpenCLCalcForcesAndEnergyKernel& kernel,
@@ -79,7 +58,7 @@ private:
 class OpenCLParallelCalcForcesAndEnergyKernel::FinishComputationTask : public OpenCLContext::WorkTask {
 public:
     FinishComputationTask(ContextImpl& context, OpenCLContext& cl, OpenCLCalcForcesAndEnergyKernel& kernel,
-            bool includeForce, bool includeEnergy, int groups, double& energy, long long& completionTime, void* pinnedMemory, bool& valid, int& numTiles) :
+            bool includeForce, bool includeEnergy, int groups, double& energy, double& completionTime, void* pinnedMemory, bool& valid, int& numTiles) :
             context(context), cl(cl), kernel(kernel), includeForce(includeForce), includeEnergy(includeEnergy), groups(groups), energy(energy),
             completionTime(completionTime), pinnedMemory(pinnedMemory), valid(valid), numTiles(numTiles) {
     }
@@ -97,7 +76,7 @@ public:
             else
                 cl.getQueue().finish();
         }
-        completionTime = getTime();
+        completionTime = getCurrentTime();
         if (cl.getNonbondedUtilities().getUsePeriodic() && numTiles > cl.getNonbondedUtilities().getInteractingTiles().getSize()) {
             valid = false;
             cl.getNonbondedUtilities().updateNeighborListSize();
@@ -110,7 +89,7 @@ private:
     bool includeForce, includeEnergy;
     int groups;
     double& energy;
-    long long& completionTime;
+    double& completionTime;
     void* pinnedMemory;
     bool& valid;
     int& numTiles;
@@ -211,334 +190,6 @@ double OpenCLParallelCalcForcesAndEnergyKernel::finishComputation(ContextImpl& c
     return energy;
 }
 
-class OpenCLParallelCalcHarmonicBondForceKernel::Task : public OpenCLContext::WorkTask {
-public:
-    Task(ContextImpl& context, CommonCalcHarmonicBondForceKernel& kernel, bool includeForce,
-            bool includeEnergy, double& energy) : context(context), kernel(kernel),
-            includeForce(includeForce), includeEnergy(includeEnergy), energy(energy) {
-    }
-    void execute() {
-        energy += kernel.execute(context, includeForce, includeEnergy);
-    }
-private:
-    ContextImpl& context;
-    CommonCalcHarmonicBondForceKernel& kernel;
-    bool includeForce, includeEnergy;
-    double& energy;
-};
-
-OpenCLParallelCalcHarmonicBondForceKernel::OpenCLParallelCalcHarmonicBondForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, const System& system) :
-        CalcHarmonicBondForceKernel(name, platform), data(data) {
-    for (int i = 0; i < (int) data.contexts.size(); i++)
-        kernels.push_back(Kernel(new CommonCalcHarmonicBondForceKernel(name, platform, *data.contexts[i], system)));
-}
-
-void OpenCLParallelCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).initialize(system, force);
-}
-
-double OpenCLParallelCalcHarmonicBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    for (int i = 0; i < (int) data.contexts.size(); i++) {
-        OpenCLContext& cl = *data.contexts[i];
-        ComputeContext::WorkThread& thread = cl.getWorkThread();
-        thread.addTask(new Task(context, getKernel(i), includeForces, includeEnergy, data.contextEnergy[i]));
-    }
-    return 0.0;
-}
-
-void OpenCLParallelCalcHarmonicBondForceKernel::copyParametersToContext(ContextImpl& context, const HarmonicBondForce& force, int firstBond, int lastBond) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).copyParametersToContext(context, force, firstBond, lastBond);
-}
-
-class OpenCLParallelCalcCustomBondForceKernel::Task : public OpenCLContext::WorkTask {
-public:
-    Task(ContextImpl& context, CommonCalcCustomBondForceKernel& kernel, bool includeForce,
-            bool includeEnergy, double& energy) : context(context), kernel(kernel),
-            includeForce(includeForce), includeEnergy(includeEnergy), energy(energy) {
-    }
-    void execute() {
-        energy += kernel.execute(context, includeForce, includeEnergy);
-    }
-private:
-    ContextImpl& context;
-    CommonCalcCustomBondForceKernel& kernel;
-    bool includeForce, includeEnergy;
-    double& energy;
-};
-
-OpenCLParallelCalcCustomBondForceKernel::OpenCLParallelCalcCustomBondForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, const System& system) :
-        CalcCustomBondForceKernel(name, platform), data(data) {
-    for (int i = 0; i < (int) data.contexts.size(); i++)
-        kernels.push_back(Kernel(new CommonCalcCustomBondForceKernel(name, platform, *data.contexts[i], system)));
-}
-
-void OpenCLParallelCalcCustomBondForceKernel::initialize(const System& system, const CustomBondForce& force) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).initialize(system, force);
-}
-
-double OpenCLParallelCalcCustomBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    for (int i = 0; i < (int) data.contexts.size(); i++) {
-        OpenCLContext& cl = *data.contexts[i];
-        ComputeContext::WorkThread& thread = cl.getWorkThread();
-        thread.addTask(new Task(context, getKernel(i), includeForces, includeEnergy, data.contextEnergy[i]));
-    }
-    return 0.0;
-}
-
-void OpenCLParallelCalcCustomBondForceKernel::copyParametersToContext(ContextImpl& context, const CustomBondForce& force, int firstBond, int lastBond) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).copyParametersToContext(context, force, firstBond, lastBond);
-}
-
-class OpenCLParallelCalcHarmonicAngleForceKernel::Task : public OpenCLContext::WorkTask {
-public:
-    Task(ContextImpl& context, CommonCalcHarmonicAngleForceKernel& kernel, bool includeForce,
-            bool includeEnergy, double& energy) : context(context), kernel(kernel),
-            includeForce(includeForce), includeEnergy(includeEnergy), energy(energy) {
-    }
-    void execute() {
-        energy += kernel.execute(context, includeForce, includeEnergy);
-    }
-private:
-    ContextImpl& context;
-    CommonCalcHarmonicAngleForceKernel& kernel;
-    bool includeForce, includeEnergy;
-    double& energy;
-};
-
-OpenCLParallelCalcHarmonicAngleForceKernel::OpenCLParallelCalcHarmonicAngleForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, const System& system) :
-        CalcHarmonicAngleForceKernel(name, platform), data(data) {
-    for (int i = 0; i < (int) data.contexts.size(); i++)
-        kernels.push_back(Kernel(new CommonCalcHarmonicAngleForceKernel(name, platform, *data.contexts[i], system)));
-}
-
-void OpenCLParallelCalcHarmonicAngleForceKernel::initialize(const System& system, const HarmonicAngleForce& force) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).initialize(system, force);
-}
-
-double OpenCLParallelCalcHarmonicAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    for (int i = 0; i < (int) data.contexts.size(); i++) {
-        OpenCLContext& cl = *data.contexts[i];
-        ComputeContext::WorkThread& thread = cl.getWorkThread();
-        thread.addTask(new Task(context, getKernel(i), includeForces, includeEnergy, data.contextEnergy[i]));
-    }
-    return 0.0;
-}
-
-void OpenCLParallelCalcHarmonicAngleForceKernel::copyParametersToContext(ContextImpl& context, const HarmonicAngleForce& force, int firstAngle, int lastAngle) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).copyParametersToContext(context, force, firstAngle, lastAngle);
-}
-
-class OpenCLParallelCalcCustomAngleForceKernel::Task : public OpenCLContext::WorkTask {
-public:
-    Task(ContextImpl& context, CommonCalcCustomAngleForceKernel& kernel, bool includeForce,
-            bool includeEnergy, double& energy) : context(context), kernel(kernel),
-            includeForce(includeForce), includeEnergy(includeEnergy), energy(energy) {
-    }
-    void execute() {
-        energy += kernel.execute(context, includeForce, includeEnergy);
-    }
-private:
-    ContextImpl& context;
-    CommonCalcCustomAngleForceKernel& kernel;
-    bool includeForce, includeEnergy;
-    double& energy;
-};
-
-OpenCLParallelCalcCustomAngleForceKernel::OpenCLParallelCalcCustomAngleForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, const System& system) :
-        CalcCustomAngleForceKernel(name, platform), data(data) {
-    for (int i = 0; i < (int) data.contexts.size(); i++)
-        kernels.push_back(Kernel(new CommonCalcCustomAngleForceKernel(name, platform, *data.contexts[i], system)));
-}
-
-void OpenCLParallelCalcCustomAngleForceKernel::initialize(const System& system, const CustomAngleForce& force) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).initialize(system, force);
-}
-
-double OpenCLParallelCalcCustomAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    for (int i = 0; i < (int) data.contexts.size(); i++) {
-        OpenCLContext& cl = *data.contexts[i];
-        ComputeContext::WorkThread& thread = cl.getWorkThread();
-        thread.addTask(new Task(context, getKernel(i), includeForces, includeEnergy, data.contextEnergy[i]));
-    }
-    return 0.0;
-}
-
-void OpenCLParallelCalcCustomAngleForceKernel::copyParametersToContext(ContextImpl& context, const CustomAngleForce& force, int firstAngle, int lastAngle) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).copyParametersToContext(context, force, firstAngle, lastAngle);
-}
-
-class OpenCLParallelCalcPeriodicTorsionForceKernel::Task : public OpenCLContext::WorkTask {
-public:
-    Task(ContextImpl& context, CommonCalcPeriodicTorsionForceKernel& kernel, bool includeForce,
-            bool includeEnergy, double& energy) : context(context), kernel(kernel),
-            includeForce(includeForce), includeEnergy(includeEnergy), energy(energy) {
-    }
-    void execute() {
-        energy += kernel.execute(context, includeForce, includeEnergy);
-    }
-private:
-    ContextImpl& context;
-    CommonCalcPeriodicTorsionForceKernel& kernel;
-    bool includeForce, includeEnergy;
-    double& energy;
-};
-
-OpenCLParallelCalcPeriodicTorsionForceKernel::OpenCLParallelCalcPeriodicTorsionForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, const System& system) :
-        CalcPeriodicTorsionForceKernel(name, platform), data(data) {
-    for (int i = 0; i < (int) data.contexts.size(); i++)
-        kernels.push_back(Kernel(new CommonCalcPeriodicTorsionForceKernel(name, platform, *data.contexts[i], system)));
-}
-
-void OpenCLParallelCalcPeriodicTorsionForceKernel::initialize(const System& system, const PeriodicTorsionForce& force) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).initialize(system, force);
-}
-
-double OpenCLParallelCalcPeriodicTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    for (int i = 0; i < (int) data.contexts.size(); i++) {
-        OpenCLContext& cl = *data.contexts[i];
-        ComputeContext::WorkThread& thread = cl.getWorkThread();
-        thread.addTask(new Task(context, getKernel(i), includeForces, includeEnergy, data.contextEnergy[i]));
-    }
-    return 0.0;
-}
-
-void OpenCLParallelCalcPeriodicTorsionForceKernel::copyParametersToContext(ContextImpl& context, const PeriodicTorsionForce& force, int firstTorsion, int lastTorsion) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).copyParametersToContext(context, force, firstTorsion, lastTorsion);
-}
-
-class OpenCLParallelCalcRBTorsionForceKernel::Task : public OpenCLContext::WorkTask {
-public:
-    Task(ContextImpl& context, CommonCalcRBTorsionForceKernel& kernel, bool includeForce,
-            bool includeEnergy, double& energy) : context(context), kernel(kernel),
-            includeForce(includeForce), includeEnergy(includeEnergy), energy(energy) {
-    }
-    void execute() {
-        energy += kernel.execute(context, includeForce, includeEnergy);
-    }
-private:
-    ContextImpl& context;
-    CommonCalcRBTorsionForceKernel& kernel;
-    bool includeForce, includeEnergy;
-    double& energy;
-};
-
-OpenCLParallelCalcRBTorsionForceKernel::OpenCLParallelCalcRBTorsionForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, const System& system) :
-        CalcRBTorsionForceKernel(name, platform), data(data) {
-    for (int i = 0; i < (int) data.contexts.size(); i++)
-        kernels.push_back(Kernel(new CommonCalcRBTorsionForceKernel(name, platform, *data.contexts[i], system)));
-}
-
-void OpenCLParallelCalcRBTorsionForceKernel::initialize(const System& system, const RBTorsionForce& force) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).initialize(system, force);
-}
-
-double OpenCLParallelCalcRBTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    for (int i = 0; i < (int) data.contexts.size(); i++) {
-        OpenCLContext& cl = *data.contexts[i];
-        ComputeContext::WorkThread& thread = cl.getWorkThread();
-        thread.addTask(new Task(context, getKernel(i), includeForces, includeEnergy, data.contextEnergy[i]));
-    }
-    return 0.0;
-}
-
-void OpenCLParallelCalcRBTorsionForceKernel::copyParametersToContext(ContextImpl& context, const RBTorsionForce& force) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).copyParametersToContext(context, force);
-}
-
-class OpenCLParallelCalcCMAPTorsionForceKernel::Task : public OpenCLContext::WorkTask {
-public:
-    Task(ContextImpl& context, CommonCalcCMAPTorsionForceKernel& kernel, bool includeForce,
-            bool includeEnergy, double& energy) : context(context), kernel(kernel),
-            includeForce(includeForce), includeEnergy(includeEnergy), energy(energy) {
-    }
-    void execute() {
-        energy += kernel.execute(context, includeForce, includeEnergy);
-    }
-private:
-    ContextImpl& context;
-    CommonCalcCMAPTorsionForceKernel& kernel;
-    bool includeForce, includeEnergy;
-    double& energy;
-};
-
-OpenCLParallelCalcCMAPTorsionForceKernel::OpenCLParallelCalcCMAPTorsionForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, const System& system) :
-        CalcCMAPTorsionForceKernel(name, platform), data(data) {
-    for (int i = 0; i < (int) data.contexts.size(); i++)
-        kernels.push_back(Kernel(new CommonCalcCMAPTorsionForceKernel(name, platform, *data.contexts[i], system)));
-}
-
-void OpenCLParallelCalcCMAPTorsionForceKernel::initialize(const System& system, const CMAPTorsionForce& force) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).initialize(system, force);
-}
-
-double OpenCLParallelCalcCMAPTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    for (int i = 0; i < (int) data.contexts.size(); i++) {
-        OpenCLContext& cl = *data.contexts[i];
-        ComputeContext::WorkThread& thread = cl.getWorkThread();
-        thread.addTask(new Task(context, getKernel(i), includeForces, includeEnergy, data.contextEnergy[i]));
-    }
-    return 0.0;
-}
-
-void OpenCLParallelCalcCMAPTorsionForceKernel::copyParametersToContext(ContextImpl& context, const CMAPTorsionForce& force) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).copyParametersToContext(context, force);
-}
-
-class OpenCLParallelCalcCustomTorsionForceKernel::Task : public OpenCLContext::WorkTask {
-public:
-    Task(ContextImpl& context, CommonCalcCustomTorsionForceKernel& kernel, bool includeForce,
-            bool includeEnergy, double& energy) : context(context), kernel(kernel),
-            includeForce(includeForce), includeEnergy(includeEnergy), energy(energy) {
-    }
-    void execute() {
-        energy += kernel.execute(context, includeForce, includeEnergy);
-    }
-private:
-    ContextImpl& context;
-    CommonCalcCustomTorsionForceKernel& kernel;
-    bool includeForce, includeEnergy;
-    double& energy;
-};
-
-OpenCLParallelCalcCustomTorsionForceKernel::OpenCLParallelCalcCustomTorsionForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, const System& system) :
-        CalcCustomTorsionForceKernel(name, platform), data(data) {
-    for (int i = 0; i < (int) data.contexts.size(); i++)
-        kernels.push_back(Kernel(new CommonCalcCustomTorsionForceKernel(name, platform, *data.contexts[i], system)));
-}
-
-void OpenCLParallelCalcCustomTorsionForceKernel::initialize(const System& system, const CustomTorsionForce& force) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).initialize(system, force);
-}
-
-double OpenCLParallelCalcCustomTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    for (int i = 0; i < (int) data.contexts.size(); i++) {
-        OpenCLContext& cl = *data.contexts[i];
-        ComputeContext::WorkThread& thread = cl.getWorkThread();
-        thread.addTask(new Task(context, getKernel(i), includeForces, includeEnergy, data.contextEnergy[i]));
-    }
-    return 0.0;
-}
-
-void OpenCLParallelCalcCustomTorsionForceKernel::copyParametersToContext(ContextImpl& context, const CustomTorsionForce& force, int firstTorsion, int lastTorsion) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).copyParametersToContext(context, force, firstTorsion, lastTorsion);
-}
-
 class OpenCLParallelCalcNonbondedForceKernel::Task : public OpenCLContext::WorkTask {
 public:
     Task(ContextImpl& context, OpenCLCalcNonbondedForceKernel& kernel, bool includeForce,
@@ -587,167 +238,3 @@ void OpenCLParallelCalcNonbondedForceKernel::getPMEParameters(double& alpha, int
 void OpenCLParallelCalcNonbondedForceKernel::getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
     dynamic_cast<const OpenCLCalcNonbondedForceKernel&>(kernels[0].getImpl()).getLJPMEParameters(alpha, nx, ny, nz);
 }
-
-class OpenCLParallelCalcCustomNonbondedForceKernel::Task : public OpenCLContext::WorkTask {
-public:
-    Task(ContextImpl& context, CommonCalcCustomNonbondedForceKernel& kernel, bool includeForce,
-            bool includeEnergy, double& energy) : context(context), kernel(kernel),
-            includeForce(includeForce), includeEnergy(includeEnergy), energy(energy) {
-    }
-    void execute() {
-        energy += kernel.execute(context, includeForce, includeEnergy);
-    }
-private:
-    ContextImpl& context;
-    CommonCalcCustomNonbondedForceKernel& kernel;
-    bool includeForce, includeEnergy;
-    double& energy;
-};
-
-OpenCLParallelCalcCustomNonbondedForceKernel::OpenCLParallelCalcCustomNonbondedForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, const System& system) :
-        CalcCustomNonbondedForceKernel(name, platform), data(data) {
-    for (int i = 0; i < (int) data.contexts.size(); i++)
-        kernels.push_back(Kernel(new CommonCalcCustomNonbondedForceKernel(name, platform, *data.contexts[i], system)));
-}
-
-void OpenCLParallelCalcCustomNonbondedForceKernel::initialize(const System& system, const CustomNonbondedForce& force) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).initialize(system, force);
-}
-
-double OpenCLParallelCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    for (int i = 0; i < (int) data.contexts.size(); i++) {
-        OpenCLContext& cl = *data.contexts[i];
-        ComputeContext::WorkThread& thread = cl.getWorkThread();
-        thread.addTask(new Task(context, getKernel(i), includeForces, includeEnergy, data.contextEnergy[i]));
-    }
-    return 0.0;
-}
-
-void OpenCLParallelCalcCustomNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force, int firstParticle, int lastParticle) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).copyParametersToContext(context, force, firstParticle, lastParticle);
-}
-
-class OpenCLParallelCalcCustomExternalForceKernel::Task : public OpenCLContext::WorkTask {
-public:
-    Task(ContextImpl& context, CommonCalcCustomExternalForceKernel& kernel, bool includeForce,
-            bool includeEnergy, double& energy) : context(context), kernel(kernel),
-            includeForce(includeForce), includeEnergy(includeEnergy), energy(energy) {
-    }
-    void execute() {
-        energy += kernel.execute(context, includeForce, includeEnergy);
-    }
-private:
-    ContextImpl& context;
-    CommonCalcCustomExternalForceKernel& kernel;
-    bool includeForce, includeEnergy;
-    double& energy;
-};
-
-OpenCLParallelCalcCustomExternalForceKernel::OpenCLParallelCalcCustomExternalForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, const System& system) :
-        CalcCustomExternalForceKernel(name, platform), data(data) {
-    for (int i = 0; i < (int) data.contexts.size(); i++)
-        kernels.push_back(Kernel(new CommonCalcCustomExternalForceKernel(name, platform, *data.contexts[i], system)));
-}
-
-void OpenCLParallelCalcCustomExternalForceKernel::initialize(const System& system, const CustomExternalForce& force) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).initialize(system, force);
-}
-
-double OpenCLParallelCalcCustomExternalForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    for (int i = 0; i < (int) data.contexts.size(); i++) {
-        OpenCLContext& cl = *data.contexts[i];
-        ComputeContext::WorkThread& thread = cl.getWorkThread();
-        thread.addTask(new Task(context, getKernel(i), includeForces, includeEnergy, data.contextEnergy[i]));
-    }
-    return 0.0;
-}
-
-void OpenCLParallelCalcCustomExternalForceKernel::copyParametersToContext(ContextImpl& context, const CustomExternalForce& force, int firstParticle, int lastParticle) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).copyParametersToContext(context, force, firstParticle, lastParticle);
-}
-
-class OpenCLParallelCalcCustomHbondForceKernel::Task : public OpenCLContext::WorkTask {
-public:
-    Task(ContextImpl& context, CommonCalcCustomHbondForceKernel& kernel, bool includeForce,
-            bool includeEnergy, double& energy) : context(context), kernel(kernel),
-            includeForce(includeForce), includeEnergy(includeEnergy), energy(energy) {
-    }
-    void execute() {
-        energy += kernel.execute(context, includeForce, includeEnergy);
-    }
-private:
-    ContextImpl& context;
-    CommonCalcCustomHbondForceKernel& kernel;
-    bool includeForce, includeEnergy;
-    double& energy;
-};
-
-OpenCLParallelCalcCustomHbondForceKernel::OpenCLParallelCalcCustomHbondForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, const System& system) :
-        CalcCustomHbondForceKernel(name, platform), data(data) {
-    for (int i = 0; i < (int) data.contexts.size(); i++)
-        kernels.push_back(Kernel(new CommonCalcCustomHbondForceKernel(name, platform, *data.contexts[i], system)));
-}
-
-void OpenCLParallelCalcCustomHbondForceKernel::initialize(const System& system, const CustomHbondForce& force) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).initialize(system, force);
-}
-
-double OpenCLParallelCalcCustomHbondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    for (int i = 0; i < (int) data.contexts.size(); i++) {
-        OpenCLContext& cl = *data.contexts[i];
-        ComputeContext::WorkThread& thread = cl.getWorkThread();
-        thread.addTask(new Task(context, getKernel(i), includeForces, includeEnergy, data.contextEnergy[i]));
-    }
-    return 0.0;
-}
-
-void OpenCLParallelCalcCustomHbondForceKernel::copyParametersToContext(ContextImpl& context, const CustomHbondForce& force) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).copyParametersToContext(context, force);
-}
-
-class OpenCLParallelCalcCustomCompoundBondForceKernel::Task : public OpenCLContext::WorkTask {
-public:
-    Task(ContextImpl& context, CommonCalcCustomCompoundBondForceKernel& kernel, bool includeForce,
-            bool includeEnergy, double& energy) : context(context), kernel(kernel),
-            includeForce(includeForce), includeEnergy(includeEnergy), energy(energy) {
-    }
-    void execute() {
-        energy += kernel.execute(context, includeForce, includeEnergy);
-    }
-private:
-    ContextImpl& context;
-    CommonCalcCustomCompoundBondForceKernel& kernel;
-    bool includeForce, includeEnergy;
-    double& energy;
-};
-
-OpenCLParallelCalcCustomCompoundBondForceKernel::OpenCLParallelCalcCustomCompoundBondForceKernel(std::string name, const Platform& platform, OpenCLPlatform::PlatformData& data, const System& system) :
-        CalcCustomCompoundBondForceKernel(name, platform), data(data) {
-    for (int i = 0; i < (int) data.contexts.size(); i++)
-        kernels.push_back(Kernel(new CommonCalcCustomCompoundBondForceKernel(name, platform, *data.contexts[i], system)));
-}
-
-void OpenCLParallelCalcCustomCompoundBondForceKernel::initialize(const System& system, const CustomCompoundBondForce& force) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).initialize(system, force);
-}
-
-double OpenCLParallelCalcCustomCompoundBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    for (int i = 0; i < (int) data.contexts.size(); i++) {
-        OpenCLContext& cl = *data.contexts[i];
-        ComputeContext::WorkThread& thread = cl.getWorkThread();
-        thread.addTask(new Task(context, getKernel(i), includeForces, includeEnergy, data.contextEnergy[i]));
-    }
-    return 0.0;
-}
-
-void OpenCLParallelCalcCustomCompoundBondForceKernel::copyParametersToContext(ContextImpl& context, const CustomCompoundBondForce& force) {
-    for (int i = 0; i < (int) kernels.size(); i++)
-        getKernel(i).copyParametersToContext(context, force);
-}

platforms/opencl/tests/TestOpenCLCustomAngleForce.cpp

@@ -32,55 +32,6 @@
 #include "OpenCLTests.h"
 #include "TestCustomAngleForce.h"
 
-void testParallelComputation() {
-    System system;
-    const int numParticles = 200;
-    for (int i = 0; i < numParticles; i++)
-        system.addParticle(1.0);
-    CustomAngleForce* force = new CustomAngleForce("k*(theta-theta0)^2");
-    force->addPerAngleParameter("k");
-    force->addPerAngleParameter("theta0");
-    for (int i = 2; i < numParticles; i++)
-        force->addAngle(i-2, i-1, i, {1.0, 1.1});
-    system.addForce(force);
-    vector<Vec3> positions(numParticles);
-    for (int i = 0; i < numParticles; i++)
-        positions[i] = Vec3(i, i%2, 0);
-    VerletIntegrator integrator1(0.01);
-    Context context1(system, integrator1, platform);
-    context1.setPositions(positions);
-    State state1 = context1.getState(State::Forces | State::Energy);
-    VerletIntegrator integrator2(0.01);
-
-    map<string, string> props;
-    string deviceIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLDeviceIndex());
-    props[OpenCLPlatform::OpenCLDeviceIndex()] = deviceIndex+","+deviceIndex;
-    string platformIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLPlatformIndex());
-    props[OpenCLPlatform::OpenCLPlatformIndex()] = platformIndex;
-
-    Context context2(system, integrator2, platform, props);
-    context2.setPositions(positions);
-    State state2 = context2.getState(State::Forces | State::Energy);
-    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
-    for (int i = 0; i < numParticles; i++)
-        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
-
-    // Try updating some parameters and see if they still match.
-
-    vector<double> params;
-    for (int i = 95; i < 102; i++) {
-        int p1, p2, p3;
-        force->getAngleParameters(i, p1, p2, p3, params);
-        force->setAngleParameters(i, p1, p2, p3, {2.0, 1.2});
-    }
-    force->updateParametersInContext(context1);
-    force->updateParametersInContext(context2);
-    State state3 = context1.getState(State::Energy);
-    State state4 = context2.getState(State::Energy);
-    ASSERT_EQUAL_TOL(state3.getPotentialEnergy(), state4.getPotentialEnergy(), 1e-5);
-    ASSERT(fabs(state1.getPotentialEnergy()-state3.getPotentialEnergy()) > 0.1);
-}
-
 void runPlatformTests() {
     testParallelComputation();
 }

platforms/opencl/tests/TestOpenCLCustomBondForce.cpp

@@ -32,55 +32,6 @@
 #include "OpenCLTests.h"
 #include "TestCustomBondForce.h"
 
-void testParallelComputation() {
-    System system;
-    const int numParticles = 200;
-    for (int i = 0; i < numParticles; i++)
-        system.addParticle(1.0);
-    CustomBondForce* force = new CustomBondForce(("k*(r-r0)^2"));
-    force->addPerBondParameter("k");
-    force->addPerBondParameter("r0");
-    for (int i = 1; i < numParticles; i++)
-        force->addBond(i-1, i, {1.0, 1.1});
-    system.addForce(force);
-    vector<Vec3> positions(numParticles);
-    for (int i = 0; i < numParticles; i++)
-        positions[i] = Vec3(i, 0, 0);
-    VerletIntegrator integrator1(0.01);
-    Context context1(system, integrator1, platform);
-    context1.setPositions(positions);
-    State state1 = context1.getState(State::Forces | State::Energy);
-    VerletIntegrator integrator2(0.01);
-
-    map<string, string> props;
-    string deviceIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLDeviceIndex());
-    props[OpenCLPlatform::OpenCLDeviceIndex()] = deviceIndex+","+deviceIndex;
-    string platformIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLPlatformIndex());
-    props[OpenCLPlatform::OpenCLPlatformIndex()] = platformIndex;
-
-    Context context2(system, integrator2, platform, props);
-    context2.setPositions(positions);
-    State state2 = context2.getState(State::Forces | State::Energy);
-    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
-    for (int i = 0; i < numParticles; i++)
-        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
-
-    // Try updating some parameters and see if they still match.
-
-    vector<double> params;
-    for (int i = 95; i < 102; i++) {
-        int p1, p2;
-        force->getBondParameters(i, p1, p2, params);
-        force->setBondParameters(i, p1, p2, {2.0, 1.2});
-    }
-    force->updateParametersInContext(context1);
-    force->updateParametersInContext(context2);
-    State state3 = context1.getState(State::Energy);
-    State state4 = context2.getState(State::Energy);
-    ASSERT_EQUAL_TOL(state3.getPotentialEnergy(), state4.getPotentialEnergy(), 1e-5);
-    ASSERT(fabs(state1.getPotentialEnergy()-state3.getPotentialEnergy()) > 0.1);
-}
-
 void runPlatformTests() {
     testParallelComputation();
 }

platforms/opencl/tests/TestOpenCLCustomCompoundBondForce.cpp

@@ -32,43 +32,6 @@
 #include "OpenCLTests.h"
 #include "TestCustomCompoundBondForce.h"
 
-void testParallelComputation() {
-    System system;
-    const int numParticles = 200;
-    for (int i = 0; i < numParticles; i++)
-        system.addParticle(1.0);
-    CustomCompoundBondForce* force = new CustomCompoundBondForce(2, ("(distance(p1,p2)-1.1)^2"));
-    vector<int> particles(2);
-    vector<double> params;
-    for (int i = 1; i < numParticles; i++) {
-        particles[0] = i-1;
-        particles[1] = i;
-        force->addBond(particles, params);
-    }
-    system.addForce(force);
-    vector<Vec3> positions(numParticles);
-    for (int i = 0; i < numParticles; i++)
-        positions[i] = Vec3(i, 0, 0);
-    VerletIntegrator integrator1(0.01);
-    Context context1(system, integrator1, platform);
-    context1.setPositions(positions);
-    State state1 = context1.getState(State::Forces | State::Energy);
-    VerletIntegrator integrator2(0.01);
-
-    map<string, string> props;
-    string deviceIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLDeviceIndex());
-    props[OpenCLPlatform::OpenCLDeviceIndex()] = deviceIndex+","+deviceIndex;
-    string platformIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLPlatformIndex());
-    props[OpenCLPlatform::OpenCLPlatformIndex()] = platformIndex;
-
-    Context context2(system, integrator2, platform, props);
-    context2.setPositions(positions);
-    State state2 = context2.getState(State::Forces | State::Energy);
-    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
-    for (int i = 0; i < numParticles; i++)
-        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
-}
-
 void runPlatformTests() {
     testParallelComputation();
 }

platforms/opencl/tests/TestOpenCLCustomExternalForce.cpp

@@ -33,41 +33,6 @@
 #include "TestCustomExternalForce.h"
 #include "sfmt/SFMT.h"
 
-void testParallelComputation() {
-    System system;
-    const int numParticles = 200;
-    for (int i = 0; i < numParticles; i++)
-        system.addParticle(1.0);
-    CustomExternalForce* force = new CustomExternalForce("x^2+y^2+z^2");
-    vector<double> params;
-    for (int i = 0; i < numParticles; i++)
-        force->addParticle(i, params);
-    system.addForce(force);
-    OpenMM_SFMT::SFMT sfmt;
-    init_gen_rand(0, sfmt);
-    vector<Vec3> positions(numParticles);
-    for (int i = 0; i < numParticles; i++)
-        positions[i] = Vec3(5*genrand_real2(sfmt), 5*genrand_real2(sfmt), 5*genrand_real2(sfmt));
-    VerletIntegrator integrator1(0.01);
-    Context context1(system, integrator1, platform);
-    context1.setPositions(positions);
-    State state1 = context1.getState(State::Forces | State::Energy);
-    VerletIntegrator integrator2(0.01);
-
-    map<string, string> props;
-    string deviceIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLDeviceIndex());
-    props[OpenCLPlatform::OpenCLDeviceIndex()] = deviceIndex+","+deviceIndex;
-    string platformIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLPlatformIndex());
-    props[OpenCLPlatform::OpenCLPlatformIndex()] = platformIndex;
-
-    Context context2(system, integrator2, platform, props);
-    context2.setPositions(positions);
-    State state2 = context2.getState(State::Forces | State::Energy);
-    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
-    for (int i = 0; i < numParticles; i++)
-        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
-}
-
 void runPlatformTests() {
     testParallelComputation();
 }

platforms/opencl/tests/TestOpenCLCustomNonbondedForce.cpp

@@ -32,60 +32,6 @@
 #include "OpenCLTests.h"
 #include "TestCustomNonbondedForce.h"
 
-void testParallelComputation() {
-    System system;
-    const int numParticles = 200;
-    for (int i = 0; i < numParticles; i++)
-        system.addParticle(1.0);
-    CustomNonbondedForce* force = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6); sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)");
-    force->addPerParticleParameter("sigma");
-    force->addPerParticleParameter("eps");
-    for (int i = 0; i < numParticles; i++)
-        force->addParticle({0.5, 1.0});
-    system.addForce(force);
-    OpenMM_SFMT::SFMT sfmt;
-    init_gen_rand(0, sfmt);
-    vector<Vec3> positions(numParticles);
-    for (int i = 0; i < numParticles; i++)
-        positions[i] = Vec3(5*genrand_real2(sfmt), 5*genrand_real2(sfmt), 5*genrand_real2(sfmt));
-    for (int i = 0; i < numParticles; ++i)
-        for (int j = 0; j < i; ++j) {
-            Vec3 delta = positions[i]-positions[j];
-            if (delta.dot(delta) < 0.1)
-                force->addExclusion(i, j);
-        }
-    VerletIntegrator integrator1(0.01);
-    Context context1(system, integrator1, platform);
-    context1.setPositions(positions);
-    State state1 = context1.getState(State::Forces | State::Energy);
-    VerletIntegrator integrator2(0.01);
-
-    map<string, string> props;
-    string deviceIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLDeviceIndex());
-    props[OpenCLPlatform::OpenCLDeviceIndex()] = deviceIndex+","+deviceIndex;
-    string platformIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLPlatformIndex());
-    props[OpenCLPlatform::OpenCLPlatformIndex()] = platformIndex;
-
-    Context context2(system, integrator2, platform, props);
-    context2.setPositions(positions);
-    State state2 = context2.getState(State::Forces | State::Energy);
-    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
-    for (int i = 0; i < numParticles; i++)
-        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
-
-    // Try updating some parameters and see if they still match.
-
-    vector<double> params;
-    for (int i = 95; i < 102; i++)
-        force->setParticleParameters(i, {0.6, 2.0});
-    force->updateParametersInContext(context1);
-    force->updateParametersInContext(context2);
-    State state3 = context1.getState(State::Energy);
-    State state4 = context2.getState(State::Energy);
-    ASSERT_EQUAL_TOL(state3.getPotentialEnergy(), state4.getPotentialEnergy(), 1e-5);
-    ASSERT(fabs(state1.getPotentialEnergy()-state3.getPotentialEnergy()) > 0.1);
-}
-
 void runPlatformTests() {
     testParallelComputation();
 }

platforms/opencl/tests/TestOpenCLCustomTorsionForce.cpp

@@ -32,55 +32,6 @@
 #include "OpenCLTests.h"
 #include "TestCustomTorsionForce.h"
 
-void testParallelComputation() {
-    System system;
-    const int numParticles = 200;
-    for (int i = 0; i < numParticles; i++)
-        system.addParticle(1.0);
-    CustomTorsionForce* force = new CustomTorsionForce("k*sin(theta-theta0)");
-    force->addPerTorsionParameter("k");
-    force->addPerTorsionParameter("theta0");
-    for (int i = 3; i < numParticles; i++)
-        force->addTorsion(i-3, i-2, i-1, i, {1.0, 1.1});
-    system.addForce(force);
-    vector<Vec3> positions(numParticles);
-    for (int i = 0; i < numParticles; i++)
-        positions[i] = Vec3(i, i%2, i%3);
-    VerletIntegrator integrator1(0.01);
-    Context context1(system, integrator1, platform);
-    context1.setPositions(positions);
-    State state1 = context1.getState(State::Forces | State::Energy);
-    VerletIntegrator integrator2(0.01);
-
-    map<string, string> props;
-    string deviceIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLDeviceIndex());
-    props[OpenCLPlatform::OpenCLDeviceIndex()] = deviceIndex+","+deviceIndex;
-    string platformIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLPlatformIndex());
-    props[OpenCLPlatform::OpenCLPlatformIndex()] = platformIndex;
-
-    Context context2(system, integrator2, platform, props);
-    context2.setPositions(positions);
-    State state2 = context2.getState(State::Forces | State::Energy);
-    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
-    for (int i = 0; i < numParticles; i++)
-        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
-
-    // Try updating some parameters and see if they still match.
-
-    vector<double> params;
-    for (int i = 95; i < 102; i++) {
-        int p1, p2, p3, p4;
-        force->getTorsionParameters(i, p1, p2, p3, p4, params);
-        force->setTorsionParameters(i, p1, p2, p3, p4, {2.0, 1.2});
-    }
-    force->updateParametersInContext(context1);
-    force->updateParametersInContext(context2);
-    State state3 = context1.getState(State::Energy);
-    State state4 = context2.getState(State::Energy);
-    ASSERT_EQUAL_TOL(state3.getPotentialEnergy(), state4.getPotentialEnergy(), 1e-5);
-    ASSERT(fabs(state1.getPotentialEnergy()-state3.getPotentialEnergy()) > 0.1);
-}
-
 void runPlatformTests() {
     testParallelComputation();
 }

platforms/opencl/tests/TestOpenCLHarmonicAngleForce.cpp

@@ -32,54 +32,6 @@
 
 #include "OpenCLTests.h"
 #include "TestHarmonicAngleForce.h"
-#include <map>
-
-void testParallelComputation() {
-    System system;
-    const int numParticles = 200;
-    for (int i = 0; i < numParticles; i++)
-        system.addParticle(1.0);
-    HarmonicAngleForce* force = new HarmonicAngleForce();
-    for (int i = 2; i < numParticles; i++)
-        force->addAngle(i-2, i-1, i, 1.1, i);
-    system.addForce(force);
-    vector<Vec3> positions(numParticles);
-    for (int i = 0; i < numParticles; i++)
-        positions[i] = Vec3(i, i%2, 0);
-    VerletIntegrator integrator1(0.01);
-    Context context1(system, integrator1, platform);
-    context1.setPositions(positions);
-    State state1 = context1.getState(State::Forces | State::Energy);
-    VerletIntegrator integrator2(0.01);
-
-    map<string, string> props;
-    string deviceIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLDeviceIndex());
-    props[OpenCLPlatform::OpenCLDeviceIndex()] = deviceIndex+","+deviceIndex;
-    string platformIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLPlatformIndex());
-    props[OpenCLPlatform::OpenCLPlatformIndex()] = platformIndex;
-
-    Context context2(system, integrator2, platform, props);
-    context2.setPositions(positions);
-    State state2 = context2.getState(State::Forces | State::Energy);
-    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
-    for (int i = 0; i < numParticles; i++)
-        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
-
-    // Try updating some parameters and see if they still match.
-
-    for (int i = 95; i < 102; i++) {
-        int p1, p2, p3;
-        double angle, k;
-        force->getAngleParameters(i, p1, p2, p3, angle, k);
-        force->setAngleParameters(i, p1, p2, p3, angle+0.1, 2*k);
-    }
-    force->updateParametersInContext(context1);
-    force->updateParametersInContext(context2);
-    State state3 = context1.getState(State::Energy);
-    State state4 = context2.getState(State::Energy);
-    ASSERT_EQUAL_TOL(state3.getPotentialEnergy(), state4.getPotentialEnergy(), 1e-5);
-    ASSERT(fabs(state1.getPotentialEnergy()-state3.getPotentialEnergy()) > 0.1);
-}
 
 void runPlatformTests() {
     testParallelComputation();

platforms/opencl/tests/TestOpenCLHarmonicBondForce.cpp

@@ -31,54 +31,6 @@
 
 #include "OpenCLTests.h"
 #include "TestHarmonicBondForce.h"
-#include <map>
-
-void testParallelComputation() {
-    System system;
-    const int numParticles = 200;
-    for (int i = 0; i < numParticles; i++)
-        system.addParticle(1.0);
-    HarmonicBondForce* force = new HarmonicBondForce();
-    for (int i = 1; i < numParticles; i++)
-        force->addBond(i-1, i, 1.1, i);
-    system.addForce(force);
-    vector<Vec3> positions(numParticles);
-    for (int i = 0; i < numParticles; i++)
-        positions[i] = Vec3(i, 0, 0);
-    VerletIntegrator integrator1(0.01);
-    Context context1(system, integrator1, platform);
-    context1.setPositions(positions);
-    State state1 = context1.getState(State::Forces | State::Energy);
-    VerletIntegrator integrator2(0.01);
-
-    map<string, string> props;
-    string deviceIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLDeviceIndex());
-    props[OpenCLPlatform::OpenCLDeviceIndex()] = deviceIndex+","+deviceIndex;
-    string platformIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLPlatformIndex());
-    props[OpenCLPlatform::OpenCLPlatformIndex()] = platformIndex;
-    
-    Context context2(system, integrator2, platform, props);
-    context2.setPositions(positions);
-    State state2 = context2.getState(State::Forces | State::Energy);
-    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
-    for (int i = 0; i < numParticles; i++)
-        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
-
-    // Try updating some parameters and see if they still match.
-
-    for (int i = 95; i < 102; i++) {
-        int p1, p2;
-        double length, k;
-        force->getBondParameters(i, p1, p2, length, k);
-        force->setBondParameters(i, p1, p2, length+0.1, 2*k);
-    }
-    force->updateParametersInContext(context1);
-    force->updateParametersInContext(context2);
-    State state3 = context1.getState(State::Energy);
-    State state4 = context2.getState(State::Energy);
-    ASSERT_EQUAL_TOL(state3.getPotentialEnergy(), state4.getPotentialEnergy(), 1e-5);
-    ASSERT(fabs(state1.getPotentialEnergy()-state3.getPotentialEnergy()) > 0.1);
-}
 
 void runPlatformTests() {
     testParallelComputation();

platforms/opencl/tests/TestOpenCLNonbondedForce.cpp

@@ -35,110 +35,6 @@
 #include "OpenCLTests.h"
 #include "TestNonbondedForce.h"
 #include "opencl.hpp"
-#include <string>
-
-void testParallelComputation(NonbondedForce::NonbondedMethod method) {
-    System system;
-    const int numParticles = 200;
-    for (int i = 0; i < numParticles; i++)
-        system.addParticle(1.0);
-    NonbondedForce* force = new NonbondedForce();
-    for (int i = 0; i < numParticles; i++)
-        force->addParticle(i%2-0.5, 0.5, 1.0);
-    force->setNonbondedMethod(method);
-    system.addForce(force);
-    system.setDefaultPeriodicBoxVectors(Vec3(5,0,0), Vec3(0,5,0), Vec3(0,0,5));
-    OpenMM_SFMT::SFMT sfmt;
-    init_gen_rand(0, sfmt);
-    vector<Vec3> positions(numParticles);
-    for (int i = 0; i < numParticles; i++)
-        positions[i] = Vec3(5*genrand_real2(sfmt), 5*genrand_real2(sfmt), 5*genrand_real2(sfmt));
-    force->addGlobalParameter("scale", 0.5);
-    for (int i = 0; i < numParticles; ++i)
-        for (int j = 0; j < i; ++j) {
-            Vec3 delta = positions[i]-positions[j];
-            if (delta.dot(delta) < 0.1) {
-                force->addException(i, j, 0, 1, 0);
-            }
-            else if (delta.dot(delta) < 0.2) {
-                int index = force->addException(i, j, 0.5, 1, 1.0);
-                force->addExceptionParameterOffset("scale", index, 0.5, 0.4, 0.3);
-            }
-        }
-    
-    // Create two contexts, one with a single device and one with two devices.
-    
-    VerletIntegrator integrator1(0.01);
-    Context context1(system, integrator1, platform);
-    context1.setPositions(positions);
-    State state1 = context1.getState(State::Forces | State::Energy);
-    VerletIntegrator integrator2(0.01);
-
-    map<string, string> props;
-    string deviceIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLDeviceIndex());
-    props[OpenCLPlatform::OpenCLDeviceIndex()] = deviceIndex+","+deviceIndex;
-    string platformIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLPlatformIndex());
-    props[OpenCLPlatform::OpenCLPlatformIndex()] = platformIndex;
-
-    Context context2(system, integrator2, platform, props);
-    context2.setPositions(positions);
-    State state2 = context2.getState(State::Forces | State::Energy);
-    
-    // See if they agree.
-    
-    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
-    for (int i = 0; i < numParticles; i++)
-        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
-    
-    // Modify some parameters and see if they still agree.
-
-    for (int i = 0; i < numParticles; i += 5) {
-        double charge, sigma, epsilon;
-        force->getParticleParameters(i, charge, sigma, epsilon);
-        force->setParticleParameters(i, 0.9*charge, sigma, epsilon);
-    }
-    for (int i = force->getNumExceptions()/2-10; i < force->getNumExceptions()/2+10; i++) {
-        int p1, p2;
-        double charge, sigma, epsilon;
-        force->getExceptionParameters(i, p1, p2, charge, sigma, epsilon);
-        if (epsilon != 0)
-            force->setExceptionParameters(i, p1, p2, charge, sigma, 2*epsilon);
-    }
-    force->updateParametersInContext(context1);
-    force->updateParametersInContext(context2);
-    state1 = context1.getState(State::Forces | State::Energy);
-    state2 = context2.getState(State::Forces | State::Energy);
-    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
-    for (int i = 0; i < numParticles; i++)
-        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
-}
-
-void testReordering() {
-    // Check that reordering of atoms doesn't alter their positions.
-    
-    const int numParticles = 200;
-    System system;
-    system.setDefaultPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(2.1, 6, 0), Vec3(-1.5, -0.5, 6));
-    NonbondedForce *nonbonded = new NonbondedForce();
-    nonbonded->setNonbondedMethod(NonbondedForce::PME);
-    system.addForce(nonbonded);
-    vector<Vec3> positions;
-    OpenMM_SFMT::SFMT sfmt;
-    init_gen_rand(0, sfmt);
-    for (int i = 0; i < numParticles; i++) {
-        system.addParticle(1.0);
-        nonbonded->addParticle(0.0, 0.0, 0.0);
-        positions.push_back(Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5)*20);
-    }
-    VerletIntegrator integrator(0.001);
-    Context context(system, integrator, platform);
-    context.setPositions(positions);
-    integrator.step(1);
-    State state = context.getState(State::Positions | State::Velocities);
-    for (int i = 0; i < numParticles; i++) {
-        ASSERT_EQUAL_VEC(positions[i], state.getPositions()[i], 1e-6);
-    }
-}
 
 bool canRunHugeTest() {
     // Create a minimal context just to see which platform and device are being used.

platforms/opencl/tests/TestOpenCLPeriodicTorsionForce.cpp

@@ -32,53 +32,6 @@
 #include "OpenCLTests.h"
 #include "TestPeriodicTorsionForce.h"
 
-void testParallelComputation() {
-    System system;
-    const int numParticles = 200;
-    for (int i = 0; i < numParticles; i++)
-        system.addParticle(1.0);
-    PeriodicTorsionForce* force = new PeriodicTorsionForce();
-    for (int i = 3; i < numParticles; i++)
-        force->addTorsion(i-3, i-2, i-1, i, 2, 1.1, i);
-    system.addForce(force);
-    vector<Vec3> positions(numParticles);
-    for (int i = 0; i < numParticles; i++)
-        positions[i] = Vec3(i, i%2, i%3);
-    VerletIntegrator integrator1(0.01);
-    Context context1(system, integrator1, platform);
-    context1.setPositions(positions);
-    State state1 = context1.getState(State::Forces | State::Energy);
-    VerletIntegrator integrator2(0.01);
-
-    map<string, string> props;
-    string deviceIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLDeviceIndex());
-    props[OpenCLPlatform::OpenCLDeviceIndex()] = deviceIndex+","+deviceIndex;
-    string platformIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLPlatformIndex());
-    props[OpenCLPlatform::OpenCLPlatformIndex()] = platformIndex;
-
-    Context context2(system, integrator2, platform, props);
-    context2.setPositions(positions);
-    State state2 = context2.getState(State::Forces | State::Energy);
-    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
-    for (int i = 0; i < numParticles; i++)
-        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
-
-    // Try updating some parameters and see if they still match.
-
-    for (int i = 95; i < 102; i++) {
-        int p1, p2, p3, p4, periodicity;
-        double phase, k;
-        force->getTorsionParameters(i, p1, p2, p3, p4, periodicity, phase, k);
-        force->setTorsionParameters(i, p1, p2, p3, p4, periodicity, phase+0.1, 2*k);
-    }
-    force->updateParametersInContext(context1);
-    force->updateParametersInContext(context2);
-    State state3 = context1.getState(State::Energy);
-    State state4 = context2.getState(State::Energy);
-    ASSERT_EQUAL_TOL(state3.getPotentialEnergy(), state4.getPotentialEnergy(), 1e-5);
-    ASSERT(fabs(state1.getPotentialEnergy()-state3.getPotentialEnergy()) > 0.1);
-}
-
 void runPlatformTests() {
     testParallelComputation();
 }

platforms/opencl/tests/TestOpenCLRBTorsionForce.cpp

@@ -32,38 +32,6 @@
 #include "OpenCLTests.h"
 #include "TestRBTorsionForce.h"
 
-void testParallelComputation() {
-    System system;
-    const int numParticles = 200;
-    for (int i = 0; i < numParticles; i++)
-        system.addParticle(1.0);
-    RBTorsionForce* force = new RBTorsionForce();
-    for (int i = 3; i < numParticles; i++)
-        force->addTorsion(i-3, i-2, i-1, i, 2, 0.1*i, 0.5*i, i, 1, 1);
-    system.addForce(force);
-    vector<Vec3> positions(numParticles);
-    for (int i = 0; i < numParticles; i++)
-        positions[i] = Vec3(i, i%2, i%3);
-    VerletIntegrator integrator1(0.01);
-    Context context1(system, integrator1, platform);
-    context1.setPositions(positions);
-    State state1 = context1.getState(State::Forces | State::Energy);
-    VerletIntegrator integrator2(0.01);
-
-    map<string, string> props;
-    string deviceIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLDeviceIndex());
-    props[OpenCLPlatform::OpenCLDeviceIndex()] = deviceIndex+","+deviceIndex;
-    string platformIndex = platform.getPropertyValue(context1, OpenCLPlatform::OpenCLPlatformIndex());
-    props[OpenCLPlatform::OpenCLPlatformIndex()] = platformIndex;
-
-    Context context2(system, integrator2, platform, props);
-    context2.setPositions(positions);
-    State state2 = context2.getState(State::Forces | State::Energy);
-    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
-    for (int i = 0; i < numParticles; i++)
-        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
-}
-
 void runPlatformTests() {
     testParallelComputation();
 }

tests/TestCustomAngleForce.h

@@ -216,6 +216,51 @@ void testEnergyParameterDerivatives() {
     }
 }
 
+void testParallelComputation() {
+    System system;
+    const int numParticles = 200;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    CustomAngleForce* force = new CustomAngleForce("k*(theta-theta0)^2");
+    force->addPerAngleParameter("k");
+    force->addPerAngleParameter("theta0");
+    for (int i = 2; i < numParticles; i++)
+        force->addAngle(i-2, i-1, i, {1.0, 1.1});
+    system.addForce(force);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numParticles; i++)
+        positions[i] = Vec3(i, i%2, 0);
+    VerletIntegrator integrator1(0.01);
+    Context context1(system, integrator1, platform);
+    context1.setPositions(positions);
+    State state1 = context1.getState(State::Forces | State::Energy);
+    VerletIntegrator integrator2(0.01);
+    string deviceIndex = platform.getPropertyValue(context1, "DeviceIndex");
+    map<string, string> props;
+    props["DeviceIndex"] = deviceIndex+","+deviceIndex;
+    Context context2(system, integrator2, platform, props);
+    context2.setPositions(positions);
+    State state2 = context2.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
+
+    // Try updating some parameters and see if they still match.
+
+    vector<double> params;
+    for (int i = 95; i < 102; i++) {
+        int p1, p2, p3;
+        force->getAngleParameters(i, p1, p2, p3, params);
+        force->setAngleParameters(i, p1, p2, p3, {2.0, 1.2});
+    }
+    force->updateParametersInContext(context1);
+    force->updateParametersInContext(context2);
+    State state3 = context1.getState(State::Energy);
+    State state4 = context2.getState(State::Energy);
+    ASSERT_EQUAL_TOL(state3.getPotentialEnergy(), state4.getPotentialEnergy(), 1e-5);
+    ASSERT(fabs(state1.getPotentialEnergy()-state3.getPotentialEnergy()) > 0.1);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {

tests/TestCustomBondForce.h

@@ -235,6 +235,51 @@ void testEnergyParameterDerivatives() {
     }
 }
 
+void testParallelComputation() {
+    System system;
+    const int numParticles = 200;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    CustomBondForce* force = new CustomBondForce(("k*(r-r0)^2"));
+    force->addPerBondParameter("k");
+    force->addPerBondParameter("r0");
+    for (int i = 1; i < numParticles; i++)
+        force->addBond(i-1, i, {1.0, 1.1});
+    system.addForce(force);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numParticles; i++)
+        positions[i] = Vec3(i, 0, 0);
+    VerletIntegrator integrator1(0.01);
+    Context context1(system, integrator1, platform);
+    context1.setPositions(positions);
+    State state1 = context1.getState(State::Forces | State::Energy);
+    VerletIntegrator integrator2(0.01);
+    string deviceIndex = platform.getPropertyValue(context1, "DeviceIndex");
+    map<string, string> props;
+    props["DeviceIndex"] = deviceIndex+","+deviceIndex;
+    Context context2(system, integrator2, platform, props);
+    context2.setPositions(positions);
+    State state2 = context2.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
+
+    // Try updating some parameters and see if they still match.
+
+    vector<double> params;
+    for (int i = 95; i < 102; i++) {
+        int p1, p2;
+        force->getBondParameters(i, p1, p2, params);
+        force->setBondParameters(i, p1, p2, {2.0, 1.2});
+    }
+    force->updateParametersInContext(context1);
+    force->updateParametersInContext(context2);
+    State state3 = context1.getState(State::Energy);
+    State state4 = context2.getState(State::Energy);
+    ASSERT_EQUAL_TOL(state3.getPotentialEnergy(), state4.getPotentialEnergy(), 1e-5);
+    ASSERT(fabs(state1.getPotentialEnergy()-state3.getPotentialEnergy()) > 0.1);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {

tests/TestCustomCompoundBondForce.h

@@ -483,6 +483,39 @@ void testEnergyParameterDerivatives() {
     }
 }
 
+void testParallelComputation() {
+    System system;
+    const int numParticles = 200;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    CustomCompoundBondForce* force = new CustomCompoundBondForce(2, ("(distance(p1,p2)-1.1)^2"));
+    vector<int> particles(2);
+    vector<double> params;
+    for (int i = 1; i < numParticles; i++) {
+        particles[0] = i-1;
+        particles[1] = i;
+        force->addBond(particles, params);
+    }
+    system.addForce(force);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numParticles; i++)
+        positions[i] = Vec3(i, 0, 0);
+    VerletIntegrator integrator1(0.01);
+    Context context1(system, integrator1, platform);
+    context1.setPositions(positions);
+    State state1 = context1.getState(State::Forces | State::Energy);
+    VerletIntegrator integrator2(0.01);
+    string deviceIndex = platform.getPropertyValue(context1, "DeviceIndex");
+    map<string, string> props;
+    props["DeviceIndex"] = deviceIndex+","+deviceIndex;
+    Context context2(system, integrator2, platform, props);
+    context2.setPositions(positions);
+    State state2 = context2.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {

tests/TestCustomExternalForce.h

@@ -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) 2008-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2024 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -35,6 +35,7 @@
 #include "openmm/System.h"
 #include "openmm/VerletIntegrator.h"
 #include "SimTKOpenMMRealType.h"
+#include "sfmt/SFMT.h"
 #include <iostream>
 #include <vector>
 
@@ -231,6 +232,37 @@ void testAtan2() {
     ASSERT_EQUAL_TOL(atan2(positions[0][0], positions[0][1]), state.getPotentialEnergy(), 1e-5);
 }
 
+void testParallelComputation() {
+    System system;
+    const int numParticles = 200;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    CustomExternalForce* force = new CustomExternalForce("x^2+y^2+z^2");
+    vector<double> params;
+    for (int i = 0; i < numParticles; i++)
+        force->addParticle(i, params);
+    system.addForce(force);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numParticles; i++)
+        positions[i] = Vec3(5*genrand_real2(sfmt), 5*genrand_real2(sfmt), 5*genrand_real2(sfmt));
+    VerletIntegrator integrator1(0.01);
+    Context context1(system, integrator1, platform);
+    context1.setPositions(positions);
+    State state1 = context1.getState(State::Forces | State::Energy);
+    VerletIntegrator integrator2(0.01);
+    string deviceIndex = platform.getPropertyValue(context1, "DeviceIndex");
+    map<string, string> props;
+    props["DeviceIndex"] = deviceIndex+","+deviceIndex;
+    Context context2(system, integrator2, platform, props);
+    context2.setPositions(positions);
+    State state2 = context2.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {

tests/TestCustomNonbondedForce.h

@@ -1508,6 +1508,56 @@ void testComputedValues(int mode) {
     }
 }
 
+void testParallelComputation() {
+    System system;
+    const int numParticles = 200;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    CustomNonbondedForce* force = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6); sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)");
+    force->addPerParticleParameter("sigma");
+    force->addPerParticleParameter("eps");
+    for (int i = 0; i < numParticles; i++)
+        force->addParticle({0.5, 1.0});
+    system.addForce(force);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numParticles; i++)
+        positions[i] = Vec3(5*genrand_real2(sfmt), 5*genrand_real2(sfmt), 5*genrand_real2(sfmt));
+    for (int i = 0; i < numParticles; ++i)
+        for (int j = 0; j < i; ++j) {
+            Vec3 delta = positions[i]-positions[j];
+            if (delta.dot(delta) < 0.1)
+                force->addExclusion(i, j);
+        }
+    VerletIntegrator integrator1(0.01);
+    Context context1(system, integrator1, platform);
+    context1.setPositions(positions);
+    State state1 = context1.getState(State::Forces | State::Energy);
+    VerletIntegrator integrator2(0.01);
+    string deviceIndex = platform.getPropertyValue(context1, "DeviceIndex");
+    map<string, string> props;
+    props["DeviceIndex"] = deviceIndex+","+deviceIndex;
+    Context context2(system, integrator2, platform, props);
+    context2.setPositions(positions);
+    State state2 = context2.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
+
+    // Try updating some parameters and see if they still match.
+
+    vector<double> params;
+    for (int i = 95; i < 102; i++)
+        force->setParticleParameters(i, {0.6, 2.0});
+    force->updateParametersInContext(context1);
+    force->updateParametersInContext(context2);
+    State state3 = context1.getState(State::Energy);
+    State state4 = context2.getState(State::Energy);
+    ASSERT_EQUAL_TOL(state3.getPotentialEnergy(), state4.getPotentialEnergy(), 1e-5);
+    ASSERT(fabs(state1.getPotentialEnergy()-state3.getPotentialEnergy()) > 0.1);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {

tests/TestCustomTorsionForce.h

@@ -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) 2008-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2024 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -259,6 +259,51 @@ void testEnergyParameterDerivatives() {
     }
 }
 
+void testParallelComputation() {
+    System system;
+    const int numParticles = 200;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    CustomTorsionForce* force = new CustomTorsionForce("k*sin(theta-theta0)");
+    force->addPerTorsionParameter("k");
+    force->addPerTorsionParameter("theta0");
+    for (int i = 3; i < numParticles; i++)
+        force->addTorsion(i-3, i-2, i-1, i, {1.0, 1.1});
+    system.addForce(force);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numParticles; i++)
+        positions[i] = Vec3(i, i%2, i%3);
+    VerletIntegrator integrator1(0.01);
+    Context context1(system, integrator1, platform);
+    context1.setPositions(positions);
+    State state1 = context1.getState(State::Forces | State::Energy);
+    VerletIntegrator integrator2(0.01);
+    string deviceIndex = platform.getPropertyValue(context1, "DeviceIndex");
+    map<string, string> props;
+    props["DeviceIndex"] = deviceIndex+","+deviceIndex;
+    Context context2(system, integrator2, platform, props);
+    context2.setPositions(positions);
+    State state2 = context2.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-5);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
+
+    // Try updating some parameters and see if they still match.
+
+    vector<double> params;
+    for (int i = 95; i < 102; i++) {
+        int p1, p2, p3, p4;
+        force->getTorsionParameters(i, p1, p2, p3, p4, params);
+        force->setTorsionParameters(i, p1, p2, p3, p4, {2.0, 1.2});
+    }
+    force->updateParametersInContext(context1);
+    force->updateParametersInContext(context2);
+    State state3 = context1.getState(State::Energy);
+    State state4 = context2.getState(State::Energy);
+    ASSERT_EQUAL_TOL(state3.getPotentialEnergy(), state4.getPotentialEnergy(), 1e-5);
+    ASSERT(fabs(state1.getPotentialEnergy()-state3.getPotentialEnergy()) > 0.1);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {