Continuing to implement new CUDA platform: HarmonicAngleForce,...

Continuing to implement new CUDA platform: HarmonicAngleForce, PeriodicTorsionForce, RBTorsionForce, CMAPTorsionForce

platforms/cuda2/src/CudaBondedUtilities.cpp

@@ -74,7 +74,7 @@ void CudaBondedUtilities::initialize(const System& system) {
         int numAtoms = forceAtoms[i][0].size();
         int startAtom = 0;
         while (startAtom < numAtoms) {
-            int width = max(numAtoms-startAtom, 4);
+            int width = min(numAtoms-startAtom, 4);
             if (width == 3)
                 width = 4;
             vector<unsigned int> indexVec(width*numBonds);
@@ -136,7 +136,7 @@ string CudaBondedUtilities::createForceSource(int forceIndex, int numBonds, int
         int atomsToLoad = min(indexWidth, numAtoms-startAtom);
         for (int j = 0; j < atomsToLoad; j++) {
             s<<"    unsigned int atom"<<(startAtom+j+1)<<" = atoms"<<i<<suffix[j]<<";\n";
-            s<<"    real4 pos"<<(j+1)<<" = posq[atom"<<(j+1)<<"];\n";
+            s<<"    real4 pos"<<(startAtom+j+1)<<" = posq[atom"<<(startAtom+j+1)<<"];\n";
         }
         startAtom += indexWidth;
     }

platforms/cuda2/src/CudaKernelFactory.cpp

@@ -80,16 +80,16 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
         return new CudaCalcHarmonicBondForceKernel(name, platform, cu, context.getSystem());
 //    if (name == CalcCustomBondForceKernel::Name())
 //        return new CudaCalcCustomBondForceKernel(name, platform, cu, context.getSystem());
-//    if (name == CalcHarmonicAngleForceKernel::Name())
-//        return new CudaCalcHarmonicAngleForceKernel(name, platform, cu, context.getSystem());
+    if (name == CalcHarmonicAngleForceKernel::Name())
+        return new CudaCalcHarmonicAngleForceKernel(name, platform, cu, context.getSystem());
 //    if (name == CalcCustomAngleForceKernel::Name())
 //        return new CudaCalcCustomAngleForceKernel(name, platform, cu, context.getSystem());
-//    if (name == CalcPeriodicTorsionForceKernel::Name())
-//        return new CudaCalcPeriodicTorsionForceKernel(name, platform, cu, context.getSystem());
-//    if (name == CalcRBTorsionForceKernel::Name())
-//        return new CudaCalcRBTorsionForceKernel(name, platform, cu, context.getSystem());
-//    if (name == CalcCMAPTorsionForceKernel::Name())
-//        return new CudaCalcCMAPTorsionForceKernel(name, platform, cu, context.getSystem());
+    if (name == CalcPeriodicTorsionForceKernel::Name())
+        return new CudaCalcPeriodicTorsionForceKernel(name, platform, cu, context.getSystem());
+    if (name == CalcRBTorsionForceKernel::Name())
+        return new CudaCalcRBTorsionForceKernel(name, platform, cu, context.getSystem());
+    if (name == CalcCMAPTorsionForceKernel::Name())
+        return new CudaCalcCMAPTorsionForceKernel(name, platform, cu, context.getSystem());
 //    if (name == CalcCustomTorsionForceKernel::Name())
 //        return new CudaCalcCustomTorsionForceKernel(name, platform, cu, context.getSystem());
 //    if (name == CalcNonbondedForceKernel::Name())

platforms/cuda2/src/CudaKernels.h

@@ -299,47 +299,47 @@ private:
 //    std::vector<std::string> globalParamNames;
 //    std::vector<cl_float> globalParamValues;
 //};
-//
-///**
-// * This kernel is invoked by HarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
-// */
-//class CudaCalcHarmonicAngleForceKernel : public CalcHarmonicAngleForceKernel {
-//public:
-//    CudaCalcHarmonicAngleForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : CalcHarmonicAngleForceKernel(name, platform),
-//            hasInitializedKernel(false), cu(cu), system(system), params(NULL) {
-//    }
-//    ~CudaCalcHarmonicAngleForceKernel();
-//    /**
-//     * Initialize the kernel.
-//     *
-//     * @param system     the System this kernel will be applied to
-//     * @param force      the HarmonicAngleForce this kernel will be used for
-//     */
-//    void initialize(const System& system, const HarmonicAngleForce& force);
-//    /**
-//     * Execute the kernel to calculate the forces and/or energy.
-//     *
-//     * @param context        the context in which to execute this kernel
-//     * @param includeForces  true if forces should be calculated
-//     * @param includeEnergy  true if the energy should be calculated
-//     * @return the potential energy due to the force
-//     */
-//    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
-//    /**
-//     * Copy changed parameters over to a context.
-//     *
-//     * @param context    the context to copy parameters to
-//     * @param force      the HarmonicAngleForce to copy the parameters from
-//     */
-//    void copyParametersToContext(ContextImpl& context, const HarmonicAngleForce& force);
-//private:
-//    int numAngles;
-//    bool hasInitializedKernel;
-//    CudaContext& cu;
-//    System& system;
-//    CudaArray<mm_float2>* params;
-//};
-//
+
+/**
+ * This kernel is invoked by HarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CudaCalcHarmonicAngleForceKernel : public CalcHarmonicAngleForceKernel {
+public:
+    CudaCalcHarmonicAngleForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : CalcHarmonicAngleForceKernel(name, platform),
+            hasInitializedKernel(false), cu(cu), system(system), params(NULL) {
+    }
+    ~CudaCalcHarmonicAngleForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the HarmonicAngleForce this kernel will be used for
+     */
+    void initialize(const System& system, const HarmonicAngleForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the HarmonicAngleForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const HarmonicAngleForce& force);
+private:
+    int numAngles;
+    bool hasInitializedKernel;
+    CudaContext& cu;
+    System& system;
+    CudaArray* params;
+};
+
 ///**
 // * This kernel is invoked by CustomAngleForce to calculate the forces acting on the system and the energy of the system.
 // */
@@ -382,122 +382,123 @@ private:
 //    std::vector<std::string> globalParamNames;
 //    std::vector<cl_float> globalParamValues;
 //};
-//
-///**
-// * This kernel is invoked by PeriodicTorsionForce to calculate the forces acting on the system and the energy of the system.
-// */
-//class CudaCalcPeriodicTorsionForceKernel : public CalcPeriodicTorsionForceKernel {
-//public:
-//    CudaCalcPeriodicTorsionForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : CalcPeriodicTorsionForceKernel(name, platform),
-//            hasInitializedKernel(false), cu(cu), system(system), params(NULL) {
-//    }
-//    ~CudaCalcPeriodicTorsionForceKernel();
-//    /**
-//     * Initialize the kernel.
-//     *
-//     * @param system     the System this kernel will be applied to
-//     * @param force      the PeriodicTorsionForce this kernel will be used for
-//     */
-//    void initialize(const System& system, const PeriodicTorsionForce& force);
-//    /**
-//     * Execute the kernel to calculate the forces and/or energy.
-//     *
-//     * @param context        the context in which to execute this kernel
-//     * @param includeForces  true if forces should be calculated
-//     * @param includeEnergy  true if the energy should be calculated
-//     * @return the potential energy due to the force
-//     */
-//    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
-//    /**
-//     * Copy changed parameters over to a context.
-//     *
-//     * @param context    the context to copy parameters to
-//     * @param force      the PeriodicTorsionForce to copy the parameters from
-//     */
-//    void copyParametersToContext(ContextImpl& context, const PeriodicTorsionForce& force);
-//private:
-//    int numTorsions;
-//    bool hasInitializedKernel;
-//    CudaContext& cu;
-//    System& system;
-//    CudaArray<mm_float4>* params;
-//};
-//
-///**
-// * This kernel is invoked by RBTorsionForce to calculate the forces acting on the system and the energy of the system.
-// */
-//class CudaCalcRBTorsionForceKernel : public CalcRBTorsionForceKernel {
-//public:
-//    CudaCalcRBTorsionForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : CalcRBTorsionForceKernel(name, platform),
-//            hasInitializedKernel(false), cu(cu), system(system), params(NULL) {
-//    }
-//    ~CudaCalcRBTorsionForceKernel();
-//    /**
-//     * Initialize the kernel.
-//     *
-//     * @param system     the System this kernel will be applied to
-//     * @param force      the RBTorsionForce this kernel will be used for
-//     */
-//    void initialize(const System& system, const RBTorsionForce& force);
-//    /**
-//     * Execute the kernel to calculate the forces and/or energy.
-//     *
-//     * @param context        the context in which to execute this kernel
-//     * @param includeForces  true if forces should be calculated
-//     * @param includeEnergy  true if the energy should be calculated
-//     * @return the potential energy due to the force
-//     */
-//    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
-//    /**
-//     * Copy changed parameters over to a context.
-//     *
-//     * @param context    the context to copy parameters to
-//     * @param force      the RBTorsionForce to copy the parameters from
-//     */
-//    void copyParametersToContext(ContextImpl& context, const RBTorsionForce& force);
-//private:
-//    int numTorsions;
-//    bool hasInitializedKernel;
-//    CudaContext& cu;
-//    System& system;
-//    CudaArray<mm_float8>* params;
-//};
-//
-///**
-// * This kernel is invoked by CMAPTorsionForce to calculate the forces acting on the system and the energy of the system.
-// */
-//class CudaCalcCMAPTorsionForceKernel : public CalcCMAPTorsionForceKernel {
-//public:
-//    CudaCalcCMAPTorsionForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : CalcCMAPTorsionForceKernel(name, platform),
-//            hasInitializedKernel(false), cu(cu), system(system), coefficients(NULL), mapPositions(NULL), torsionMaps(NULL) {
-//    }
-//    ~CudaCalcCMAPTorsionForceKernel();
-//    /**
-//     * Initialize the kernel.
-//     *
-//     * @param system     the System this kernel will be applied to
-//     * @param force      the CMAPTorsionForce this kernel will be used for
-//     */
-//    void initialize(const System& system, const CMAPTorsionForce& force);
-//    /**
-//     * Execute the kernel to calculate the forces and/or energy.
-//     *
-//     * @param context        the context in which to execute this kernel
-//     * @param includeForces  true if forces should be calculated
-//     * @param includeEnergy  true if the energy should be calculated
-//     * @return the potential energy due to the force
-//     */
-//    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
-//private:
-//    int numTorsions;
-//    bool hasInitializedKernel;
-//    CudaContext& cu;
-//    System& system;
-//    CudaArray<mm_float4>* coefficients;
-//    CudaArray<mm_int2>* mapPositions;
-//    CudaArray<cl_int>* torsionMaps;
-//};
-//
+
+/**
+ * This kernel is invoked by PeriodicTorsionForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CudaCalcPeriodicTorsionForceKernel : public CalcPeriodicTorsionForceKernel {
+public:
+    CudaCalcPeriodicTorsionForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : CalcPeriodicTorsionForceKernel(name, platform),
+            hasInitializedKernel(false), cu(cu), system(system), params(NULL) {
+    }
+    ~CudaCalcPeriodicTorsionForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the PeriodicTorsionForce this kernel will be used for
+     */
+    void initialize(const System& system, const PeriodicTorsionForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the PeriodicTorsionForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const PeriodicTorsionForce& force);
+private:
+    int numTorsions;
+    bool hasInitializedKernel;
+    CudaContext& cu;
+    System& system;
+    CudaArray* params;
+};
+
+/**
+ * This kernel is invoked by RBTorsionForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CudaCalcRBTorsionForceKernel : public CalcRBTorsionForceKernel {
+public:
+    CudaCalcRBTorsionForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : CalcRBTorsionForceKernel(name, platform),
+            hasInitializedKernel(false), cu(cu), system(system), params1(NULL), params2(NULL) {
+    }
+    ~CudaCalcRBTorsionForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the RBTorsionForce this kernel will be used for
+     */
+    void initialize(const System& system, const RBTorsionForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the RBTorsionForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const RBTorsionForce& force);
+private:
+    int numTorsions;
+    bool hasInitializedKernel;
+    CudaContext& cu;
+    System& system;
+    CudaArray* params1;
+    CudaArray* params2;
+};
+
+/**
+ * This kernel is invoked by CMAPTorsionForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CudaCalcCMAPTorsionForceKernel : public CalcCMAPTorsionForceKernel {
+public:
+    CudaCalcCMAPTorsionForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : CalcCMAPTorsionForceKernel(name, platform),
+            hasInitializedKernel(false), cu(cu), system(system), coefficients(NULL), mapPositions(NULL), torsionMaps(NULL) {
+    }
+    ~CudaCalcCMAPTorsionForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CMAPTorsionForce this kernel will be used for
+     */
+    void initialize(const System& system, const CMAPTorsionForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+private:
+    int numTorsions;
+    bool hasInitializedKernel;
+    CudaContext& cu;
+    System& system;
+    CudaArray* coefficients;
+    CudaArray* mapPositions;
+    CudaArray* torsionMaps;
+};
+
 ///**
 // * This kernel is invoked by CustomTorsionForce to calculate the forces acting on the system and the energy of the system.
 // */

