Created CUDA implementation of CustomTorsionForce

openmmapi/include/OpenMM.h

@@ -36,6 +36,8 @@
 #include "openmm/BrownianIntegrator.h"
 #include "openmm/CMMotionRemover.h"
 #include "openmm/CustomBondForce.h"
+#include "openmm/CustomAngleForce.h"
+#include "openmm/CustomTorsionForce.h"
 #include "openmm/CustomExternalForce.h"
 #include "openmm/CustomGBForce.h"
 #include "openmm/CustomNonbondedForce.h"

platforms/cuda/src/CudaKernelFactory.cpp

@@ -49,6 +49,8 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
         return new CudaCalcPeriodicTorsionForceKernel(name, platform, data, context.getSystem());
     if (name == CalcRBTorsionForceKernel::Name())
         return new CudaCalcRBTorsionForceKernel(name, platform, data, context.getSystem());
+    if (name == CalcCustomTorsionForceKernel::Name())
+        return new CudaCalcCustomTorsionForceKernel(name, platform, data, context.getSystem());
     if (name == CalcNonbondedForceKernel::Name())
         return new CudaCalcNonbondedForceKernel(name, platform, data, context.getSystem());
     if (name == CalcCustomNonbondedForceKernel::Name())

platforms/cuda/src/CudaKernels.cpp

@@ -392,6 +392,55 @@ double CudaCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
     return 0.0;
 }
 
+CudaCalcCustomTorsionForceKernel::~CudaCalcCustomTorsionForceKernel() {
+}
+
+void CudaCalcCustomTorsionForceKernel::initialize(const System& system, const CustomTorsionForce& force) {
+    numTorsions = force.getNumTorsions();
+    vector<int> particle1(numTorsions);
+    vector<int> particle2(numTorsions);
+    vector<int> particle3(numTorsions);
+    vector<int> particle4(numTorsions);
+    vector<vector<double> > params(numTorsions);
+    for (int i = 0; i < numTorsions; i++)
+        force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], params[i]);
+    vector<string> paramNames;
+    for (int i = 0; i < force.getNumPerTorsionParameters(); i++)
+        paramNames.push_back(force.getPerTorsionParameterName(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);
+    }
+    gpuSetCustomTorsionParameters(data.gpu, particle1, particle2, particle3, particle4, params, force.getEnergyFunction(), paramNames, globalParamNames);
+    if (globalParamValues.size() > 0)
+        SetCustomTorsionGlobalParams(globalParamValues);
+}
+
+void CudaCalcCustomTorsionForceKernel::executeForces(ContextImpl& context) {
+    updateGlobalParams(context);
+    kCalculateCustomTorsionForces(data.gpu);
+}
+
+double CudaCalcCustomTorsionForceKernel::executeEnergy(ContextImpl& context) {
+    updateGlobalParams(context);
+    kCalculateCustomTorsionForces(data.gpu);
+    return 0.0;
+}
+
+void CudaCalcCustomTorsionForceKernel::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)
+        SetCustomTorsionGlobalParams(globalParamValues);
+}
+
 CudaCalcNonbondedForceKernel::~CudaCalcNonbondedForceKernel() {
 }
 

platforms/cuda/src/CudaKernels.h

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

platforms/cuda/src/CudaPlatform.cpp

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

platforms/cuda/src/kernels/cudaKernels.h

@@ -44,6 +44,7 @@ extern void kCalculateCDLJGBVIForces1(gpuContext gpu);
 extern void kCalculateCDLJForces(gpuContext gpu);
 extern void kCalculateCustomBondForces(gpuContext gpu);
 extern void kCalculateCustomAngleForces(gpuContext gpu);
+extern void kCalculateCustomTorsionForces(gpuContext gpu);
 extern void kCalculateCustomExternalForces(gpuContext gpu);
 extern void kCalculateCustomNonbondedForces(gpuContext gpu, bool neighborListValid);
 extern void kReduceObcGbsaBornForces(gpuContext gpu);
@@ -80,6 +81,8 @@ extern void SetCalculateCustomBondForcesSim(gpuContext gpu);
 extern void GetCalculateCustomBondForcesSim(gpuContext gpu);
 extern void SetCalculateCustomAngleForcesSim(gpuContext gpu);
 extern void GetCalculateCustomAngleForcesSim(gpuContext gpu);
