Began CUDA version of CustomCentroidBondForce

platforms/cuda/include/CudaKernels.h

@@ -922,6 +922,58 @@ private:
     CUfunction donorKernel, acceptorKernel;
 };
 
+/**
+ * This kernel is invoked by CustomCentroidBondForce to calculate the forces acting on the system.
+ */
+class CudaCalcCustomCentroidBondForceKernel : public CalcCustomCentroidBondForceKernel {
+public:
+    CudaCalcCustomCentroidBondForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system) : CalcCustomCentroidBondForceKernel(name, platform),
+            cu(cu), params(NULL), globals(NULL), groupParticles(NULL), groupWeights(NULL), groupOffsets(NULL), groupForces(NULL), bondGroups(NULL), centerPositions(NULL), system(system) {
+    }
+    ~CudaCalcCustomCentroidBondForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomCentroidBondForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomCentroidBondForce& 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 CustomCentroidBondForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomCentroidBondForce& force);
+
+private:
+    int numGroups, numBonds;
+    CudaContext& cu;
+    CudaParameterSet* params;
+    CudaArray* globals;
+    CudaArray* groupParticles;
+    CudaArray* groupWeights;
+    CudaArray* groupOffsets;
+    CudaArray* groupForces;
+    CudaArray* bondGroups;
+    CudaArray* centerPositions;
+    std::vector<std::string> globalParamNames;
+    std::vector<float> globalParamValues;
+    std::vector<CudaArray*> tabulatedFunctions;
+    std::vector<void*> groupForcesArgs;
+    CUfunction computeCentersKernel, groupForcesKernel, applyForcesKernel;
+    const System& system;
+};
+
 /**
  * This kernel is invoked by CustomCompoundBondForce to calculate the forces acting on the system.
  */

platforms/cuda/src/CudaKernelFactory.cpp

@@ -104,6 +104,8 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
         return new CudaCalcCustomExternalForceKernel(name, platform, cu, context.getSystem());
     if (name == CalcCustomHbondForceKernel::Name())
         return new CudaCalcCustomHbondForceKernel(name, platform, cu, context.getSystem());
+    if (name == CalcCustomCentroidBondForceKernel::Name())
+        return new CudaCalcCustomCentroidBondForceKernel(name, platform, cu, context.getSystem());
     if (name == CalcCustomCompoundBondForceKernel::Name())
         return new CudaCalcCustomCompoundBondForceKernel(name, platform, cu, context.getSystem());
     if (name == CalcCustomManyParticleForceKernel::Name())

platforms/cuda/src/CudaPlatform.cpp