platforms/cuda2/src/kernels/angleForce.cu

+real3 v0 = make_real3(pos2.x-pos1.x, pos2.y-pos1.y, pos2.z-pos1.z);
+real3 v1 = make_real3(pos2.x-pos3.x, pos2.y-pos3.y, pos2.z-pos3.z);
+real3 cp = cross(v0, v1);
+real rp = cp.x*cp.x + cp.y*cp.y + cp.z*cp.z;
+rp = max(SQRT(rp), (real) 1.0e-06f);
+real r21 = v0.x*v0.x + v0.y*v0.y + v0.z*v0.z;
+real r23 = v1.x*v1.x + v1.y*v1.y + v1.z*v1.z;
+real dot = v0.x*v1.x + v0.y*v1.y + v0.z*v1.z;
+real cosine = min(max(dot*RSQRT(r21*r23), (real) -1), (real) 1);
+real theta = ACOS(cosine);
+COMPUTE_FORCE
+real3 force1 = cross(v0, cp)*(dEdAngle/(r21*rp));
+real3 force3 = cross(cp, v1)*(dEdAngle/(r23*rp));
+real3 force2 = -force1-force3;

platforms/cuda2/src/kernels/cmapTorsionForce.cu

+const real PI = (real) 3.14159265358979323846;
+
+// Compute the first angle.
+
+real3 v0a = make_real3(pos1.x-pos2.x, pos1.y-pos2.y, pos1.z-pos2.z);
+real3 v1a = make_real3(pos3.x-pos2.x, pos3.y-pos2.y, pos3.z-pos2.z);
+real3 v2a = make_real3(pos3.x-pos4.x, pos3.y-pos4.y, pos3.z-pos4.z);
+real3 cp0a = cross(v0a, v1a);
+real3 cp1a = cross(v1a, v2a);
+real cosangle = dot(normalize(cp0a), normalize(cp1a));
+real angleA;
+if (cosangle > 0.99f || cosangle < -0.99f) {
+    // We're close to the singularity in acos(), so take the cross product and use asin() instead.
+
+    real3 cross_prod = cross(cp0a, cp1a);
+    real scale = dot(cp0a, cp0a)*dot(cp1a, cp1a);
+    angleA = ASIN(SQRT(dot(cross_prod, cross_prod)/scale));
+    if (cosangle < 0.0f)
+        angleA = PI-angleA;
+}
+else
+   angleA = ACOS(cosangle);
+angleA = (dot(v0a, cp1a) >= 0 ? angleA : -angleA);
+angleA = fmod(angleA+2.0f*PI, 2.0f*PI);
+
+// Compute the second angle.
+
+real3 v0b = make_real3(pos5.x-pos6.x, pos5.y-pos6.y, pos5.z-pos6.z);
+real3 v1b = make_real3(pos7.x-pos6.x, pos7.y-pos6.y, pos7.z-pos6.z);
+real3 v2b = make_real3(pos7.x-pos8.x, pos7.y-pos8.y, pos7.z-pos8.z);
+real3 cp0b = cross(v0b, v1b);
+real3 cp1b = cross(v1b, v2b);
+cosangle = dot(normalize(cp0b), normalize(cp1b));
+real angleB;
+if (cosangle > 0.99f || cosangle < -0.99f) {
+    // We're close to the singularity in acos(), so take the cross product and use asin() instead.
+
+    real3 cross_prod = cross(cp0b, cp1b);
+    real scale = dot(cp0b, cp0b)*dot(cp1b, cp1b);
+    angleB = ASIN(SQRT(dot(cross_prod, cross_prod)/scale));
+    if (cosangle < 0.0f)
+        angleB = PI-angleB;
+}
+else
+   angleB = ACOS(cosangle);
+angleB = (dot(v0b, cp1b) >= 0 ? angleB : -angleB);
+angleB = fmod(angleB+2.0f*PI, 2.0f*PI);
+
+// Identify which patch this is in.
+
+int2 pos = MAP_POS[MAPS[index]];
+int size = pos.y;
+real delta = 2*PI/size;
+int s = (int) (angleA/delta);
+int t = (int) (angleB/delta);
+float4 c[4];
+int coeffIndex = pos.x+4*(s+size*t);
+c[0] = COEFF[coeffIndex];
+c[1] = COEFF[coeffIndex+1];
+c[2] = COEFF[coeffIndex+2];
+c[3] = COEFF[coeffIndex+3];
+real da = angleA/delta-s;
+real db = angleB/delta-t;
+
+// Evaluate the spline to determine the energy and gradients.
+
+real torsionEnergy = 0.0f;
+real dEdA = 0.0f;
+real dEdB = 0.0f;
+torsionEnergy = da*torsionEnergy + ((c[3].w*db + c[3].z)*db + c[3].y)*db + c[3].x;
+dEdA = db*dEdA + (3.0f*c[3].w*da + 2.0f*c[2].w)*da + c[1].w;
+dEdB = da*dEdB + (3.0f*c[3].w*db + 2.0f*c[3].z)*db + c[3].y;
+torsionEnergy = da*torsionEnergy + ((c[2].w*db + c[2].z)*db + c[2].y)*db + c[2].x;
+dEdA = db*dEdA + (3.0f*c[3].z*da + 2.0f*c[2].z)*da + c[1].z;
+dEdB = da*dEdB + (3.0f*c[2].w*db + 2.0f*c[2].z)*db + c[2].y;
+torsionEnergy = da*torsionEnergy + ((c[1].w*db + c[1].z)*db + c[1].y)*db + c[1].x;
+dEdA = db*dEdA + (3.0f*c[3].y*da + 2.0f*c[2].y)*da + c[1].y;
+dEdB = da*dEdB + (3.0f*c[1].w*db + 2.0f*c[1].z)*db + c[1].y;
+torsionEnergy = da*torsionEnergy + ((c[0].w*db + c[0].z)*db + c[0].y)*db + c[0].x;
+dEdA = db*dEdA + (3.0f*c[3].x*da + 2.0f*c[2].x)*da + c[1].x;
+dEdB = da*dEdB + (3.0f*c[0].w*db + 2.0f*c[0].z)*db + c[0].y;
+dEdA /= delta;
+dEdB /= delta;
+energy += torsionEnergy;
+
+// Apply the force to the first torsion.
+
+real normCross1 = dot(cp0a, cp0a);
+real normSqrBC = dot(v1a, v1a);
+real normBC = SQRT(normSqrBC);
+real normCross2 = dot(cp1a, cp1a);
+real dp = RECIP(normSqrBC);
+real4 ff = make_real4((-dEdA*normBC)/normCross1, dot(v0a, v1a)*dp, dot(v2a, v1a)*dp, (dEdA*normBC)/normCross2);
+real3 force1 = ff.x*cp0a;
+real3 force4 = ff.w*cp1a;
+real3 d = ff.y*force1 - ff.z*force4;
+real3 force2 = d-force1;
+real3 force3 = -d-force4;
+
+// Apply the force to the second torsion.
+
+normCross1 = dot(cp0b, cp0b);
+normSqrBC = dot(v1b, v1b);
+normBC = SQRT(normSqrBC);
+normCross2 = dot(cp1b, cp1b);
+dp = RECIP(normSqrBC);
+ff = make_real4((-dEdB*normBC)/normCross1, dot(v0b, v1b)*dp, dot(v2b, v1b)*dp, (dEdB*normBC)/normCross2);
+real3 force5 = ff.x*cp0b;
+real3 force8 = ff.w*cp1b;
+d = ff.y*force5 - ff.z*force8;
+real3 force6 = d-force5;
+real3 force7 = -d-force8;