+extern void SetCalculateCustomTorsionForcesSim(gpuContext gpu);
+extern void GetCalculateCustomTorsionForcesSim(gpuContext gpu);
 extern void SetCalculateCustomExternalForcesSim(gpuContext gpu);
 extern void GetCalculateCustomExternalForcesSim(gpuContext gpu);
 extern void SetCalculateCustomNonbondedForcesSim(gpuContext gpu);
@@ -120,6 +123,9 @@ 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 SetCustomTorsionForceExpression(const Expression<256>& expression);
+extern void SetCustomTorsionEnergyExpression(const Expression<256>& expression);
+extern void SetCustomTorsionGlobalParams(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

@@ -363,6 +363,11 @@ struct cudaGmxSimulation {
     float4*         pCustomAngleParams;             // Parameters for custom angles
     unsigned int    customAngles;                   // Number of custom angles
     unsigned int    customAngleParameters;          // Number of parameters for custom angles
+    int4*           pCustomTorsionID1;              // Atom indices for custom torsions
+    int4*           pCustomTorsionID2;              // Atom indices for custom torsions
+    float4*         pCustomTorsionParams;           // Parameters for custom torsions
+    unsigned int    customTorsions;                 // Number of custom torsions
+    unsigned int    customTorsionParameters;        // Number of parameters for custom torsions
     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

@@ -112,6 +112,9 @@ struct Molecule {
     vector<int> rbTorsions;
     vector<int> constraints;
     vector<int> lj14s;
+    vector<int> customBonds;
+    vector<int> customAngles;
+    vector<int> customTorsions;
 };
 
 static const float dielectricOffset         =    0.009f;
@@ -753,6 +756,60 @@ void gpuSetCustomAngleParameters(gpuContext gpu, const vector<int>& angleAtom1,
     SetCustomAngleForceExpression(createExpression<256>(gpu, energyExp, Lepton::Parser::parse(energyExp).differentiate("theta").optimize().createProgram(), variables, globalParamNames, gpu->sim.customExpressionStackSize));
 }
 
+extern "C"
+void gpuSetCustomTorsionParameters(gpuContext gpu, const vector<int>& torsionAtom1, const vector<int>& torsionAtom2, const vector<int>& torsionAtom3, const vector<int>& torsionAtom4, const vector<vector<double> >& torsionParams,
+            const string& energyExp, const vector<string>& paramNames, const vector<string>& globalParamNames)
+{
+    if (paramNames.size() > 4)
+        throw OpenMMException("CudaPlatform only supports four per-torsion parameters for custom torsion forces");
+    if (globalParamNames.size() > 8)
+        throw OpenMMException("CudaPlatform only supports eight global parameters for custom torsion forces");
+    if (gpu->psCustomTorsionID1 != NULL)
+        throw OpenMMException("CudaPlatform only supports a single CustomTorsionForce per System");
+    gpu->sim.customTorsions = torsionAtom1.size();
+    gpu->sim.customTorsionParameters = paramNames.size();
+    gpu->psCustomTorsionID1 = new CUDAStream<int4>(gpu->sim.customTorsions, 1, "CustomTorsionId1");
+    gpu->sim.pCustomTorsionID1 = gpu->psCustomTorsionID1->_pDevData;
+    gpu->psCustomTorsionID2 = new CUDAStream<int4>(gpu->sim.customTorsions, 1, "CustomTorsionId2");
+    gpu->sim.pCustomTorsionID2 = gpu->psCustomTorsionID2->_pDevData;
+    gpu->psCustomTorsionParams = new CUDAStream<float4>(gpu->sim.customTorsions, 1, "CustomTorsionParams");
+    gpu->sim.pCustomTorsionParams = gpu->psCustomTorsionParams->_pDevData;
+    vector<int> forceBufferCounter(gpu->natoms, 0);
+    for (int i = 0; i < (int) torsionAtom1.size(); i++) {
+        (*gpu->psCustomTorsionID1)[i].x = torsionAtom1[i];
+        (*gpu->psCustomTorsionID1)[i].y = torsionAtom2[i];
+        (*gpu->psCustomTorsionID1)[i].z = torsionAtom3[i];
+        (*gpu->psCustomTorsionID1)[i].w = torsionAtom4[i];
+        (*gpu->psCustomTorsionID2)[i].x = forceBufferCounter[torsionAtom1[i]]++;
+        (*gpu->psCustomTorsionID2)[i].y = forceBufferCounter[torsionAtom2[i]]++;
+        (*gpu->psCustomTorsionID2)[i].z = forceBufferCounter[torsionAtom3[i]]++;
+        (*gpu->psCustomTorsionID2)[i].w = forceBufferCounter[torsionAtom4[i]]++;
+        if (torsionParams[i].size() > 0)
+            (*gpu->psCustomTorsionParams)[i].x = (float) torsionParams[i][0];
+        if (torsionParams[i].size() > 1)
+            (*gpu->psCustomTorsionParams)[i].y = (float) torsionParams[i][1];
+        if (torsionParams[i].size() > 2)
+            (*gpu->psCustomTorsionParams)[i].z = (float) torsionParams[i][2];
+        if (torsionParams[i].size() > 3)
+            (*gpu->psCustomTorsionParams)[i].w = (float) torsionParams[i][3];
+    }
+    gpu->psCustomTorsionID1->Upload();
+    gpu->psCustomTorsionID2->Upload();
+    gpu->psCustomTorsionParams->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]);
+    SetCustomTorsionEnergyExpression(createExpression<256>(gpu, energyExp, Lepton::Parser::parse(energyExp).optimize().createProgram(), variables, globalParamNames, gpu->sim.customExpressionStackSize));
+    SetCustomTorsionForceExpression(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)
@@ -1886,6 +1943,9 @@ void* gpuInit(int numAtoms, unsigned int device, bool useBlockingSync)
     gpu->psCustomAngleID1           = NULL;
     gpu->psCustomAngleID2           = NULL;
     gpu->psCustomAngleParams        = NULL;
+    gpu->psCustomTorsionID1         = NULL;
+    gpu->psCustomTorsionID2         = NULL;
+    gpu->psCustomTorsionParams      = NULL;
     gpu->psCustomExternalID         = NULL;
     gpu->psCustomExternalParams     = NULL;
     gpu->psEwaldCosSinSum           = NULL;
@@ -1927,6 +1987,8 @@ void* gpuInit(int numAtoms, unsigned int device, bool useBlockingSync)
         gpu->tabulatedFunctions[i].coefficients = NULL;
     gpu->sim.customExpressionStackSize = 0;
     gpu->sim.customBonds = 0;
+    gpu->sim.customAngles = 0;
+    gpu->sim.customTorsions = 0;
     
     // Initialize output buffer before reading parameters
     gpu->pOutputBufferCounter       = new unsigned int[gpu->sim.paddedNumberOfAtoms];
@@ -2063,6 +2125,11 @@ void gpuShutDown(gpuContext gpu)
         delete gpu->psCustomAngleID2;
         delete gpu->psCustomAngleParams;
     }
