Created CUDA version of CustomAngleForce

platforms/cuda/src/CudaKernelFactory.cpp

@@ -43,6 +43,8 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
         return new CudaCalcCustomBondForceKernel(name, platform, data, context.getSystem());
     if (name == CalcHarmonicAngleForceKernel::Name())
         return new CudaCalcHarmonicAngleForceKernel(name, platform, data, context.getSystem());
+    if (name == CalcCustomAngleForceKernel::Name())
+        return new CudaCalcCustomAngleForceKernel(name, platform, data, context.getSystem());
     if (name == CalcPeriodicTorsionForceKernel::Name())
         return new CudaCalcPeriodicTorsionForceKernel(name, platform, data, context.getSystem());
     if (name == CalcRBTorsionForceKernel::Name())

platforms/cuda/src/CudaKernels.cpp

@@ -278,6 +278,54 @@ double CudaCalcHarmonicAngleForceKernel::executeEnergy(ContextImpl& context) {
     return 0.0;
 }
 
+CudaCalcCustomAngleForceKernel::~CudaCalcCustomAngleForceKernel() {
+}
+
+void CudaCalcCustomAngleForceKernel::initialize(const System& system, const CustomAngleForce& force) {
+    numAngles = force.getNumAngles();
+    vector<int> particle1(numAngles);
+    vector<int> particle2(numAngles);
+    vector<int> particle3(numAngles);
+    vector<vector<double> > params(numAngles);
+    for (int i = 0; i < numAngles; i++)
+        force.getAngleParameters(i, particle1[i], particle2[i], particle3[i], params[i]);
+    vector<string> paramNames;
+    for (int i = 0; i < force.getNumPerAngleParameters(); i++)
+        paramNames.push_back(force.getPerAngleParameterName(i));
+    globalParamNames.resize(force.getNumGlobalParameters());
+    globalParamValues.resize(force.getNumGlobalParameters());
+    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
+        globalParamNames[i] = force.getGlobalParameterName(i);
+        globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
+    }
+    gpuSetCustomAngleParameters(data.gpu, particle1, particle2, particle3, params, force.getEnergyFunction(), paramNames, globalParamNames);
+    if (globalParamValues.size() > 0)
+        SetCustomAngleGlobalParams(globalParamValues);
+}
+
+void CudaCalcCustomAngleForceKernel::executeForces(ContextImpl& context) {
+    updateGlobalParams(context);
+    kCalculateCustomAngleForces(data.gpu);
+}
+
+double CudaCalcCustomAngleForceKernel::executeEnergy(ContextImpl& context) {
+    updateGlobalParams(context);
+    kCalculateCustomAngleForces(data.gpu);
+    return 0.0;
+}
+
+void CudaCalcCustomAngleForceKernel::updateGlobalParams(ContextImpl& context) {
+    bool changed = false;
+    for (int i = 0; i < (int) globalParamNames.size(); i++) {
+        float value = (float) context.getParameter(globalParamNames[i]);
+        if (value != globalParamValues[i])
+            changed = true;
+        globalParamValues[i] = value;
+    }
+    if (changed)
+        SetCustomAngleGlobalParams(globalParamValues);
+}
+
 CudaCalcPeriodicTorsionForceKernel::~CudaCalcPeriodicTorsionForceKernel() {
 }
 

platforms/cuda/src/CudaKernels.h

@@ -256,6 +256,44 @@ private:
     System& system;
 };
 
+/**
+ * This kernel is invoked by CustomAngleForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CudaCalcCustomAngleForceKernel : public CalcCustomAngleForceKernel {
+public:
+    CudaCalcCustomAngleForceKernel(std::string name, const Platform& platform, CudaPlatform::PlatformData& data, System& system) : CalcCustomAngleForceKernel(name, platform),
+            data(data), system(system) {
+    }
+    ~CudaCalcCustomAngleForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomAngleForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomAngleForce& 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 CustomAngleForce
+     */
+    double executeEnergy(ContextImpl& context);
+private:
+    void updateGlobalParams(ContextImpl& context);
+    int numAngles;
+    CudaPlatform::PlatformData& data;
+    std::vector<std::string> globalParamNames;
+    std::vector<float> globalParamValues;
+    System& system;
+};
+
 /**
  * This kernel is invoked by PeriodicTorsionForce to calculate the forces acting on the system and the energy of the system.
  */