platforms/cuda2/src/kernels/harmonicAngleForce.cu

+float2 angleParams = PARAMS[index];
+real deltaIdeal = theta-angleParams.x;
+energy += 0.5f*angleParams.y*deltaIdeal*deltaIdeal;
+real dEdAngle = angleParams.y*deltaIdeal;

platforms/cuda2/src/kernels/periodicTorsionForce.cu

+real4 torsionParams = PARAMS[index];
+real deltaAngle = torsionParams.z*theta-torsionParams.y;
+energy += torsionParams.x*(1.0f+COS(deltaAngle));
+real sinDeltaAngle = SIN(deltaAngle);
+real dEdAngle = -torsionParams.x*torsionParams.z*sinDeltaAngle;

platforms/cuda2/src/kernels/rbTorsionForce.cu

+real4 torsionParams1 = PARAMS1[index];
+real2 torsionParams2 = PARAMS2[index];
+if (theta < 0)
+    theta += PI;
+else
+    theta -= PI;
+cosangle = -cosangle;
+real cosFactor = cosangle;
+real dEdAngle = -torsionParams1.y;
+real rbEnergy = torsionParams1.x;
+rbEnergy += torsionParams1.y*cosFactor;
+dEdAngle -= 2.0f*torsionParams1.z*cosFactor;
+cosFactor *= cosangle;
+dEdAngle -= 3.0f*torsionParams1.w*cosFactor;
+rbEnergy += torsionParams1.z*cosFactor;
+cosFactor *= cosangle;
+dEdAngle -= 4.0f*torsionParams2.x*cosFactor;
+rbEnergy += torsionParams1.w*cosFactor;
+cosFactor *= cosangle;
+dEdAngle -= 5.0f*torsionParams2.y*cosFactor;
+rbEnergy += torsionParams2.x*cosFactor;
+rbEnergy += torsionParams2.y*cosFactor*cosangle;
+energy += rbEnergy;
+dEdAngle *= SIN(theta);