@@ -80,6 +80,7 @@ CudaPlatform::CudaPlatform() {
     registerKernelFactory(CalcCustomGBForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomExternalForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomHbondForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomCentroidBondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCompoundBondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);

platforms/cuda/src/kernels/customCentroidBond.cu

+/**
+ * Compute the center of each group.
+ */
+extern "C" __global__ void computeGroupCenters(const real4* __restrict__ posq, const int* __restrict__ groupParticles,
+        const real* __restrict__ groupWeights, const int* __restrict__ groupOffsets, real4* __restrict__ centerPositions) {
+    __shared__ volatile real3 temp[64];
+    for (int group = blockIdx.x; group < NUM_GROUPS; group += gridDim.x) {
+        // The threads in this block work together to compute the center one group.
+        
+        int firstIndex = groupOffsets[group];
+        int lastIndex = groupOffsets[group+1];
+        real3 center = make_real3(0, 0, 0);
+        for (int index = threadIdx.x; index < lastIndex-firstIndex; index += blockDim.x) {
+            int atom = groupParticles[firstIndex+index];
+            real weight = groupWeights[firstIndex+index];
+            real4 pos = posq[atom];
+            center.x += weight*pos.x;
+            center.y += weight*pos.y;
+            center.z += weight*pos.z;
+        }
+        
+        // Sum the values.
+        
+        int thread = threadIdx.x;
+        temp[thread].x = center.x;
+        temp[thread].y = center.y;
+        temp[thread].z = center.z;
+        __syncthreads();
+        if (thread < 32) {
+            temp[thread].x += temp[thread+32].x;
+            temp[thread].y += temp[thread+32].y;
+            temp[thread].z += temp[thread+32].z;
+            if (thread < 16) {
+                temp[thread].x += temp[thread+16].x;
+                temp[thread].y += temp[thread+16].y;
+                temp[thread].z += temp[thread+16].z;
+            }
+            if (thread < 8) {
+                temp[thread].x += temp[thread+8].x;
+                temp[thread].y += temp[thread+8].y;
+                temp[thread].z += temp[thread+8].z;
+            }
+            if (thread < 4) {
+                temp[thread].x += temp[thread+4].x;
+                temp[thread].y += temp[thread+4].y;
+                temp[thread].z += temp[thread+4].z;
+            }
+            if (thread < 2) {
+                temp[thread].x += temp[thread+2].x;
+                temp[thread].y += temp[thread+2].y;
+                temp[thread].z += temp[thread+2].z;
+            }
+        }
+        if (thread == 0)
+            centerPositions[group] = make_real4(temp[0].x+temp[1].x, temp[0].y+temp[1].y, temp[0].z+temp[1].z, 0);
+    }
+}
+
+/**
+ * Convert a real4 to a real3 by removing its last element.
+ */
+inline __device__ real3 trim(real4 v) {
+    return make_real3(v.x, v.y, v.z);
+}
+
+/**
+ * Compute the difference between two vectors, setting the fourth component to the squared magnitude.
+ */
+inline __device__ real4 delta(real4 vec1, real4 vec2) {
+    real4 result = make_real4(vec1.x-vec2.x, vec1.y-vec2.y, vec1.z-vec2.z, 0);
+    result.w = result.x*result.x + result.y*result.y + result.z*result.z;
+    return result;
+}
+
+/**
+ * Compute the angle between two vectors.  The w component of each vector should contain the squared magnitude.
+ */
+__device__ real computeAngle(real4 vec1, real4 vec2) {
+    real dotProduct = vec1.x*vec2.x + vec1.y*vec2.y + vec1.z*vec2.z;
+    real cosine = dotProduct*RSQRT(vec1.w*vec2.w);
+    real angle;
+    if (cosine > 0.99f || cosine < -0.99f) {
+        // We're close to the singularity in acos(), so take the cross product and use asin() instead.
+
+        real3 crossProduct = cross(vec1, vec2);
+        real scale = vec1.w*vec2.w;
+        angle = ASIN(SQRT(dot(crossProduct, crossProduct)/scale));
+        if (cosine < 0.0f)
+            angle = M_PI-angle;
+    }
+    else
+       angle = ACOS(cosine);
+    return angle;
+}
+
+/**
+ * Compute the cross product of two vectors, setting the fourth component to the squared magnitude.
+ */
+inline __device__ real4 computeCross(real4 vec1, real4 vec2) {
+    real3 cp = cross(vec1, vec2);
+    return make_real4(cp.x, cp.y, cp.z, cp.x*cp.x+cp.y*cp.y+cp.z*cp.z);
+}
+
+/**
+ * Compute the forces on groups based on the bonds.
+ */
+extern "C" __global__ void computeGroupForces(unsigned long long* __restrict__ groupForce, real* __restrict__ energyBuffer, const real4* __restrict__ centerPositions,
+        const int* __restrict__ bondGroups
+        EXTRA_ARGS) {
+    extern __shared__ real4 posBuffer[];
+    real energy = 0;
+    for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_BONDS; index += blockDim.x*gridDim.x) {
+        COMPUTE_FORCE
+    }
+    energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] += energy;
+}
+
+/**
+ * Apply the forces from the group centers to the individual atoms.
+ */
+extern "C" __global__ void applyForcesToAtoms(const int* __restrict__ groupParticles, const real* __restrict__ groupWeights, const int* __restrict__ groupOffsets,
+        const long long* __restrict__ groupForce, unsigned long long* __restrict__ atomForce) {
+    for (int group = blockIdx.x; group < NUM_GROUPS; group += gridDim.x) {
+        long long fx = groupForce[group];
+        long long fy = groupForce[group+NUM_GROUPS];
+        long long fz = groupForce[group+NUM_GROUPS*2];
+        int firstIndex = groupOffsets[group];
+        int lastIndex = groupOffsets[group+1];
+        for (int index = threadIdx.x; index < lastIndex-firstIndex; index += blockDim.x) {
+            int atom = groupParticles[firstIndex+index];
+            real weight = groupWeights[firstIndex+index];
+            atomicAdd(&atomForce[atom], static_cast<unsigned long long>((long long) (fx*weight)));
+            atomicAdd(&atomForce[atom+PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (fy*weight)));
+            atomicAdd(&atomForce[atom+2*PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (fz*weight)));
+        }
+    }
+}

