Continuing to implement new CUDA platform: GBSAOBCForce

platforms/cuda2/src/CudaContext.cpp

@@ -179,6 +179,7 @@ CudaContext::CudaContext(const System& system, int deviceIndex, bool useBlocking
     compilationDefines["RECIP"] = useDoublePrecision ? "1.0/" : "1.0f/";
     compilationDefines["EXP"] = useDoublePrecision ? "exp" : "expf";
     compilationDefines["LOG"] = useDoublePrecision ? "log" : "logf";
+    compilationDefines["POW"] = useDoublePrecision ? "pow" : "powf";
     compilationDefines["COS"] = useDoublePrecision ? "cos" : "cosf";
     compilationDefines["SIN"] = useDoublePrecision ? "sin" : "sinf";
     compilationDefines["TAN"] = useDoublePrecision ? "tan" : "tanf";

platforms/cuda2/src/CudaKernelFactory.cpp

@@ -96,8 +96,8 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
         return new CudaCalcNonbondedForceKernel(name, platform, cu, context.getSystem());
     if (name == CalcCustomNonbondedForceKernel::Name())
         return new CudaCalcCustomNonbondedForceKernel(name, platform, cu, context.getSystem());
-//    if (name == CalcGBSAOBCForceKernel::Name())
-//        return new CudaCalcGBSAOBCForceKernel(name, platform, cu);
+    if (name == CalcGBSAOBCForceKernel::Name())
+        return new CudaCalcGBSAOBCForceKernel(name, platform, cu);
 //    if (name == CalcCustomGBForceKernel::Name())
 //        return new CudaCalcCustomGBForceKernel(name, platform, cu, context.getSystem());
     if (name == CalcCustomExternalForceKernel::Name())

platforms/cuda2/src/CudaKernels.h

@@ -666,57 +666,55 @@ private:
     System& system;
 };
 
-///**
-// * This kernel is invoked by GBSAOBCForce to calculate the forces acting on the system.
-// */
-//class CudaCalcGBSAOBCForceKernel : public CalcGBSAOBCForceKernel {
-//public:
-//    CudaCalcGBSAOBCForceKernel(std::string name, const Platform& platform, CudaContext& cu) : CalcGBSAOBCForceKernel(name, platform), cu(cu),
-//            hasCreatedKernels(false), params(NULL), bornSum(NULL), longBornSum(NULL), bornRadii(NULL), bornForce(NULL),
-//            longBornForce(NULL), obcChain(NULL) {
-//    }
-//    ~CudaCalcGBSAOBCForceKernel();
-//    /**
-//     * Initialize the kernel.
-//     *
-//     * @param system     the System this kernel will be applied to
-//     * @param force      the GBSAOBCForce this kernel will be used for
-//     */
-//    void initialize(const System& system, const GBSAOBCForce& force);
-//    /**
-//     * Execute the kernel to calculate the forces and/or energy.
-//     *
-//     * @param context        the context in which to execute this kernel
-//     * @param includeForces  true if forces should be calculated
-//     * @param includeEnergy  true if the energy should be calculated
-//     * @return the potential energy due to the force
-//     */
-//    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
-//    /**
-//     * Copy changed parameters over to a context.
-//     *
-//     * @param context    the context to copy parameters to
-//     * @param force      the GBSAOBCForce to copy the parameters from
-//     */
-//    void copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force);
-//private:
-//    double prefactor;
-//    bool hasCreatedKernels;
-//    int maxTiles;
-//    CudaContext& cu;
-//    CudaArray<mm_float2>* params;
-//    CudaArray<cl_float>* bornSum;
-//    CudaArray<cl_long>* longBornSum;
-//    CudaArray<cl_float>* bornRadii;
-//    CudaArray<cl_float>* bornForce;
-//    CudaArray<cl_long>* longBornForce;
-//    CudaArray<cl_float>* obcChain;
-//    CUfunction computeBornSumKernel;
-//    CUfunction reduceBornSumKernel;
-//    CUfunction force1Kernel;
-//    CUfunction reduceBornForceKernel;
-//};
-//
+/**
+ * This kernel is invoked by GBSAOBCForce to calculate the forces acting on the system.
+ */
+class CudaCalcGBSAOBCForceKernel : public CalcGBSAOBCForceKernel {
+public:
+    CudaCalcGBSAOBCForceKernel(std::string name, const Platform& platform, CudaContext& cu) : CalcGBSAOBCForceKernel(name, platform), cu(cu),
+            hasCreatedKernels(false), params(NULL), bornSum(NULL), bornRadii(NULL), bornForce(NULL), obcChain(NULL) {
+    }
+    ~CudaCalcGBSAOBCForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the GBSAOBCForce this kernel will be used for
+     */
+    void initialize(const System& system, const GBSAOBCForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the GBSAOBCForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force);
+private:
+    double prefactor;
+    bool hasCreatedKernels;
+    int maxTiles;
+    CudaContext& cu;
+    CudaArray* params;
+    CudaArray* bornSum;
+    CudaArray* bornRadii;
+    CudaArray* bornForce;
+    CudaArray* obcChain;
+    CUfunction computeBornSumKernel;
+    CUfunction reduceBornSumKernel;
+    CUfunction force1Kernel;
+    CUfunction reduceBornForceKernel;
+    std::vector<void*> computeSumArgs, force1Args;
+};
+
 ///**
 // * This kernel is invoked by CustomGBForce to calculate the forces acting on the system.
 // */