platforms/cuda2/src/kernels/torsionForce.cu

+const real PI = (real) 3.14159265358979323846;
+real3 v0 = make_real3(pos1.x-pos2.x, pos1.y-pos2.y, pos1.z-pos2.z);
+real3 v1 = make_real3(pos3.x-pos2.x, pos3.y-pos2.y, pos3.z-pos2.z);
+real3 v2 = make_real3(pos3.x-pos4.x, pos3.y-pos4.y, pos3.z-pos4.z);
+real3 cp0 = cross(v0, v1);
+real3 cp1 = cross(v1, v2);
+real cosangle = dot(normalize(cp0), normalize(cp1));
+real theta;
+if (cosangle > 0.99f || cosangle < -0.99f) {
+    // We're close to the singularity in acos(), so take the cross product and use asin() instead.
+
+    real3 cross_prod = cross(cp0, cp1);
+    real scale = dot(cp0, cp0)*dot(cp1, cp1);
+    theta = ASIN(SQRT(dot(cross_prod, cross_prod)/scale));
+    if (cosangle < 0)
+        theta = PI-theta;
+}
+else
+   theta = acos(cosangle);
+theta = (dot(v0, cp1) >= 0 ? theta : -theta);
+COMPUTE_FORCE
+real normCross1 = dot(cp0, cp0);
+real normSqrBC = dot(v1, v1);
+real normBC = sqrt(normSqrBC);
+real normCross2 = dot(cp1, cp1);
+real dp = RECIP(normSqrBC);
+real4 ff = make_real4((-dEdAngle*normBC)/normCross1, dot(v0, v1)*dp, dot(v2, v1)*dp, (dEdAngle*normBC)/normCross2);
+real3 force1 = ff.x*cp0;
+real3 force4 = ff.w*cp1;
+real3 s = ff.y*force1 - ff.z*force4;
+real3 force2 = s-force1;
+real3 force3 = -s-force4;