+    if (gpu->psCustomTorsionParams != NULL) {
+        delete gpu->psCustomTorsionID1;
+        delete gpu->psCustomTorsionID2;
+        delete gpu->psCustomTorsionParams;
+    }
     if (gpu->psCustomExternalParams != NULL) {
         delete gpu->psCustomExternalID;
         delete gpu->psCustomExternalParams;
@@ -2436,6 +2503,7 @@ int gpuSetConstants(gpuContext gpu)
     SetCalculateCustomNonbondedForcesSim(gpu);
     SetCalculateCustomBondForcesSim(gpu);
     SetCalculateCustomAngleForcesSim(gpu);
+    SetCalculateCustomTorsionForcesSim(gpu);
     SetCalculateCustomExternalForcesSim(gpu);
     SetCalculateLocalForcesSim(gpu);
     SetCalculateObcGbsaBornSumSim(gpu);
@@ -2569,6 +2637,21 @@ static void findMoleculeGroups(gpuContext gpu)
         int atom1 = (*gpu->psLJ14ID)[i].x;
         molecules[atomMolecule[atom1]].lj14s.push_back(i);
     }
+    for (int i = 0; i < (int)gpu->sim.customBonds; i++)
+    {
+        int atom1 = (*gpu->psCustomBondID)[i].x;
+        molecules[atomMolecule[atom1]].customBonds.push_back(i);
+    }
+    for (int i = 0; i < (int)gpu->sim.customAngles; i++)
+    {
+        int atom1 = (*gpu->psCustomAngleID1)[i].x;
+        molecules[atomMolecule[atom1]].customAngles.push_back(i);
+    }
+    for (int i = 0; i < (int)gpu->sim.customTorsions; i++)
+    {
+        int atom1 = (*gpu->psCustomTorsionID1)[i].x;
+        molecules[atomMolecule[atom1]].customTorsions.push_back(i);
+    }
 
     // Sort them into groups of identical molecules.
 