platforms/cuda/src/CudaPlatform.cpp

@@ -50,6 +50,7 @@ CudaPlatform::CudaPlatform() {
     registerKernelFactory(CalcHarmonicBondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomBondForceKernel::Name(), factory);
     registerKernelFactory(CalcHarmonicAngleForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomAngleForceKernel::Name(), factory);
     registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);

platforms/cuda/src/kernels/cudaKernels.h

@@ -43,6 +43,7 @@ extern void kCalculateCDLJObcGbsaForces1(gpuContext gpu);
 extern void kCalculateCDLJGBVIForces1(gpuContext gpu);
 extern void kCalculateCDLJForces(gpuContext gpu);
 extern void kCalculateCustomBondForces(gpuContext gpu);
+extern void kCalculateCustomAngleForces(gpuContext gpu);
 extern void kCalculateCustomExternalForces(gpuContext gpu);
 extern void kCalculateCustomNonbondedForces(gpuContext gpu, bool neighborListValid);
 extern void kReduceObcGbsaBornForces(gpuContext gpu);
@@ -77,6 +78,8 @@ extern void SetCalculateCDLJForcesSim(gpuContext gpu);
 extern void GetCalculateCDLJForcesSim(gpuContext gpu);
 extern void SetCalculateCustomBondForcesSim(gpuContext gpu);
 extern void GetCalculateCustomBondForcesSim(gpuContext gpu);
+extern void SetCalculateCustomAngleForcesSim(gpuContext gpu);
+extern void GetCalculateCustomAngleForcesSim(gpuContext gpu);
 extern void SetCalculateCustomExternalForcesSim(gpuContext gpu);
 extern void GetCalculateCustomExternalForcesSim(gpuContext gpu);
 extern void SetCalculateCustomNonbondedForcesSim(gpuContext gpu);
@@ -114,6 +117,9 @@ extern void GetRandomSim(gpuContext gpu);
 extern void SetCustomBondForceExpression(const Expression<256>& expression);
 extern void SetCustomBondEnergyExpression(const Expression<256>& expression);
 extern void SetCustomBondGlobalParams(const std::vector<float>&  paramValues);
+extern void SetCustomAngleForceExpression(const Expression<256>& expression);
+extern void SetCustomAngleEnergyExpression(const Expression<256>& expression);
+extern void SetCustomAngleGlobalParams(const std::vector<float>&  paramValues);
 extern void SetCustomExternalForceExpressions(const Expression<256>& expressionX, const Expression<256>& expressionY, const Expression<256>& expressionZ);
 extern void SetCustomExternalEnergyExpression(const Expression<256>& expression);
 extern void SetCustomExternalGlobalParams(const std::vector<float>& paramValues);

platforms/cuda/src/kernels/cudatypes.h

