Created OpenCL implementation of CMAPTorsionForce

platforms/opencl/src/OpenCLContext.h

@@ -74,6 +74,16 @@ struct mm_float8 {
         s0(s0), s1(s1), s2(s2), s3(s3), s4(s4), s5(s5), s6(s6), s7(s7) {
     }
 };
+struct mm_float16 {
+    cl_float s0, s1, s2, s3, s4, s5, s6, s7, s8, s9, s10, s11, s12, s13, s14, s15;
+    mm_float16() {
+    }
+    mm_float16(cl_float s0, cl_float s1, cl_float s2, cl_float s3, cl_float s4, cl_float s5, cl_float s6, cl_float s7,
+            cl_float s8, cl_float s9, cl_float s10, cl_float s11, cl_float s12, cl_float s13, cl_float s14, cl_float s15) :
+        s0(s0), s1(s1), s2(s2), s3(s3), s4(s4), s5(s5), s6(s6), s7(s7),
+        s8(s8), s9(s9), s10(s10), s11(s11), s12(s12), s13(s13), s14(s14), s15(15) {
+    }
+};
 struct mm_int2 {
     cl_int x, y;
     mm_int2() {
@@ -96,6 +106,16 @@ struct mm_int8 {
         s0(s0), s1(s1), s2(s2), s3(s3), s4(s4), s5(s5), s6(s6), s7(s7) {
     }
 };
+struct mm_int16 {
+    cl_int s0, s1, s2, s3, s4, s5, s6, s7, s8, s9, s10, s11, s12, s13, s14, s15;
+    mm_int16() {
+    }
+    mm_int16(cl_int s0, cl_int s1, cl_int s2, cl_int s3, cl_int s4, cl_int s5, cl_int s6, cl_int s7,
+            cl_int s8, cl_int s9, cl_int s10, cl_int s11, cl_int s12, cl_int s13, cl_int s14, cl_int s15) :
+        s0(s0), s1(s1), s2(s2), s3(s3), s4(s4), s5(s5), s6(s6), s7(s7),
+        s8(s8), s9(s9), s10(s10), s11(s11), s12(s12), s13(s13), s14(s14), s15(15) {
+    }
+};
 
 /**
  * This class contains the information associated with a Context by the OpenCL Platform.

platforms/opencl/src/OpenCLKernelFactory.cpp

@@ -52,6 +52,8 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo
         return new OpenCLCalcPeriodicTorsionForceKernel(name, platform, cl, context.getSystem());
     if (name == CalcRBTorsionForceKernel::Name())
         return new OpenCLCalcRBTorsionForceKernel(name, platform, cl, context.getSystem());
+    if (name == CalcCMAPTorsionForceKernel::Name())
+        return new OpenCLCalcCMAPTorsionForceKernel(name, platform, cl, context.getSystem());
     if (name == CalcCustomTorsionForceKernel::Name())
         return new OpenCLCalcCustomTorsionForceKernel(name, platform, cl, context.getSystem());
     if (name == CalcNonbondedForceKernel::Name())

platforms/opencl/src/OpenCLKernels.cpp

@@ -28,6 +28,7 @@
 #include "OpenCLForceInfo.h"
 #include "openmm/LangevinIntegrator.h"
 #include "openmm/Context.h"
+#include "openmm/internal/CMAPTorsionForceImpl.h"
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/internal/CustomHbondForceImpl.h"
 #include "openmm/internal/NonbondedForceImpl.h"
@@ -683,7 +684,7 @@ public:
         int particle1, particle2, particle3, particle4, periodicity1, periodicity2;
         double phase1, phase2, k1, k2;
         force.getTorsionParameters(group1, particle1, particle2, particle3, particle4, periodicity1, phase1, k1);
-        force.getTorsionParameters(group1, particle1, particle2, particle3, particle4, periodicity2, phase2, k2);
+        force.getTorsionParameters(group2, particle1, particle2, particle3, particle4, periodicity2, phase2, k2);
         return (periodicity1 == periodicity2 && phase1 == phase2 && k1 == k2);
     }
 private:
@@ -767,7 +768,7 @@ public:
         int particle1, particle2, particle3, particle4;
         double c0a, c0b, c1a, c1b, c2a, c2b, c3a, c3b, c4a, c4b, c5a, c5b;
         force.getTorsionParameters(group1, particle1, particle2, particle3, particle4, c0a, c1a, c2a, c3a, c4a, c5a);
-        force.getTorsionParameters(group1, particle1, particle2, particle3, particle4, c0b, c1b, c2b, c3b, c4b, c5b);
+        force.getTorsionParameters(group2, particle1, particle2, particle3, particle4, c0b, c1b, c2b, c3b, c4b, c5b);
         return (c0a == c0b && c1a == c1b && c2a == c2b && c3a == c3b && c4a == c4b && c5a == c5b);
     }
 private:
@@ -830,6 +831,124 @@ double OpenCLCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
     return 0.0;
 }
 
+class OpenCLCMAPTorsionForceInfo : public OpenCLForceInfo {
+public:
+    OpenCLCMAPTorsionForceInfo(int requiredBuffers, const CMAPTorsionForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
+    }
+    int getNumParticleGroups() {
+        return force.getNumTorsions();
+    }
+    void getParticlesInGroup(int index, std::vector<int>& particles) {
+        int map, a1, a2, a3, a4, b1, b2, b3, b4;
+        force.getTorsionParameters(index, map, a1, a2, a3, a4, b1, b2, b3, b4);
+        particles.resize(8);
+        particles[0] = a1;
+        particles[1] = a2;
+        particles[2] = a3;
+        particles[3] = a4;
+        particles[4] = b1;
+        particles[5] = b2;
+        particles[6] = b3;
+        particles[7] = b4;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        int map1, map2, a1, a2, a3, a4, b1, b2, b3, b4;
+        force.getTorsionParameters(group1, map1, a1, a2, a3, a4, b1, b2, b3, b4);
+        force.getTorsionParameters(group2, map2, a1, a2, a3, a4, b1, b2, b3, b4);
+        return (map1 == map2);
+    }
+private:
+    const CMAPTorsionForce& force;
+};
+
+OpenCLCalcCMAPTorsionForceKernel::~OpenCLCalcCMAPTorsionForceKernel() {
+    if (coefficients != NULL)
+        delete coefficients;
+    if (mapPositions != NULL)
+        delete mapPositions;
+    if (torsionMaps != NULL)
+        delete torsionMaps;
+    if (torsionIndices != NULL)
+        delete torsionIndices;
+}
+
+void OpenCLCalcCMAPTorsionForceKernel::initialize(const System& system, const CMAPTorsionForce& force) {
+    numTorsions = force.getNumTorsions();
+    if (numTorsions == 0)
+        return;
+    int numMaps = force.getNumMaps();
+    vector<mm_float4> coeffVec;
+    vector<mm_int2> mapPositionsVec(numMaps);
+    vector<double> energy;
+    vector<vector<double> > c;
+    int currentPosition = 0;
+    for (int i = 0; i < numMaps; i++) {
+        int size;
+        force.getMapParameters(i, size, energy);
+        CMAPTorsionForceImpl::calcMapDerivatives(size, energy, c);
+        mapPositionsVec[i] = mm_int2(currentPosition, size);
+        currentPosition += 4*size*size;
+        for (int j = 0; j < size*size; j++) {
+            coeffVec.push_back(mm_float4(c[j][0], c[j][1], c[j][2], c[j][3]));
+            coeffVec.push_back(mm_float4(c[j][4], c[j][5], c[j][6], c[j][7]));
+            coeffVec.push_back(mm_float4(c[j][8], c[j][9], c[j][10], c[j][11]));
+            coeffVec.push_back(mm_float4(c[j][12], c[j][13], c[j][14], c[j][15]));
+        }
+    }
+    vector<int> forceBufferCounter(system.getNumParticles(), 0);
+    vector<cl_int> torsionMapsVec(numTorsions);
+    vector<mm_int16> torsionIndicesVec(numTorsions);
+    for (int i = 0; i < numTorsions; i++) {
+        mm_int16& ind = torsionIndicesVec[i];
+        force.getTorsionParameters(i, torsionMapsVec[i], ind.s0, ind.s1, ind.s2, ind.s3, ind.s4, ind.s5, ind.s6, ind.s7);
+        ind.s8 = forceBufferCounter[ind.s0]++;
+        ind.s9 = forceBufferCounter[ind.s1]++;
+        ind.s10 = forceBufferCounter[ind.s2]++;
+        ind.s11 = forceBufferCounter[ind.s3]++;
+        ind.s12 = forceBufferCounter[ind.s4]++;
+        ind.s13 = forceBufferCounter[ind.s5]++;
+        ind.s14 = forceBufferCounter[ind.s6]++;
+        ind.s15 = forceBufferCounter[ind.s7]++;
+    }
+    coefficients = new OpenCLArray<mm_float4>(cl, coeffVec.size(), "cmapTorsionCoefficients");
+    mapPositions = new OpenCLArray<mm_int2>(cl, numMaps, "cmapTorsionMapPositions");
+    torsionMaps = new OpenCLArray<cl_int>(cl, numTorsions, "cmapTorsionMaps");
+    torsionIndices = new OpenCLArray<mm_int16>(cl, numTorsions, "cmapTorsionIndices");
+    coefficients->upload(coeffVec);
+    mapPositions->upload(mapPositionsVec);
+    torsionMaps->upload(torsionMapsVec);
+    torsionIndices->upload(torsionIndicesVec);
+    int maxBuffers = 1;
+    for (int i = 0; i < (int) forceBufferCounter.size(); i++)
+        maxBuffers = max(maxBuffers, forceBufferCounter[i]);
+    cl.addForce(new OpenCLCMAPTorsionForceInfo(maxBuffers, force));
+    cl::Program program = cl.createProgram(OpenCLKernelSources::cmapTorsionForce);
+    kernel = cl::Kernel(program, "computeCMAPTorsionForces");
+}
+
+void OpenCLCalcCMAPTorsionForceKernel::executeForces(ContextImpl& context) {
+    if (numTorsions == 0)
+        return;
+    if (!hasInitializedKernel) {
+        hasInitializedKernel = true;
+        kernel.setArg<cl_int>(0, cl.getPaddedNumAtoms());
+        kernel.setArg<cl_int>(1, numTorsions);
+        kernel.setArg<cl::Buffer>(2, cl.getForceBuffers().getDeviceBuffer());
+        kernel.setArg<cl::Buffer>(3, cl.getEnergyBuffer().getDeviceBuffer());
+        kernel.setArg<cl::Buffer>(4, cl.getPosq().getDeviceBuffer());
+        kernel.setArg<cl::Buffer>(5, coefficients->getDeviceBuffer());
+        kernel.setArg<cl::Buffer>(6, mapPositions->getDeviceBuffer());
+        kernel.setArg<cl::Buffer>(7, torsionIndices->getDeviceBuffer());
+        kernel.setArg<cl::Buffer>(8, torsionMaps->getDeviceBuffer());
+    }
+    cl.executeKernel(kernel, numTorsions);
+}
+
+double OpenCLCalcCMAPTorsionForceKernel::executeEnergy(ContextImpl& context) {
+    executeForces(context);
+    return 0.0;
+}
+
 class OpenCLCustomTorsionForceInfo : public OpenCLForceInfo {
 public:
     OpenCLCustomTorsionForceInfo(int requiredBuffers, const CustomTorsionForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {

platforms/opencl/src/OpenCLKernels.h

@@ -356,7 +356,7 @@ private:
 class OpenCLCalcPeriodicTorsionForceKernel : public CalcPeriodicTorsionForceKernel {
 public:
     OpenCLCalcPeriodicTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcPeriodicTorsionForceKernel(name, platform),
-            hasInitializedKernel(false), cl(cl), system(system) {
+            hasInitializedKernel(false), cl(cl), system(system), params(NULL), indices(NULL) {
     }
     ~OpenCLCalcPeriodicTorsionForceKernel();
     /**
@@ -395,7 +395,7 @@ private:
 class OpenCLCalcRBTorsionForceKernel : public CalcRBTorsionForceKernel {
 public:
     OpenCLCalcRBTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcRBTorsionForceKernel(name, platform),
-            hasInitializedKernel(false), cl(cl), system(system) {
+            hasInitializedKernel(false), cl(cl), system(system), params(NULL), indices(NULL) {
     }
     ~OpenCLCalcRBTorsionForceKernel();
     /**
@@ -428,6 +428,47 @@ private:
     cl::Kernel kernel;
 };
 
+/**
+ * This kernel is invoked by CMAPTorsionForce to calculate the forces acting on the system and the energy of the system.
+ */
+class OpenCLCalcCMAPTorsionForceKernel : public CalcCMAPTorsionForceKernel {
+public:
+    OpenCLCalcCMAPTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCMAPTorsionForceKernel(name, platform),
+            hasInitializedKernel(false), cl(cl), system(system), coefficients(NULL), mapPositions(NULL), torsionIndices(NULL), torsionMaps(NULL) {
+    }
+    ~OpenCLCalcCMAPTorsionForceKernel();
+    /**
+     * 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.
+     *
+     * @param context    the context in which to execute this kernel
+     */
+    void executeForces(ContextImpl& context);
+    /**
+     * Execute the kernel to calculate the energy.
+     *
+     * @param context    the context in which to execute this kernel
+     * @return the potential energy due to the CMAPTorsionForce
+     */
+    double executeEnergy(ContextImpl& context);
+private:
+    int numTorsions;
+    bool hasInitializedKernel;
+    OpenCLContext& cl;
+    System& system;
+    OpenCLArray<mm_float4>* coefficients;
+    OpenCLArray<mm_int2>* mapPositions;
+    OpenCLArray<mm_int16>* torsionIndices;
+    OpenCLArray<cl_int>* torsionMaps;
+    cl::Kernel kernel;
+};
+
 /**
  * This kernel is invoked by CustomTorsionForce to calculate the forces acting on the system and the energy of the system.
  */

platforms/opencl/src/OpenCLPlatform.cpp

@@ -53,6 +53,7 @@ OpenCLPlatform::OpenCLPlatform() {
     registerKernelFactory(CalcCustomAngleForceKernel::Name(), factory);
     registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory);
+    registerKernelFactory(CalcCMAPTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomNonbondedForceKernel::Name(), factory);

platforms/opencl/src/kernels/cmapTorsionForce.cl

+/**
+ * Compute CNAP torsion forces.
+ */
+
+__kernel void computeCMAPTorsionForces(int numAtoms, int numTorsions, __global float4* forceBuffers, __global float* energyBuffer,
+        __global float4* posq, __global float4* coeff, __global int2* mapPositions, __global int16* indices, __global int* maps) {
+    const float PI = 3.14159265358979323846f;
+    float energy = 0.0f;
+    for (int index = get_global_id(0); index < numTorsions; index += get_global_size(0)) {
+        int16 atoms = indices[index];
+        float4 a1 = posq[atoms.s0];
+        float4 a2 = posq[atoms.s1];
+        float4 a3 = posq[atoms.s2];
+        float4 a4 = posq[atoms.s3];
+        float4 b1 = posq[atoms.s4];
+        float4 b2 = posq[atoms.s5];
+        float4 b3 = posq[atoms.s6];
+        float4 b4 = posq[atoms.s7];
+
+        // Compute the first angle.
+
+        float4 v0a = (float4) (a1.xyz-a2.xyz, 0.0f);
+        float4 v1a = (float4) (a3.xyz-a2.xyz, 0.0f);
+        float4 v2a = (float4) (a3.xyz-a4.xyz, 0.0f);
+        float4 cp0a = cross(v0a, v1a);
+        float4 cp1a = cross(v1a, v2a);
+        float cosangle = dot(normalize(cp0a), normalize(cp1a));
+        float 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.
+
+            float4 cross_prod = cross(cp0a, cp1a);
+            float 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.
+
+        float4 v0b = (float4) (b1.xyz-b2.xyz, 0.0f);
+        float4 v1b = (float4) (b3.xyz-b2.xyz, 0.0f);
+        float4 v2b = (float4) (b3.xyz-b4.xyz, 0.0f);
+        float4 cp0b = cross(v0b, v1b);
+        float4 cp1b = cross(v1b, v2b);
+        cosangle = dot(normalize(cp0b), normalize(cp1b));
+        float 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.
+
+            float4 cross_prod = cross(cp0b, cp1b);
+            float 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 = mapPositions[maps[index]];
+        int size = pos.y;
+        float delta = 2*PI/size;
+        int s = (int) (angleA/delta);
+        int t = (int) (angleB/delta);
+        float4 c[4];
+        int coeffIndex = 4*(pos.x+s+size*t);
+        c[0] = coeff[coeffIndex];
+        c[1] = coeff[coeffIndex+1];
+        c[2] = coeff[coeffIndex+2];
+        c[3] = coeff[coeffIndex+3];
+        float da = angleA/delta-s;
+        float db = angleB/delta-t;
+
+        // Evaluate the spline to determine the energy and gradients.
+
+        float torsionEnergy = 0.0f;
+        float dEdA = 0.0f;
+        float 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.
+
+        float normCross1 = dot(cp0a, cp0a);
+        float normSqrBC = dot(v1a, v1a);
+        float normBC = sqrt(normSqrBC);
+        float normCross2 = dot(cp1a, cp1a);
+        float dp = 1.0f/normSqrBC;
+        float4 ff = (float4) ((-dEdA*normBC)/normCross1, dot(v0a, v1a)*dp, dot(v2a, v1a)*dp, (dEdA*normBC)/normCross2);
+        float4 internalF0 = ff.x*cp0a;
+        float4 internalF3 = ff.w*cp1a;
+        float4 d = ff.y*internalF0 - ff.z*internalF3;
+        int4 offset = atoms.lo.lo + numAtoms*atoms.hi.lo;
+        float4 forceA = forceBuffers[offset.x];
+        float4 forceB = forceBuffers[offset.y];
+        float4 forceC = forceBuffers[offset.z];
+        float4 forceD = forceBuffers[offset.w];
+        forceA.xyz += internalF0.xyz;
+        forceB.xyz += d.xyz-internalF0.xyz;
+        forceC.xyz += -d.xyz-internalF3.xyz;
+        forceD.xyz += internalF3.xyz;
+        forceBuffers[offset.x] = forceA;
+        forceBuffers[offset.y] = forceB;
+        forceBuffers[offset.z] = forceC;
+        forceBuffers[offset.w] = forceD;
+
+        // Apply the force to the second torsion.
+
+        normCross1 = dot(cp0b, cp0b);
+        normSqrBC = dot(v1b, v1b);
+        normBC = sqrt(normSqrBC);
+        normCross2 = dot(cp1b, cp1b);
+        dp = 1.0f/normSqrBC;
+        ff = (float4) ((-dEdB*normBC)/normCross1, dot(v0b, v1b)*dp, dot(v2b, v1b)*dp, (dEdB*normBC)/normCross2);
+        internalF0 = ff.x*cp0b;
+        internalF3 = ff.w*cp1b;
+        d = ff.y*internalF0 - ff.z*internalF3;
+        offset = atoms.lo.hi + numAtoms*atoms.hi.hi;
+        forceA = forceBuffers[offset.x];
+        forceB = forceBuffers[offset.y];
+        forceC = forceBuffers[offset.z];
+        forceD = forceBuffers[offset.w];
+        forceA.xyz += internalF0.xyz;
+        forceB.xyz += d.xyz-internalF0.xyz;
+        forceC.xyz += -d.xyz-internalF3.xyz;
+        forceD.xyz += internalF3.xyz;
+        forceBuffers[offset.x] = forceA;
+        forceBuffers[offset.y] = forceB;
+        forceBuffers[offset.z] = forceC;
+        forceBuffers[offset.w] = forceD;
+    }
+    energyBuffer[get_global_id(0)] += energy;
+}

platforms/opencl/tests/TestOpenCLCMAPTorsionForce.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 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 OpenCL implementation of CMAPTorsionForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "OpenCLPlatform.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.
+
+    OpenCLPlatform 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/reference/tests/TestReferenceCMAPTorsionForce.cpp

@@ -30,7 +30,7 @@
  * -------------------------------------------------------------------------- */
 
 /**
- * This tests all the different force terms in the reference implementation of CMAPTorsionForce.
+ * This tests the reference implementation of CMAPTorsionForce.
  */
 
 #include "../../../tests/AssertionUtilities.h"