@@ -2588,7 +2671,8 @@ static void findMoleculeGroups(gpuContext gpu)
             if (mol.atoms.size() != mol2.atoms.size() || mol.bonds.size() != mol2.bonds.size()
                     || mol.angles.size() != mol2.angles.size() || mol.periodicTorsions.size() != mol2.periodicTorsions.size()
                     || mol.rbTorsions.size() != mol2.rbTorsions.size() || mol.constraints.size() != mol2.constraints.size()
-                    || mol.lj14s.size() != mol2.lj14s.size())
+                    || mol.lj14s.size() != mol2.lj14s.size() || mol.customBonds.size() != mol2.customBonds.size()
+                    || mol.customAngles.size() != mol2.customAngles.size() || mol.customTorsions.size() != mol2.customTorsions.size())
                 identical = false;
             int atomOffset = mol2.atoms[0]-mol.atoms[0];
             float4* posq = gpu->psPosq4->_pSysData;
@@ -2652,6 +2736,45 @@ static void findMoleculeGroups(gpuContext gpu)
                         lj14Param[mol.lj14s[i]].x != lj14Param[mol2.lj14s[i]].x || lj14Param[mol.lj14s[i]].y != lj14Param[mol2.lj14s[i]].y ||
                         lj14Param[mol.lj14s[i]].z != lj14Param[mol2.lj14s[i]].z)
                     identical = false;