@@ -352,14 +352,17 @@ struct cudaGmxSimulation {
     float4*         pGBVIData;                      // Pointer to fixed Born data for GB/VI algorithm
     float2*         pAttr;                          // Pointer to additional atom attributes (sig, eps)
     float4*         pCustomParams;                  // Pointer to atom parameters for custom nonbonded force
-    int4*           pCustomExceptionID;             // Atom indices for custom nonbonded exceptions
-    float4*         pCustomExceptionParams;         // Parameters for custom nonbonded exceptions
     unsigned int    customExceptions;               // Number of custom nonbonded exceptions
     unsigned int    customParameters;               // Number of parameters for custom nonbonded interactions
     int4*           pCustomBondID;                  // Atom indices for custom bonds
     float4*         pCustomBondParams;              // Parameters for custom bonds
     unsigned int    customBonds;                    // Number of custom bonds
     unsigned int    customBondParameters;           // Number of parameters for custom bonds
+    int4*           pCustomAngleID1;                // Atom indices for custom angles
+    int2*           pCustomAngleID2;                // Atom indices for custom angles
+    float4*         pCustomAngleParams;             // Parameters for custom angles
+    unsigned int    customAngles;                   // Number of custom angles
+    unsigned int    customAngleParameters;          // Number of parameters for custom angles
     int*            pCustomExternalID;              // Atom indices for custom external force
     float4*         pCustomExternalParams;          // Parameters for custom external force
     unsigned int    customExternals;                // Number of particles for custom external force

platforms/cuda/src/kernels/gpu.cpp

@@ -701,6 +701,58 @@ void gpuSetCustomBondParameters(gpuContext gpu, const vector<int>& bondAtom1, co
     SetCustomBondForceExpression(createExpression<256>(gpu, energyExp, Lepton::Parser::parse(energyExp).differentiate("r").optimize().createProgram(), variables, globalParamNames, gpu->sim.customExpressionStackSize));
 }
 
+extern "C"
+void gpuSetCustomAngleParameters(gpuContext gpu, const vector<int>& angleAtom1, const vector<int>& angleAtom2, const vector<int>& angleAtom3, const vector<vector<double> >& angleParams,
+            const string& energyExp, const vector<string>& paramNames, const vector<string>& globalParamNames)
+{
+    if (paramNames.size() > 4)
+        throw OpenMMException("CudaPlatform only supports four per-angle parameters for custom angle forces");
+    if (globalParamNames.size() > 8)
+        throw OpenMMException("CudaPlatform only supports eight global parameters for custom angle forces");
+    if (gpu->psCustomAngleID1 != NULL)
+        throw OpenMMException("CudaPlatform only supports a single CustomAngleForce per System");
+    gpu->sim.customAngles = angleAtom1.size();
+    gpu->sim.customAngleParameters = paramNames.size();
+    gpu->psCustomAngleID1 = new CUDAStream<int4>(gpu->sim.customAngles, 1, "CustomAngleId1");
+    gpu->sim.pCustomAngleID1 = gpu->psCustomAngleID1->_pDevData;
+    gpu->psCustomAngleID2 = new CUDAStream<int2>(gpu->sim.customAngles, 1, "CustomAngleId2");
+    gpu->sim.pCustomAngleID2 = gpu->psCustomAngleID2->_pDevData;
+    gpu->psCustomAngleParams = new CUDAStream<float4>(gpu->sim.customAngles, 1, "CustomAngleParams");
+    gpu->sim.pCustomAngleParams = gpu->psCustomAngleParams->_pDevData;
+    vector<int> forceBufferCounter(gpu->natoms, 0);
+    for (int i = 0; i < (int) angleAtom1.size(); i++) {
+        (*gpu->psCustomAngleID1)[i].x = angleAtom1[i];
+        (*gpu->psCustomAngleID1)[i].y = angleAtom2[i];
+        (*gpu->psCustomAngleID1)[i].z = angleAtom3[i];
+        (*gpu->psCustomAngleID1)[i].w = forceBufferCounter[angleAtom1[i]]++;
+        (*gpu->psCustomAngleID2)[i].x = forceBufferCounter[angleAtom2[i]]++;
+        (*gpu->psCustomAngleID2)[i].y = forceBufferCounter[angleAtom2[i]]++;
+        if (angleParams[i].size() > 0)
+            (*gpu->psCustomAngleParams)[i].x = (float) angleParams[i][0];
+        if (angleParams[i].size() > 1)
+            (*gpu->psCustomAngleParams)[i].y = (float) angleParams[i][1];
+        if (angleParams[i].size() > 2)
+            (*gpu->psCustomAngleParams)[i].z = (float) angleParams[i][2];
+        if (angleParams[i].size() > 3)
+            (*gpu->psCustomAngleParams)[i].w = (float) angleParams[i][3];
+    }
+    gpu->psCustomAngleID1->Upload();
+    gpu->psCustomAngleID2->Upload();
+    gpu->psCustomAngleParams->Upload();
+    for (int i = 0; i < (int) forceBufferCounter.size(); i++)
+        if (forceBufferCounter[i] > (int) gpu->pOutputBufferCounter[i])
+            gpu->pOutputBufferCounter[i] = forceBufferCounter[i];
+
+    // Create the Expressions.
+
+    vector<string> variables;
+    variables.push_back("theta");
+    for (int i = 0; i < (int) paramNames.size(); i++)
+        variables.push_back(paramNames[i]);
+    SetCustomAngleEnergyExpression(createExpression<256>(gpu, energyExp, Lepton::Parser::parse(energyExp).optimize().createProgram(), variables, globalParamNames, gpu->sim.customExpressionStackSize));
+    SetCustomAngleForceExpression(createExpression<256>(gpu, energyExp, Lepton::Parser::parse(energyExp).differentiate("theta").optimize().createProgram(), variables, globalParamNames, gpu->sim.customExpressionStackSize));
+}
+
 extern "C"
 void gpuSetCustomExternalParameters(gpuContext gpu, const vector<int>& atomIndex, const vector<vector<double> >& atomParams,
             const string& energyExp, const vector<string>& paramNames, const vector<string>& globalParamNames)
@@ -1829,10 +1881,11 @@ void* gpuInit(int numAtoms, unsigned int device, bool useBlockingSync)
     gpu->psLJ14ID                   = NULL;
     gpu->psLJ14Parameter            = NULL;
     gpu->psCustomParams             = NULL;
-    gpu->psCustomExceptionID        = NULL;
-    gpu->psCustomExceptionParams    = NULL;
     gpu->psCustomBondID             = NULL;
     gpu->psCustomBondParams         = NULL;
+    gpu->psCustomAngleID1           = NULL;
+    gpu->psCustomAngleID2           = NULL;
+    gpu->psCustomAngleParams        = NULL;
     gpu->psCustomExternalID         = NULL;
     gpu->psCustomExternalParams     = NULL;
     gpu->psEwaldCosSinSum           = NULL;
@@ -1999,15 +2052,17 @@ void gpuShutDown(gpuContext gpu)
     delete gpu->psxVector4;
     delete gpu->psvVector4;
     delete gpu->psSigEps2;
-    if (gpu->psCustomParams != NULL) {
+    if (gpu->psCustomParams != NULL)
         delete gpu->psCustomParams;
-        delete gpu->psCustomExceptionID;
-        delete gpu->psCustomExceptionParams;
-    }
     if (gpu->psCustomBondParams != NULL) {
         delete gpu->psCustomBondID;
         delete gpu->psCustomBondParams;
     }
+    if (gpu->psCustomAngleParams != NULL) {
+        delete gpu->psCustomAngleID1;
+        delete gpu->psCustomAngleID2;
+        delete gpu->psCustomAngleParams;
+    }
     if (gpu->psCustomExternalParams != NULL) {
         delete gpu->psCustomExternalID;
         delete gpu->psCustomExternalParams;
@@ -2380,6 +2435,7 @@ int gpuSetConstants(gpuContext gpu)
     SetCalculateCDLJObcGbsaForces1Sim(gpu);
     SetCalculateCustomNonbondedForcesSim(gpu);
     SetCalculateCustomBondForcesSim(gpu);
+    SetCalculateCustomAngleForcesSim(gpu);
     SetCalculateCustomExternalForcesSim(gpu);
     SetCalculateLocalForcesSim(gpu);
     SetCalculateObcGbsaBornSumSim(gpu);

platforms/cuda/src/kernels/gputypes.h

@@ -102,10 +102,11 @@ struct _gpuContext {
     CUDAStream<float4>* psvVector4;
     CUDAStream<float2>* psSigEps2; 
     CUDAStream<float4>* psCustomParams;     // Atom parameters for custom nonbonded force
-    CUDAStream<int4>* psCustomExceptionID;  // Atom indices for custom nonbonded exceptions
-    CUDAStream<float4>* psCustomExceptionParams; // Parameters for custom nonbonded exceptions
     CUDAStream<int4>* psCustomBondID;             // Atom indices for custom bonds
     CUDAStream<float4>* psCustomBondParams;       // Parameters for custom bonds
+    CUDAStream<int4>* psCustomAngleID1;           // Atom indices for custom angles
+    CUDAStream<int2>* psCustomAngleID2;           // Atom indices for custom angles
+    CUDAStream<float4>* psCustomAngleParams;      // Parameters for custom angles
     CUDAStream<int>* psCustomExternalID;          // Atom indices for custom external force
     CUDAStream<float4>* psCustomExternalParams;   // Parameters for custom external force
     CUDAStream<float4>* psTabulatedFunctionParams; // The min, max, and spacing for each tabulated function
@@ -214,6 +215,10 @@ extern "C"
 void gpuSetCustomBondParameters(gpuContext gpu, const std::vector<int>& bondAtom1, const std::vector<int>& bondAtom2, const std::vector<std::vector<double> >& bondParams,
             const std::string& energyExp, const std::vector<std::string>& paramNames, const std::vector<std::string>& globalParamNames);
 
+extern "C"
+void gpuSetCustomAngleParameters(gpuContext gpu, const std::vector<int>& bondAtom1, const std::vector<int>& bondAtom2, const std::vector<int>& bondAtom3, const std::vector<std::vector<double> >& bondParams,
+            const std::string& energyExp, const std::vector<std::string>& paramNames, const std::vector<std::string>& globalParamNames);
+
 extern "C"
 void gpuSetCustomExternalParameters(gpuContext gpu, const std::vector<int>& atomIndex, const std::vector<std::vector<double> >& atomParams,
             const std::string& energyExp, const std::vector<std::string>& paramNames, const std::vector<std::string>& globalParamNames);

platforms/cuda/src/kernels/kCalculateCustomAngleForces.cu

+/* -------------------------------------------------------------------------- *
+ *                                   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: Scott Le Grand, Peter Eastman                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include <stdio.h>
+#include <cuda.h>
+#include <vector_functions.h>
+#include <cstdlib>
+#include <string>
+#include <iostream>
+#include <fstream>
+using namespace std;
+
+#include "gputypes.h"
+#include "cudatypes.h"
+
+static __constant__ cudaGmxSimulation cSim;
+static __constant__ Expression<256> forceExp;
+static __constant__ Expression<256> energyExp;
+
+#include "kEvaluateExpression.h"
+
+void SetCalculateCustomAngleForcesSim(gpuContext gpu)
+{
+    cudaError_t status;
+    status = cudaMemcpyToSymbol(cSim, &gpu->sim, sizeof(cudaGmxSimulation));
+    RTERROR(status, "cudaMemcpyToSymbol: SetSim copy to cSim failed");
+}
+
+void GetCalculateCustomAngleForcesSim(gpuContext gpu)
+{
+    cudaError_t status;
+    status = cudaMemcpyFromSymbol(&gpu->sim, cSim, sizeof(cudaGmxSimulation));
+    RTERROR(status, "cudaMemcpyFromSymbol: SetSim copy from cSim failed");
+}
+
+void SetCustomAngleForceExpression(const Expression<256>& expression)
+{
+    cudaError_t status;
+    status = cudaMemcpyToSymbol(forceExp, &expression, sizeof(forceExp));
+    RTERROR(status, "SetCustomAngleForceExpression: cudaMemcpyToSymbol failed");
+}
+
+void SetCustomAngleEnergyExpression(const Expression<256>& expression)
+{
+    cudaError_t status;
+    status = cudaMemcpyToSymbol(energyExp, &expression, sizeof(energyExp));
+    RTERROR(status, "SetCustomAngleEnergyExpression: cudaMemcpyToSymbol failed");
+}
+
+void SetCustomAngleGlobalParams(const vector<float>& paramValues)
+{
+    cudaError_t status;
+    status = cudaMemcpyToSymbol(globalParams, &paramValues[0], paramValues.size()*sizeof(float));
+    RTERROR(status, "SetCustomAngleGlobalParams: cudaMemcpyToSymbol failed");
+}
+
+#define DOT3(v1, v2) (v1.x*v2.x + v1.y*v2.y + v1.z*v2.z)
+
+#define CROSS_PRODUCT(v1, v2) make_float3(v1.y*v2.z - v1.z*v2.y, v1.z*v2.x - v1.x*v2.z, v1.x*v2.y - v1.y*v2.x)
+
+__global__ void kCalculateCustomAngleForces_kernel()
+{
+    extern __shared__ float stack[];
+    float* variables = (float*) &stack[cSim.customExpressionStackSize*blockDim.x];
+    unsigned int pos = blockIdx.x * blockDim.x + threadIdx.x;
+    float totalEnergy = 0.0f;
+
+    while (pos < cSim.customAngles)
+    {
+        int4 atom       = cSim.pCustomAngleID1[pos];
+        int2 atom2      = cSim.pCustomAngleID2[pos];
+        float4 params   = cSim.pCustomAngleParams[pos];
+        float4 a1       = cSim.pPosq[atom.x];
+        float4 a2       = cSim.pPosq[atom.y];
+        float4 a3       = cSim.pPosq[atom.z];
+        float3 v0 = make_float3(a2.x-a1.x, a2.y-a1.y, a2.z-a1.z);
+        float3 v1 = make_float3(a2.x-a3.x, a2.y-a3.y, a2.z-a3.z);
+        float3 cp = CROSS_PRODUCT(v0, v1);
+        float rp = DOT3(cp, cp);
+        rp = max(sqrt(rp), 1.0e-06f);
+        float r21 = DOT3(v0, v0);
+        float r23 = DOT3(v1, v1);
+        float dot = DOT3(v0, v1);
+        float cosine = dot/sqrt(r21*r23);
+        VARIABLE(0) = acos(cosine);
+        VARIABLE(1) = params.x;
+        VARIABLE(2) = params.y;
+        VARIABLE(3) = params.z;
+        VARIABLE(4) = params.w;
+        float dEdR = kEvaluateExpression_kernel(&forceExp, stack, variables);
+        totalEnergy += kEvaluateExpression_kernel(&energyExp, stack, variables);
+        float termA =  dEdR/(r21*rp);
+        float termC = -dEdR/(r23*rp);
+        float3 c21 = CROSS_PRODUCT(v0, cp);
+        float3 c23 = CROSS_PRODUCT(v1, cp);
+        c21.x *= termA;
+        c21.y *= termA;
+        c21.z *= termA;
+        c23.x *= termC;
+        c23.y *= termC;
+        c23.z *= termC;
+        unsigned int offsetA  = atom.x + atom.w * cSim.stride;
+        unsigned int offsetB  = atom.y + atom2.x * cSim.stride;
+        unsigned int offsetC  = atom.z + atom2.y * cSim.stride;
+        float4 forceA         = cSim.pForce4[offsetA];
+        float4 forceB         = cSim.pForce4[offsetB];
+        float4 forceC         = cSim.pForce4[offsetC];
+        forceA.x             += c21.x;
+        forceA.y             += c21.y;
+        forceA.z             += c21.z;
+        forceB.x             -= c21.x+c23.x;
+        forceB.y             -= c21.y+c23.y;
+        forceB.z             -= c21.z+c23.z;
+        forceC.x             += c23.x;
+        forceC.y             += c23.y;
+        forceC.z             += c23.z;
+        cSim.pForce4[offsetA] = forceA;
+        cSim.pForce4[offsetB] = forceB;
+        cSim.pForce4[offsetC] = forceC;
+        pos                  += blockDim.x * gridDim.x;
+    }
+    cSim.pEnergy[blockIdx.x * blockDim.x + threadIdx.x] += totalEnergy;
+}
+
+void kCalculateCustomAngleForces(gpuContext gpu)
+{
+//    printf("kCalculateCustomAngleForces\n");
+    kCalculateCustomAngleForces_kernel<<<gpu->sim.blocks, gpu->sim.localForces_threads_per_block,
+            (gpu->sim.customExpressionStackSize+9)*sizeof(float)*gpu->sim.localForces_threads_per_block>>>();
+    LAUNCHERROR("kCalculateCustomAngleForces");
+}

platforms/cuda/tests/TestCudaCustomAngleForce.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-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 CUDA implementation of CustomAngleForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "CudaPlatform.h"
+#include "openmm/CustomAngleForce.h"
+#include "openmm/HarmonicAngleForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "../src/sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testAngles() {
+    CudaPlatform platform;
+
+    // Create a system using a CustomAngleForce.
+
+    System customSystem;
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    CustomAngleForce* custom = new CustomAngleForce("scale*k*(theta-theta0)^2");
+    custom->addPerAngleParameter("theta0");
+    custom->addPerAngleParameter("k");
+    custom->addGlobalParameter("scale", 0.5);
+    vector<double> parameters(2);
+    parameters[0] = 1.5;
+    parameters[1] = 0.8;
+    custom->addAngle(0, 1, 2, parameters);
+    parameters[0] = 2.0;
+    parameters[1] = 0.5;
+    custom->addAngle(1, 2, 3, parameters);
+    customSystem.addForce(custom);
+
+    // Create an identical system using a HarmonicAngleForce.
+
+    System harmonicSystem;
+    harmonicSystem.addParticle(1.0);
+    harmonicSystem.addParticle(1.0);
+    harmonicSystem.addParticle(1.0);
+    harmonicSystem.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    HarmonicAngleForce* harmonic = new HarmonicAngleForce();
+    harmonic->addAngle(0, 1, 2, 1.5, 0.8);
+    harmonic->addAngle(1, 2, 3, 2.0, 0.5);
+    harmonicSystem.addForce(harmonic);
+
+    // Set the atoms in various positions, and verify that both systems give identical forces and energy.
+
+    init_gen_rand(0);
+    vector<Vec3> positions(4);
+    VerletIntegrator integrator1(0.01);
+    VerletIntegrator integrator2(0.01);
+    for (int i = 0; i < 10; i++) {
+        for (int j = 0; j < (int) positions.size(); j++)
+            positions[j] = Vec3(5.0*genrand_real2(), 5.0*genrand_real2(), 5.0*genrand_real2());
+        double energy1, energy2;
+        vector<Vec3> forces1, forces2;
+        {
+            Context c(customSystem, integrator1, platform);
+            c.setPositions(positions);
+            State s = c.getState(State::Forces | State::Energy);
+            energy1 = s.getPotentialEnergy();
+            forces1 = s.getForces();
+        }
+        {
+            Context c(harmonicSystem, integrator1, platform);
+            c.setPositions(positions);
+            State s = c.getState(State::Forces | State::Energy);
+            energy2 = s.getPotentialEnergy();
+            forces2 = s.getForces();
+        }
+        ASSERT_EQUAL_VEC(forces2[0], forces1[0], TOL);
+        ASSERT_EQUAL_VEC(forces2[1], forces1[1], TOL);
+        ASSERT_EQUAL_VEC(forces2[2], forces1[2], TOL);
+        ASSERT_EQUAL_TOL(energy2, energy1, TOL);
+    }
+}
+
+int main() {
+    try {
+        testAngles();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+
+
+

platforms/reference/tests/TestReferenceCustomAngleForce.cpp

@@ -58,6 +58,7 @@ void testAngles() {
     customSystem.addParticle(1.0);
     customSystem.addParticle(1.0);
     customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
     CustomAngleForce* custom = new CustomAngleForce("scale*k*(theta-theta0)^2");
     custom->addPerAngleParameter("theta0");
     custom->addPerAngleParameter("k");
@@ -66,6 +67,9 @@ void testAngles() {
     parameters[0] = 1.5;
     parameters[1] = 0.8;
     custom->addAngle(0, 1, 2, parameters);
+    parameters[0] = 2.0;
+    parameters[1] = 0.5;
+    custom->addAngle(1, 2, 3, parameters);
     customSystem.addForce(custom);
 
     // Create an identical system using a HarmonicAngleForce.
@@ -74,23 +78,24 @@ void testAngles() {
     harmonicSystem.addParticle(1.0);
     harmonicSystem.addParticle(1.0);
     harmonicSystem.addParticle(1.0);
+    harmonicSystem.addParticle(1.0);
     VerletIntegrator integrator(0.01);
     HarmonicAngleForce* harmonic = new HarmonicAngleForce();
     harmonic->addAngle(0, 1, 2, 1.5, 0.8);
+    harmonic->addAngle(1, 2, 3, 2.0, 0.5);
     harmonicSystem.addForce(harmonic);
 
     // Set the atoms in various positions, and verify that both systems give identical forces and energy.
 
     init_gen_rand(0);
-    vector<Vec3> positions(3);
+    vector<Vec3> positions(4);
     VerletIntegrator integrator1(0.01);
     VerletIntegrator integrator2(0.01);
     for (int i = 0; i < 10; i++) {
         Context c1(customSystem, integrator1, platform);
         Context c2(harmonicSystem, integrator2, platform);
-        positions[0] = Vec3(5.0*genrand_real2(), 5.0*genrand_real2(), 5.0*genrand_real2());
-        positions[1] = Vec3(5.0*genrand_real2(), 5.0*genrand_real2(), 5.0*genrand_real2());
-        positions[2] = Vec3(5.0*genrand_real2(), 5.0*genrand_real2(), 5.0*genrand_real2());
+        for (int j = 0; j < (int) positions.size(); j++)
+            positions[j] = Vec3(5.0*genrand_real2(), 5.0*genrand_real2(), 5.0*genrand_real2());
         c1.setPositions(positions);
         c2.setPositions(positions);
         State s1 = c1.getState(State::Forces | State::Energy);