platforms/cuda2/src/kernels/gbsaObc2.cu

+{
+    real invRSquaredOver4 = 0.25f*invR*invR;
+    real rScaledRadiusJ = r+obcParams2.y;
+    real rScaledRadiusI = r+obcParams1.y;
+    real l_ijJ = RECIP(max(obcParams1.x, fabs(r-obcParams2.y)));
+    real l_ijI = RECIP(max(obcParams2.x, fabs(r-obcParams1.y)));
+    real u_ijJ = RECIP(rScaledRadiusJ);
+    real u_ijI = RECIP(rScaledRadiusI);
+    real l_ij2J = l_ijJ*l_ijJ;
+    real l_ij2I = l_ijI*l_ijI;
+    real u_ij2J = u_ijJ*u_ijJ;
+    real u_ij2I = u_ijI*u_ijI;
+    real t1J = LOG(u_ijJ*RECIP(l_ijJ));
+    real t1I = LOG(u_ijI*RECIP(l_ijI));
+    real t2J = (l_ij2J-u_ij2J);
+    real t2I = (l_ij2I-u_ij2I);
+    real term1 = (0.5f*(0.25f+obcParams2.y*obcParams2.y*invRSquaredOver4)*t2J + t1J*invRSquaredOver4)*invR;
+    real term2 = (0.5f*(0.25f+obcParams1.y*obcParams1.y*invRSquaredOver4)*t2I + t1I*invRSquaredOver4)*invR;
+    real tempdEdR = (obcParams1.x < rScaledRadiusJ ? bornForce1*term1/0xFFFFFFFF : 0);
+    tempdEdR += (obcParams2.x < rScaledRadiusI ? bornForce2*term2/0xFFFFFFFF : 0);
+#ifdef USE_CUTOFF
+    unsigned int includeInteraction = (atom1 < NUM_ATOMS && atom2 < NUM_ATOMS && atom1 != atom2 && r2 < CUTOFF_SQUARED);
+#else
+    unsigned int includeInteraction = (atom1 < NUM_ATOMS && atom2 < NUM_ATOMS && atom1 != atom2);
+#endif
+    dEdR += (includeInteraction ? tempdEdR : 0);
+}

platforms/cuda2/src/kernels/gbsaObcReductions.cu

+#define DIELECTRIC_OFFSET 0.009f
+#define PROBE_RADIUS 0.14f
+#define SURFACE_AREA_FACTOR -170.351730667551f //-6.0f*3.14159265358979323846f*0.0216f*1000.0f*0.4184f;
+
+/**
+ * Reduce the Born sums to compute the Born radii.
+ */
+
+extern "C" __global__ void reduceBornSum(float alpha, float beta, float gamma, const long long* __restrict__ bornSum,
+            const float2* __restrict__ params, real* __restrict__ bornRadii, real* __restrict__ obcChain) {
+    for (unsigned int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_ATOMS; index += blockDim.x*gridDim.x) {
+        // Get summed Born data
+
+        real sum = RECIP(0xFFFFFFFF)*bornSum[index];
+
+        // Now calculate Born radius and OBC term.
+
+        float offsetRadius = params[index].x;
+        sum *= 0.5f*offsetRadius;
+        real sum2 = sum*sum;
+        real sum3 = sum*sum2;
+        real tanhSum = tanh(alpha*sum - beta*sum2 + gamma*sum3);
+        real nonOffsetRadius = offsetRadius + DIELECTRIC_OFFSET;
+        real radius = RECIP(RECIP(offsetRadius) - tanhSum/nonOffsetRadius);
+        real chain = offsetRadius*(alpha - 2.0f*beta*sum + 3.0f*gamma*sum2);
+        chain = (1-tanhSum*tanhSum)*chain / nonOffsetRadius;
+        bornRadii[index] = radius;
+        obcChain[index] = chain;
+    }
+}
+
+/**
+ * Reduce the Born force.
+ */
+
+extern "C" __global__ void reduceBornForce(long long* __restrict__ bornForce, real* __restrict__ energyBuffer,
+        const float2* __restrict__ params, const real* __restrict__ bornRadii, const real* __restrict__ obcChain) {
+    real energy = 0;
+    for (unsigned int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_ATOMS; index += blockDim.x*gridDim.x) {
+        // Get summed Born force
+
+        real force = RECIP(0xFFFFFFFF)*bornForce[index];
+
+        // Now calculate the actual force
+
+        float offsetRadius = params[index].x;
+        real bornRadius = bornRadii[index];
+        real r = offsetRadius+DIELECTRIC_OFFSET+PROBE_RADIUS;
+        real ratio6 = POW((offsetRadius+DIELECTRIC_OFFSET)/bornRadius, 6);
+        real saTerm = SURFACE_AREA_FACTOR*r*r*ratio6;
+        force += saTerm/bornRadius;
+        energy += saTerm;
+        force *= bornRadius*bornRadius*obcChain[index];
+        bornForce[index] = (long long) (force*0xFFFFFFFF);
+    }
+    energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] += energy/-6;
+}
\ No newline at end of file