+            if (mol.customBonds.size() > 0) {
+                int4* customBondID = gpu->psCustomBondID->_pSysData;
+                float4* customBondParam = gpu->psCustomBondParams->_pSysData;
+                for (int i = 0; i < (int)mol.customBonds.size() && identical; i++)
+                    if (customBondID[mol.customBonds[i]].x != customBondID[mol2.customBonds[i]].x-atomOffset ||
+                            customBondID[mol.customBonds[i]].y != customBondID[mol2.customBonds[i]].y-atomOffset ||
+                            (customBondParam[mol.customBonds[i]].x != customBondParam[mol2.customBonds[i]].x && gpu->sim.customBondParameters > 0) ||
+                            (customBondParam[mol.customBonds[i]].y != customBondParam[mol2.customBonds[i]].y && gpu->sim.customBondParameters > 1) ||
+                            (customBondParam[mol.customBonds[i]].z != customBondParam[mol2.customBonds[i]].z && gpu->sim.customBondParameters > 2) ||
+                            (customBondParam[mol.customBonds[i]].w != customBondParam[mol2.customBonds[i]].w && gpu->sim.customBondParameters > 3))
+                        identical = false;
+            }
+            if (mol.customAngles.size() > 0) {
+                int4* customAngleID = gpu->psCustomAngleID1->_pSysData;
+                float4* customAngleParam = gpu->psCustomAngleParams->_pSysData;
+                for (int i = 0; i < (int)mol.customAngles.size() && identical; i++)
+                    if (customAngleID[mol.customAngles[i]].x != customAngleID[mol2.customAngles[i]].x-atomOffset ||
+                            customAngleID[mol.customAngles[i]].y != customAngleID[mol2.customAngles[i]].y-atomOffset ||
+                            customAngleID[mol.customAngles[i]].z != customAngleID[mol2.customAngles[i]].z-atomOffset ||
+                            (customAngleParam[mol.customAngles[i]].x != customAngleParam[mol2.customAngles[i]].x && gpu->sim.customAngleParameters > 0) ||
+                            (customAngleParam[mol.customAngles[i]].y != customAngleParam[mol2.customAngles[i]].y && gpu->sim.customAngleParameters > 1) ||
+                            (customAngleParam[mol.customAngles[i]].z != customAngleParam[mol2.customAngles[i]].z && gpu->sim.customAngleParameters > 2) ||
+                            (customAngleParam[mol.customAngles[i]].w != customAngleParam[mol2.customAngles[i]].w && gpu->sim.customAngleParameters > 3))
+                        identical = false;
+            }
+            if (mol.customTorsions.size() > 0) {
+                int4* customTorsionID = gpu->psCustomTorsionID1->_pSysData;
+                float4* customTorsionParam = gpu->psCustomTorsionParams->_pSysData;
+                for (int i = 0; i < (int)mol.customTorsions.size() && identical; i++)
+                    if (customTorsionID[mol.customTorsions[i]].x != customTorsionID[mol2.customTorsions[i]].x-atomOffset ||
+                            customTorsionID[mol.customTorsions[i]].y != customTorsionID[mol2.customTorsions[i]].y-atomOffset ||
+                            customTorsionID[mol.customTorsions[i]].z != customTorsionID[mol2.customTorsions[i]].z-atomOffset ||
+                            customTorsionID[mol.customTorsions[i]].w != customTorsionID[mol2.customTorsions[i]].w-atomOffset ||
+                            (customTorsionParam[mol.customTorsions[i]].x != customTorsionParam[mol2.customTorsions[i]].x && gpu->sim.customTorsionParameters > 0) ||
+                            (customTorsionParam[mol.customTorsions[i]].y != customTorsionParam[mol2.customTorsions[i]].y && gpu->sim.customTorsionParameters > 1) ||
+                            (customTorsionParam[mol.customTorsions[i]].z != customTorsionParam[mol2.customTorsions[i]].z && gpu->sim.customTorsionParameters > 2) ||
+                            (customTorsionParam[mol.customTorsions[i]].w != customTorsionParam[mol2.customTorsions[i]].w && gpu->sim.customTorsionParameters > 3))
+                        identical = false;
+            }
             if (identical)
             {
                 moleculeInstances[j].push_back(mol.atoms[0]);

platforms/cuda/src/kernels/gputypes.h

@@ -107,6 +107,9 @@ struct _gpuContext {
     CUDAStream<int4>* psCustomAngleID1;           // Atom indices for custom angles
     CUDAStream<int2>* psCustomAngleID2;           // Atom indices for custom angles
     CUDAStream<float4>* psCustomAngleParams;      // Parameters for custom angles
+    CUDAStream<int4>* psCustomTorsionID1;           // Atom indices for custom torsions
+    CUDAStream<int4>* psCustomTorsionID2;           // Atom indices for custom torsions
+    CUDAStream<float4>* psCustomTorsionParams;      // Parameters for custom torsions
     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
@@ -216,7 +219,11 @@ void gpuSetCustomBondParameters(gpuContext gpu, const std::vector<int>& bondAtom
             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,
+void gpuSetCustomAngleParameters(gpuContext gpu, const std::vector<int>& angleAtom1, const std::vector<int>& angleAtom2, const std::vector<int>& angleAtom3, const std::vector<std::vector<double> >& angleParams,
+            const std::string& energyExp, const std::vector<std::string>& paramNames, const std::vector<std::string>& globalParamNames);
+
+extern "C"
+void gpuSetCustomTorsionParameters(gpuContext gpu, const std::vector<int>& torsionAtom1, const std::vector<int>& torsionAtom2, const std::vector<int>& torsionAtom3, const std::vector<int>& torsionAtom4, const std::vector<std::vector<double> >& torsionParams,
             const std::string& energyExp, const std::vector<std::string>& paramNames, const std::vector<std::string>& globalParamNames);
 
 extern "C"

platforms/cuda/src/kernels/kCalculateCustomTorsionForces.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 SetCalculateCustomTorsionForcesSim(gpuContext gpu)
+{
+    cudaError_t status;
+    status = cudaMemcpyToSymbol(cSim, &gpu->sim, sizeof(cudaGmxSimulation));
+    RTERROR(status, "cudaMemcpyToSymbol: SetSim copy to cSim failed");
+}
+
+void GetCalculateCustomTorsionForcesSim(gpuContext gpu)
+{
+    cudaError_t status;
+    status = cudaMemcpyFromSymbol(&gpu->sim, cSim, sizeof(cudaGmxSimulation));
+    RTERROR(status, "cudaMemcpyFromSymbol: SetSim copy from cSim failed");
+}
+
+void SetCustomTorsionForceExpression(const Expression<256>& expression)
+{
+    cudaError_t status;
+    status = cudaMemcpyToSymbol(forceExp, &expression, sizeof(forceExp));
+    RTERROR(status, "SetCustomTorsionForceExpression: cudaMemcpyToSymbol failed");
+}
+
+void SetCustomTorsionEnergyExpression(const Expression<256>& expression)
+{
+    cudaError_t status;
+    status = cudaMemcpyToSymbol(energyExp, &expression, sizeof(energyExp));
+    RTERROR(status, "SetCustomTorsionEnergyExpression: cudaMemcpyToSymbol failed");
+}
+
+void SetCustomTorsionGlobalParams(const vector<float>& paramValues)
+{
+    cudaError_t status;
+    status = cudaMemcpyToSymbol(globalParams, &paramValues[0], paramValues.size()*sizeof(float));
+    RTERROR(status, "SetCustomTorsionGlobalParams: cudaMemcpyToSymbol failed");
+}
+
+#define LOCAL_HACK_PI 3.1415926535897932384626433832795
+
+#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)
+
+#define GETNORMEDDOTPRODUCT(v1, v2, dp) \
+{ \
+    dp          = DOT3(v1, v2); \
+    float norm1 = DOT3(v1, v1); \
+    float norm2 = DOT3(v2, v2); \
+    dp /= sqrt(norm1 * norm2); \
+    dp = min(dp, 1.0f); \
+    dp = max(dp, -1.0f); \
+}
+
+#define GETANGLEBETWEENTWOVECTORS(v1, v2, angle) \
+{ \
+    float dp; \
+    GETNORMEDDOTPRODUCT(v1, v2, dp); \
+    if (dp > 0.99f || dp < -0.99f) { \
+        float3 cross = CROSS_PRODUCT(v1, v2); \
+        float scale = DOT3(v1, v1)*DOT3(v2, v2); \
+        angle = asin(sqrt(DOT3(cross, cross)/scale)); \
+        if (dp < 0.0f) \
+            angle = LOCAL_HACK_PI-angle; \
+    } \
+    else { \
+        angle = acos(dp); \
+    } \
+}
+
+#define GETDIHEDRALANGLEBETWEENTHREEVECTORS(vector1, vector2, vector3, signVector, cp0, cp1, angle) \
+{ \
+    cp0 = CROSS_PRODUCT(vector1, vector2); \
+    cp1 = CROSS_PRODUCT(vector2, vector3); \
+    GETANGLEBETWEENTWOVECTORS(cp0, cp1, angle); \
+    float dp = DOT3(signVector, cp1); \
+    angle = (dp >= 0) ? angle : -angle; \
+}
+
+__global__ void kCalculateCustomTorsionForces_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.customTorsions)
+    {
+        int4 atom       = cSim.pCustomTorsionID1[pos];
+        int4 atom2      = cSim.pCustomTorsionID2[pos];
+        float4 params   = cSim.pCustomTorsionParams[pos];
+        float4 a1       = cSim.pPosq[atom.x];
+        float4 a2       = cSim.pPosq[atom.y];
+        float4 a3       = cSim.pPosq[atom.z];
+        float4 a4       = cSim.pPosq[atom.w];
+        float3 v0 = make_float3(a1.x-a2.x, a1.y-a2.y, a1.z-a2.z);
+        float3 v1 = make_float3(a3.x-a2.x, a3.y-a2.y, a3.z-a2.z);
+        float3 v2 = make_float3(a3.x-a4.x, a3.y-a4.y, a3.z-a4.z);
+        float3 cp0, cp1;
+        float dihedralAngle;
+        GETDIHEDRALANGLEBETWEENTHREEVECTORS(v0, v1, v2, v0, cp0, cp1, dihedralAngle);
+        VARIABLE(0) = dihedralAngle;
+        VARIABLE(1) = params.x;
+        VARIABLE(2) = params.y;
+        VARIABLE(3) = params.z;
+        VARIABLE(4) = params.w;
+        float dEdAngle = kEvaluateExpression_kernel(&forceExp, stack, variables);
+        totalEnergy += kEvaluateExpression_kernel(&energyExp, stack, variables);
+        float normBC = sqrt(DOT3(v1, v1));
+        float dp = 1.0f / DOT3(v1, v1);
+        float4 ff = make_float4((-dEdAngle*normBC)/DOT3(cp0, cp0), DOT3(v0, v1)*dp, DOT3(v2, v1)*dp, (dEdAngle*normBC)/DOT3(cp1, cp1));
+        float3 internalF0 = make_float3(ff.x*cp0.x, ff.x*cp0.y, ff.x*cp0.z);
+        float3 internalF3 = make_float3(ff.w*cp1.x, ff.w*cp1.y, ff.w*cp1.z);
+        float3 s = make_float3(ff.y*internalF0.x - ff.z*internalF3.x,
+                               ff.y*internalF0.y - ff.z*internalF3.y,
+                               ff.y*internalF0.z - ff.z*internalF3.z);
+        unsigned int offsetA  = atom.x + atom2.x * cSim.stride;
+        unsigned int offsetB  = atom.y + atom2.y * cSim.stride;
+        unsigned int offsetC  = atom.z + atom2.z * cSim.stride;
+        unsigned int offsetD  = atom.w + atom2.w * cSim.stride;
+        float4 forceA         = cSim.pForce4[offsetA];
+        float4 forceB         = cSim.pForce4[offsetB];
+        float4 forceC         = cSim.pForce4[offsetC];
+        float4 forceD         = cSim.pForce4[offsetD];
+        forceA.x             += internalF0.x;
+        forceA.y             += internalF0.y;
+        forceA.z             += internalF0.z;
+        forceB.x             += -internalF0.x + s.x;
+        forceB.y             += -internalF0.y + s.y;
+        forceB.z             += -internalF0.z + s.z;
+        forceC.x             += -internalF3.x - s.x;
+        forceC.y             += -internalF3.y - s.y;
+        forceC.z             += -internalF3.z - s.z;
+        forceD.x             += internalF3.x;
+        forceD.y             += internalF3.y;
+        forceD.z             += internalF3.z;
+        cSim.pForce4[offsetA] = forceA;
+        cSim.pForce4[offsetB] = forceB;
+        cSim.pForce4[offsetC] = forceC;
+        cSim.pForce4[offsetD] = forceD;
+        pos                  += blockDim.x * gridDim.x;
+    }
+    cSim.pEnergy[blockIdx.x * blockDim.x + threadIdx.x] += totalEnergy;
+}
+
+void kCalculateCustomTorsionForces(gpuContext gpu)
+{
+//    printf("kCalculateCustomTorsionForces\n");
+    kCalculateCustomTorsionForces_kernel<<<gpu->sim.blocks, gpu->sim.localForces_threads_per_block,
+            (gpu->sim.customExpressionStackSize+9)*sizeof(float)*gpu->sim.localForces_threads_per_block>>>();
+    LAUNCHERROR("kCalculateCustomTorsionForces");
+}

platforms/cuda/tests/TestCudaCustomTorsionForce.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 CustomTorsionForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "CudaPlatform.h"
+#include "openmm/CustomTorsionForce.h"
+#include "openmm/PeriodicTorsionForce.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 testTorsions() {
+    CudaPlatform platform;
+
+    // Create a system using a CustomTorsionForce.
+
+    System customSystem;
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    CustomTorsionForce* custom = new CustomTorsionForce("k*(1+cos(n*theta-theta0))");
+    custom->addPerTorsionParameter("theta0");
+    custom->addPerTorsionParameter("n");
+    custom->addGlobalParameter("k", 0.5);
+    vector<double> parameters(2);
+    parameters[0] = 1.5;
+    parameters[1] = 1;
+    custom->addTorsion(0, 1, 2, 3, parameters);
+    parameters[0] = 2.0;
+    parameters[1] = 2;
+    custom->addTorsion(1, 2, 3, 4, parameters);
+    customSystem.addForce(custom);
+
+    // Create an identical system using a PeriodicTorsionForce.
+
+    System harmonicSystem;
+    harmonicSystem.addParticle(1.0);
+    harmonicSystem.addParticle(1.0);
+    harmonicSystem.addParticle(1.0);
+    harmonicSystem.addParticle(1.0);
+    harmonicSystem.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    PeriodicTorsionForce* periodic = new PeriodicTorsionForce();
+    periodic->addTorsion(0, 1, 2, 3, 1, 1.5, 0.5);
+    periodic->addTorsion(1, 2, 3, 4, 2, 2.0, 0.5);
+    harmonicSystem.addForce(periodic);
+
+    // Set the atoms in various positions, and verify that both systems give identical forces and energy.
+
+    init_gen_rand(0);
+    vector<Vec3> positions(5);
+    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();
+        }
+        for (int i = 0; i < customSystem.getNumParticles(); i++)
+            ASSERT_EQUAL_VEC(forces2[i], forces1[i], TOL);
+        ASSERT_EQUAL_TOL(energy2, energy1, TOL);
+    }
+}
+
+int main() {
+    try {
+        testTorsions();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+