platforms/cuda2/src/kernels/vectorOps.cu

@@ -493,3 +493,57 @@ inline __device__ void operator*=(double3& a, double b) {
 inline __device__ void operator*=(double4& a, double b) {
     a.x *= b; a.y *= b; a.z *= b; a.w *= b;
 }
+
+// Dot product
+
+inline __device__ float dot(float3 a, float3 b) {
+    return a.x*b.x+a.y*b.y+a.z*b.z;
+}
+
+inline __device__ double dot(double3 a, double3 b) {
+    return a.x*b.x+a.y*b.y+a.z*b.z;
+}
+
+// Cross product
+
+inline __device__ float3 cross(float3 a, float3 b) {
+    return make_float3(a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x);
+}
+
+inline __device__ float3 cross(float4 a, float4 b) {
+    return make_float3(a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x);
+}
+
+inline __device__ double3 cross(double3 a, double3 b) {
+    return make_double3(a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x);
+}
+
+inline __device__ double3 cross(double4 a, double4 b) {
+    return make_double3(a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x);
+}
+
+// Normalize a vector
+
+inline __device__ float2 normalize(float2 a) {
+    return a*rsqrtf(a.x*a.x+a.y*a.y);
+}
+
+inline __device__ float3 normalize(float3 a) {
+    return a*rsqrtf(a.x*a.x+a.y*a.y+a.z*a.z);
+}
+
+inline __device__ float4 normalize(float4 a) {
+    return a*rsqrtf(a.x*a.x+a.y*a.y+a.z*a.z+a.w*a.w);
+}
+
+inline __device__ double2 normalize(double2 a) {
+    return a*rsqrt(a.x*a.x+a.y*a.y);
+}
+
+inline __device__ double3 normalize(double3 a) {
+    return a*rsqrt(a.x*a.x+a.y*a.y+a.z*a.z);
+}
+
+inline __device__ double4 normalize(double4 a) {
+    return a*rsqrt(a.x*a.x+a.y*a.y+a.z*a.z+a.w*a.w);
+}

platforms/cuda2/tests/TestCudaCMAPTorsionForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2010-2012 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the CUDA implementation of CMAPTorsionForce.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "CudaPlatform.h"
+#include "openmm/CMAPTorsionForce.h"
+#include "openmm/PeriodicTorsionForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testCMAPTorsions() {
+    const int mapSize = 36;
+
+    // Create two systems: one with a pair of periodic torsions, and one with a CMAP torsion
+    // that approximates the same force.
+
+    CudaPlatform platform;
+    System system1;
+    for (int i = 0; i < 5; i++)
+        system1.addParticle(1.0);
+    PeriodicTorsionForce* periodic = new PeriodicTorsionForce();
+    periodic->addTorsion(0, 1, 2, 3, 2, M_PI/4, 1.5);
+    periodic->addTorsion(1, 2, 3, 4, 3, M_PI/3, 2.0);
+    system1.addForce(periodic);
+    System system2;
+    for (int i = 0; i < 5; i++)
+        system2.addParticle(1.0);
+    CMAPTorsionForce* cmap = new CMAPTorsionForce();
+    vector<double> mapEnergy(mapSize*mapSize);
+    for (int i = 0; i < mapSize; i++) {
+        double angle1 = i*2*M_PI/mapSize;
+        double energy1 = 1.5*(1+cos(2*angle1-M_PI/4));
+        for (int j = 0; j < mapSize; j++) {
+            double angle2 = j*2*M_PI/mapSize;
+            double energy2 = 2.0*(1+cos(3*angle2-M_PI/3));
+            mapEnergy[i+j*mapSize] = energy1+energy2;
+        }
+    }
+    cmap->addMap(mapSize, mapEnergy);
+    cmap->addTorsion(0, 0, 1, 2, 3, 1, 2, 3, 4);
+    system2.addForce(cmap);
+
+    // Set the atoms in various positions, and verify that both systems give equal forces and energy.
+
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(5);
+    VerletIntegrator integrator1(0.01);
+    VerletIntegrator integrator2(0.01);
+    Context c1(system1, integrator1, platform);
+    Context c2(system2, integrator2, platform);
+    for (int i = 0; i < 50; i++) {
+        for (int j = 0; j < (int) positions.size(); j++)
+            positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
+        c1.setPositions(positions);
+        c2.setPositions(positions);
+        State s1 = c1.getState(State::Forces | State::Energy);
+        State s2 = c2.getState(State::Forces | State::Energy);
+        for (int i = 0; i < system1.getNumParticles(); i++)
+            ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 0.05);
+        ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), 1e-3);
+    }
+}
+
+int main() {
+    try {
+        testCMAPTorsions();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/cuda2/tests/TestCudaHarmonicAngleForce.cpp

+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the CUDA implementation of HarmonicAngleForce.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "CudaPlatform.h"
+#include "openmm/HarmonicAngleForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testAngles() {
+    CudaPlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    HarmonicAngleForce* forceField = new HarmonicAngleForce();
+    forceField->addAngle(0, 1, 2, PI_M/3, 1.1);
+    forceField->addAngle(1, 2, 3, PI_M/2, 1.2);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(4);
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    positions[3] = Vec3(2, 1, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    {
+        const vector<Vec3>& forces = state.getForces();
+        double torque1 = 1.1*PI_M/6;
+        double torque2 = 1.2*PI_M/4;
+        ASSERT_EQUAL_VEC(Vec3(torque1, 0, 0), forces[0], TOL);
+        ASSERT_EQUAL_VEC(Vec3(-0.5*torque2, 0.5*torque2, 0), forces[3], TOL); // reduced by sqrt(2) due to the bond length, another sqrt(2) due to the angle
+        ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
+        ASSERT_EQUAL_TOL(0.5*1.1*(PI_M/6)*(PI_M/6) + 0.5*1.2*(PI_M/4)*(PI_M/4), state.getPotentialEnergy(), TOL);
+    }
+    
+    // Try changing the angle parameters and make sure it's still correct.
+    
+    forceField->setAngleParameters(0, 0, 1, 2, PI_M/3.1, 1.3);
+    forceField->setAngleParameters(1, 1, 2, 3, PI_M/2.1, 1.4);
+    forceField->updateParametersInContext(context);
+    state = context.getState(State::Forces | State::Energy);
+    {
+        const vector<Vec3>& forces = state.getForces();
+        double dtheta1 = (PI_M/2)-(PI_M/3.1);
+        double dtheta2 = (3*PI_M/4)-(PI_M/2.1);
+        double torque1 = 1.3*dtheta1;
+        double torque2 = 1.4*dtheta2;
+        ASSERT_EQUAL_VEC(Vec3(torque1, 0, 0), forces[0], TOL);
+        ASSERT_EQUAL_VEC(Vec3(-0.5*torque2, 0.5*torque2, 0), forces[3], TOL); // reduced by sqrt(2) due to the bond length, another sqrt(2) due to the angle
+        ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
+        ASSERT_EQUAL_TOL(0.5*1.3*dtheta1*dtheta1 + 0.5*1.4*dtheta2*dtheta2, state.getPotentialEnergy(), TOL);
+    }
+}
+
+void testParallelComputation() {
+    CudaPlatform platform;
+    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);
+    string deviceIndex = platform.getPropertyValue(context1, CudaPlatform::CudaDeviceIndex());
+    map<string, string> props;
+    props[CudaPlatform::CudaDeviceIndex()] = 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);
+}
+
+int main() {
+    try {
+        testAngles();
+//        testParallelComputation();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/cuda2/tests/TestCudaPeriodicTorsionForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the CUDA implementation of PeriodicTorsionForce.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "CudaPlatform.h"
+#include "openmm/PeriodicTorsionForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testPeriodicTorsions() {
+    CudaPlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    PeriodicTorsionForce* forceField = new PeriodicTorsionForce();
+    forceField->addTorsion(0, 1, 2, 3, 2, PI_M/3, 1.1);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(4);
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    positions[3] = Vec3(1, 0, 2);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    {
+        const vector<Vec3>& forces = state.getForces();
+        double torque = -2*1.1*std::sin(2*PI_M/3);
+        ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
+        ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, 0), forces[3], TOL);
+        ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
+        ASSERT_EQUAL_TOL(1.1*(1+std::cos(2*PI_M/3)), state.getPotentialEnergy(), TOL);
+    }
+    
+    // Try changing the torsion parameters and make sure it's still correct.
+    
+    forceField->setTorsionParameters(0, 0, 1, 2, 3, 3, PI_M/3.2, 1.3);
+    forceField->updateParametersInContext(context);
+    state = context.getState(State::Forces | State::Energy);
+    {
+        const vector<Vec3>& forces = state.getForces();
+        double dtheta = (3*PI_M/2)-(PI_M/3.2);
+        double torque = -3*1.3*std::sin(dtheta);
+        ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
+        ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, 0), forces[3], TOL);
+        ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
+        ASSERT_EQUAL_TOL(1.3*(1+std::cos(dtheta)), state.getPotentialEnergy(), TOL);
+    }
+}
+
+void testParallelComputation() {
+    CudaPlatform platform;
+    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);
+    string deviceIndex = platform.getPropertyValue(context1, CudaPlatform::CudaDeviceIndex());
+    map<string, string> props;
+    props[CudaPlatform::CudaDeviceIndex()] = 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);
+}
+
+int main() {
+    try {
+        testPeriodicTorsions();
+//        testParallelComputation();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+

platforms/cuda2/tests/TestCudaRBTorsionForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the CUDA implementation of RBTorsionForce.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "CudaPlatform.h"
+#include "openmm/RBTorsionForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testRBTorsions() {
+    CudaPlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    RBTorsionForce* forceField = new RBTorsionForce();
+    forceField->addTorsion(0, 1, 2, 3, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(4);
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    positions[3] = Vec3(1, 1, 1);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    {
+        const vector<Vec3>& forces = state.getForces();
+        double psi = 0.25*PI_M - PI_M;
+        double torque = 0.0;
+        for (int i = 1; i < 6; ++i) {
+            double c = 0.1*(i+1);
+            torque += -c*i*std::pow(std::cos(psi), i-1)*std::sin(psi);
+        }
+        ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
+        ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, -0.5*torque), forces[3], TOL);
+        ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
+        double energy = 0.0;
+        for (int i = 0; i < 6; ++i) {
+            double c = 0.1*(i+1);
+            energy += c*std::pow(std::cos(psi), i);
+        }
+        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
+    }
+    
+    // Try changing the torsion parameters and make sure it's still correct.
+    
+    forceField->setTorsionParameters(0, 0, 1, 2, 3, 0.11, 0.22, 0.33, 0.44, 0.55, 0.66);
+    forceField->updateParametersInContext(context);
+    state = context.getState(State::Forces | State::Energy);
+    {
+        const vector<Vec3>& forces = state.getForces();
+        double psi = 0.25*PI_M - PI_M;
+        double torque = 0.0;
+        for (int i = 1; i < 6; ++i) {
+            double c = 0.11*(i+1);
+            torque += -c*i*std::pow(std::cos(psi), i-1)*std::sin(psi);
+        }
+        ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
+        ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, -0.5*torque), forces[3], TOL);
+        ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
+        double energy = 0.0;
+        for (int i = 0; i < 6; ++i) {
+            double c = 0.11*(i+1);
+            energy += c*std::pow(std::cos(psi), i);
+        }
+        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
+    }
+}
+
+void testParallelComputation() {
+    CudaPlatform platform;
+    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);
+    string deviceIndex = platform.getPropertyValue(context1, CudaPlatform::CudaDeviceIndex());
+    map<string, string> props;
+    props[CudaPlatform::CudaDeviceIndex()] = 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);
+}
+
+int main() {
+    try {
+        testRBTorsions();
+//        testParallelComputation();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+