platforms/cuda2/tests/TestCudaGBSAOBCForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the CUDA implementation of GBSAOBCForce.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "CudaPlatform.h"
+#include "ReferencePlatform.h"
+#include "openmm/GBSAOBCForce.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include "sfmt/SFMT.h"
+#include "openmm/NonbondedForce.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testSingleParticle() {
+    CudaPlatform platform;
+    System system;
+    system.addParticle(2.0);
+    LangevinIntegrator integrator(0, 0.1, 0.01);
+    GBSAOBCForce* gbsa = new GBSAOBCForce();
+    NonbondedForce* nonbonded = new NonbondedForce();
+    gbsa->addParticle(0.5, 0.15, 1);
+    nonbonded->addParticle(0.5, 1, 0);
+    system.addForce(gbsa);
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(1);
+    positions[0] = Vec3(0, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Energy);
+    double bornRadius = 0.15-0.009; // dielectric offset
+    double eps0 = EPSILON0;
+    double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/gbsa->getSoluteDielectric()-1.0/gbsa->getSolventDielectric())/bornRadius;
+    double extendedRadius = bornRadius+0.14; // probe radius
+    double nonpolarEnergy = CAL2JOULE*PI_M*0.0216*(10*extendedRadius)*(10*extendedRadius)*std::pow(0.15/bornRadius, 6.0); // Where did this formula come from?  Just copied it from CpuImplicitSolvent.cpp
+    ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
+    
+    // Change the parameters and see if it is still correct.
+    
+    gbsa->setParticleParameters(0, 0.4, 0.25, 1);
+    gbsa->updateParametersInContext(context);
+    state = context.getState(State::Energy);
+    bornRadius = 0.25-0.009; // dielectric offset
+    bornEnergy = (-0.4*0.4/(8*PI_M*eps0))*(1.0/gbsa->getSoluteDielectric()-1.0/gbsa->getSolventDielectric())/bornRadius;
+    extendedRadius = bornRadius+0.14;
+    nonpolarEnergy = CAL2JOULE*PI_M*0.0216*(10*extendedRadius)*(10*extendedRadius)*std::pow(0.25/bornRadius, 6.0);
+    ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
+}
+
+void testCutoffAndPeriodic() {
+    CudaPlatform cl;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    LangevinIntegrator integrator(0, 0.1, 0.01);
+    GBSAOBCForce* gbsa = new GBSAOBCForce();
+    NonbondedForce* nonbonded = new NonbondedForce();
+    gbsa->addParticle(-1, 0.15, 1);
+    nonbonded->addParticle(-1, 1, 0);
+    gbsa->addParticle(1, 0.15, 1);
+    nonbonded->addParticle(1, 1, 0);
+    const double cutoffDistance = 3.0;
+    const double boxSize = 10.0;
+    nonbonded->setCutoffDistance(cutoffDistance);
+    gbsa->setCutoffDistance(cutoffDistance);
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.addForce(gbsa);
+    system.addForce(nonbonded);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(2, 0, 0);
+
+    // Calculate the forces for both cutoff and periodic with two different atom positions.
+
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
+    Context context(system, integrator, cl);
+    context.setPositions(positions);
+    State state1 = context.getState(State::Forces);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
+    context.reinitialize();
+    context.setPositions(positions);
+    State state2 = context.getState(State::Forces);
+    positions[1][0]+= boxSize;
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
+    context.reinitialize();
+    context.setPositions(positions);
+    State state3 = context.getState(State::Forces);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
+    context.reinitialize();
+    context.setPositions(positions);
+    State state4 = context.getState(State::Forces);
+
+    // All forces should be identical, exception state3 which should be zero.
+
+    ASSERT_EQUAL_VEC(state1.getForces()[0], state2.getForces()[0], 0.01);
+    ASSERT_EQUAL_VEC(state1.getForces()[1], state2.getForces()[1], 0.01);
+    ASSERT_EQUAL_VEC(state1.getForces()[0], state4.getForces()[0], 0.01);
+    ASSERT_EQUAL_VEC(state1.getForces()[1], state4.getForces()[1], 0.01);
+    ASSERT_EQUAL_VEC(state3.getForces()[0], Vec3(0, 0, 0), 0.01);
+    ASSERT_EQUAL_VEC(state3.getForces()[1], Vec3(0, 0, 0), 0.01);
+}
+
+void testForce(int numParticles, NonbondedForce::NonbondedMethod method, GBSAOBCForce::NonbondedMethod method2) {
+    CudaPlatform cl;
+    ReferencePlatform reference;
+    System system;
+    GBSAOBCForce* gbsa = new GBSAOBCForce();
+    NonbondedForce* nonbonded = new NonbondedForce();
+    for (int i = 0; i < numParticles; ++i) {
+        system.addParticle(1.0);
+        double charge = i%2 == 0 ? -1 : 1;
+        gbsa->addParticle(charge, 0.15, 1);
+        nonbonded->addParticle(charge, 1, 0);
+    }
+    nonbonded->setNonbondedMethod(method);
+    gbsa->setNonbondedMethod(method2);
+    nonbonded->setCutoffDistance(3.0);
+    gbsa->setCutoffDistance(3.0);
+    int grid = (int) floor(0.5+pow(numParticles, 1.0/3.0));
+    if (method == NonbondedForce::CutoffPeriodic) {
+        double boxSize = (grid+1)*1.1;
+        system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    }
+    system.addForce(gbsa);
+    system.addForce(nonbonded);
+    LangevinIntegrator integrator1(0, 0.1, 0.01);
+    LangevinIntegrator integrator2(0, 0.1, 0.01);
+    Context context(system, integrator1, cl);
+    Context refContext(system, integrator2, reference);
+
+    // Set random (but uniformly distributed) positions for all the particles.
+
+    vector<Vec3> positions(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+
+    for (int i = 0; i < grid; i++)
+        for (int j = 0; j < grid; j++)
+            for (int k = 0; k < grid; k++)
+                positions[i*grid*grid+j*grid+k] = Vec3(i*1.1, j*1.1, k*1.1);
+    for (int i = 0; i < numParticles; ++i)
+        positions[i] = positions[i] + Vec3(0.5*genrand_real2(sfmt), 0.5*genrand_real2(sfmt), 0.5*genrand_real2(sfmt));
+    context.setPositions(positions);
+    refContext.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    State refState = refContext.getState(State::Forces | State::Energy);
+
+    // Make sure the CUDA and Reference platforms agree.
+
+    double norm = 0.0;
+    double diff = 0.0;
+    for (int i = 0; i < numParticles; ++i) {
+        Vec3 f = state.getForces()[i];
+        norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
+        Vec3 delta = f-refState.getForces()[i];
+        diff += delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2];
+    }
+    norm = std::sqrt(norm);
+    diff = std::sqrt(diff);
+    ASSERT_EQUAL_TOL(0.0, diff, 0.001*norm);
+    ASSERT_EQUAL_TOL(state.getPotentialEnergy(), refState.getPotentialEnergy(), 1e-3);
+
+    // Take a small step in the direction of the energy gradient.  (This doesn't work with cutoffs, since the energy
+    // changes discontinuously at the cutoff distance.)
+
+    if (method == NonbondedForce::NoCutoff)
+    {
+        const double delta = 1e-2;
+        double step = delta/norm;
+        for (int i = 0; i < numParticles; ++i) {
+            Vec3 p = positions[i];
+            Vec3 f = state.getForces()[i];
+            positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
+        }
+        context.setPositions(positions);
+
+        // See whether the potential energy changed by the expected amount.
+
+        State state2 = context.getState(State::Energy);
+        ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state.getPotentialEnergy())/delta, 1e-3*abs(state.getPotentialEnergy()));
+    }
+}
+
+int main() {
+    try {
+        testSingleParticle();
+        testCutoffAndPeriodic();
+        for (int i = 5; i < 11; i++) {
+            testForce(i*i*i, NonbondedForce::NoCutoff, GBSAOBCForce::NoCutoff);
+            testForce(i*i*i, NonbondedForce::CutoffNonPeriodic, GBSAOBCForce::CutoffNonPeriodic);
+            testForce(i*i*i, NonbondedForce::CutoffPeriodic, GBSAOBCForce::CutoffPeriodic);
+        }
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+