platforms/cuda/tests/TestCudaCustomCentroidBondForce.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) 2015 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 reference implementation of CustomCompoundBondForce.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "CudaPlatform.h"
+#include "openmm/CustomCentroidBondForce.h"
+#include "openmm/CustomCompoundBondForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+CudaPlatform platform;
+
+const double TOL = 1e-5;
+
+void testHarmonicBond() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(2.0);
+    system.addParticle(3.0);
+    system.addParticle(4.0);
+    system.addParticle(5.0);
+    CustomCentroidBondForce* force = new CustomCentroidBondForce(2, "k*distance(g1,g2)^2");
+    force->addPerBondParameter("k");
+    vector<int> particles1;
+    particles1.push_back(0);
+    particles1.push_back(1);
+    vector<int> particles2;
+    particles2.push_back(2);
+    particles2.push_back(3);
+    particles2.push_back(4);
+    force->addGroup(particles1);
+    force->addGroup(particles2);
+    vector<int> groups;
+    groups.push_back(0);
+    groups.push_back(1);
+    vector<double> parameters;
+    parameters.push_back(1.0);
+    force->addBond(groups, parameters);
+    system.addForce(force);
+    ASSERT(!system.usesPeriodicBoundaryConditions());
+
+    // The center of mass of group 0 is (1.5, 0, 0).
+
+    vector<Vec3> positions(5);
+    positions[0] = Vec3(2.5, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+
+    // The center of mass of group 1 is (-1, 0, 0).
+
+    positions[2] = Vec3(-6, 0, 0);
+    positions[3] = Vec3(-1, 0, 0);
+    positions[4] = Vec3(2, 0, 0);
+
+    // Check the forces and energy.
+
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(2.5*2.5, state.getPotentialEnergy(), TOL);
+    ASSERT_EQUAL_VEC(Vec3(-2*2.5*(1.0/3.0), 0, 0), state.getForces()[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-2*2.5*(2.0/3.0), 0, 0), state.getForces()[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(2*2.5*(3.0/12.0), 0, 0), state.getForces()[2], TOL);
+    ASSERT_EQUAL_VEC(Vec3(2*2.5*(4.0/12.0), 0, 0), state.getForces()[3], TOL);
+    ASSERT_EQUAL_VEC(Vec3(2*2.5*(5.0/12.0), 0, 0), state.getForces()[4], TOL);
+
+    // Update the per-bond parameter and see if the results change.
+
+    parameters[0] = 2.0;
+    force->setBondParameters(0, groups, parameters);
+    force->updateParametersInContext(context);
+    state = context.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(2*2.5*2.5, state.getPotentialEnergy(), TOL);
+    ASSERT_EQUAL_VEC(Vec3(-4*2.5*(1.0/3.0), 0, 0), state.getForces()[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-4*2.5*(2.0/3.0), 0, 0), state.getForces()[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(4*2.5*(3.0/12.0), 0, 0), state.getForces()[2], TOL);
+    ASSERT_EQUAL_VEC(Vec3(4*2.5*(4.0/12.0), 0, 0), state.getForces()[3], TOL);
+    ASSERT_EQUAL_VEC(Vec3(4*2.5*(5.0/12.0), 0, 0), state.getForces()[4], TOL);
+}
+
+void testComplexFunction() {
+    int numParticles = 4;
+    System system;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(2.0);
+
+    // When every group contains only one particle, a CustomCentroidBondForce is identical to a
+    // CustomCompoundBondForce.  Use that to test a complicated energy function with lots of terms.
+
+    CustomCompoundBondForce* compound = new CustomCompoundBondForce(4, "x1+y2+z4+distance(p1,p2)*angle(p3,p2,p4)+0.5*dihedral(p2,p1,p4,p3)");
+    CustomCentroidBondForce* centroid = new CustomCentroidBondForce(4, "x1+y2+z4+distance(g1,g2)*angle(g3,g2,g4)+0.5*dihedral(g2,g1,g4,g3)");
+
+    // Add a single bond to the CustomCompoundBondForce.
+
+    vector<int> particles(4);
+    vector<double> parameters;
+    particles[0] = 0;
+    particles[1] = 1;
+    particles[2] = 2;
+    particles[3] = 3;
+    compound->addBond(particles, parameters);
+
+    // Add an identical bond to the CustomCentroidBondForce.  As a stronger test, make sure that
+    // group number is different from particle number.
+
+    vector<int> groupMembers(1);
+    groupMembers[0] = 3;
+    centroid->addGroup(groupMembers);
+    groupMembers[0] = 0;
+    centroid->addGroup(groupMembers);
+    groupMembers[0] = 1;
+    centroid->addGroup(groupMembers);
+    groupMembers[0] = 2;
+    centroid->addGroup(groupMembers);
+    vector<int> groups(4);
+    groups[0] = 1;
+    groups[1] = 2;
+    groups[2] = 3;
+    groups[3] = 0;
+    centroid->addBond(groups, parameters);
+
+    // Add both forces as different force groups, and create a context.
+
+    centroid->setForceGroup(1);
+    system.addForce(compound);
+    system.addForce(centroid);
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+
+    // Evaluate the force and energy for various positions and see if they match.
+
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < 10; i++) {
+        for (int j = 0; j < numParticles; j++)
+            positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
+        context.setPositions(positions);
+        State state1 = context.getState(State::Forces | State::Energy, false, 1<<0);
+        State state2 = context.getState(State::Forces | State::Energy, false, 1<<1);
+        ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), TOL);
+        for (int i = 0; i < numParticles; i++)
+            ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], TOL);
+    }
+}
+
+void testCustomWeights() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(2.0);
+    system.addParticle(3.0);
+    system.addParticle(4.0);
+    CustomCentroidBondForce* force = new CustomCentroidBondForce(2, "distance(g1,g2)^2");
+    vector<int> particles(2);
+    vector<double> weights(2);
+    particles[0] = 0;
+    particles[1] = 1;
+    weights[0] = 0.5;
+    weights[1] = 1.5;
+    force->addGroup(particles, weights);
+    particles[0] = 2;
+    particles[1] = 3;
+    weights[0] = 2.0;
+    weights[1] = 1.0;
+    force->addGroup(particles, weights);
+    vector<int> groups;
+    groups.push_back(0);
+    groups.push_back(1);
+    vector<double> parameters;
+    force->addBond(groups, parameters);
+    system.addForce(force);
+
+    // The center of mass of group 0 is (0, 1, 0).
+
+    vector<Vec3> positions(4);
+    positions[0] = Vec3(0, 4, 0);
+    positions[1] = Vec3(0, 0, 0);
+
+    // The center of mass of group 1 is (0, 10, 0).
+
+    positions[2] = Vec3(0, 9, 0);
+    positions[3] = Vec3(0, 12, 0);
+
+    // Check the forces and energy.
+
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(9*9, state.getPotentialEnergy(), TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 2*9*(0.5/2.0), 0), state.getForces()[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 2*9*(1.5/2.0), 0), state.getForces()[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, -2*9*(2.0/3.0), 0), state.getForces()[2], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, -2*9*(1.0/3.0), 0), state.getForces()[3], TOL);
+}
+
+int main(int argc, char* argv[]) {
+    try {
+        if (argc > 1)
+            platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
+        testHarmonicBond();
+        testComplexFunction();
+        testCustomWeights();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}