Merge pull request #1784 from peastman/hbond

Fixed bug accessing parameters for CustomHbondForce

platforms/cuda/src/CudaKernels.cpp

@@ -4521,12 +4521,12 @@ void CudaCalcCustomHbondForceKernel::initialize(const System& system, const Cust
     for (int i = 0; i < (int) donorParams->getBuffers().size(); i++) {
         CudaNonbondedUtilities::ParameterInfo& buffer = donorParams->getBuffers()[i];
         extraArgs << ", const "+buffer.getType()+"* __restrict__ donor"+buffer.getName();
-        addDonorAndAcceptorCode(computeDonor, computeAcceptor, buffer.getType()+" donorParams"+cu.intToString(i+1)+" = donor"+buffer.getName()+"[index];\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, buffer.getType()+" donorParams"+cu.intToString(i+1)+" = donor"+buffer.getName()+"[donorIndex];\n");
     }
     for (int i = 0; i < (int) acceptorParams->getBuffers().size(); i++) {
         CudaNonbondedUtilities::ParameterInfo& buffer = acceptorParams->getBuffers()[i];
         extraArgs << ", const "+buffer.getType()+"* __restrict__ acceptor"+buffer.getName();
-        addDonorAndAcceptorCode(computeDonor, computeAcceptor, buffer.getType()+" acceptorParams"+cu.intToString(i+1)+" = acceptor"+buffer.getName()+"[index];\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, buffer.getType()+" acceptorParams"+cu.intToString(i+1)+" = acceptor"+buffer.getName()+"[acceptorIndex];\n");
     }
 
     // Now evaluate the expressions.

platforms/cuda/src/kernels/customHbondForce.cu

@@ -105,8 +105,8 @@ extern "C" __global__ void computeDonorForces(unsigned long long* __restrict__ f
             __syncthreads();
             if (donorIndex < NUM_DONORS) {
                 for (int index = 0; index < blockSize; index++) {
-#ifdef USE_EXCLUSIONS
                     int acceptorIndex = acceptorStart+index;
+#ifdef USE_EXCLUSIONS
                     if (acceptorIndex == exclusionIndices.x || acceptorIndex == exclusionIndices.y || acceptorIndex == exclusionIndices.z || acceptorIndex == exclusionIndices.w)
                         continue;
 #endif
@@ -193,8 +193,8 @@ extern "C" __global__ void computeAcceptorForces(unsigned long long* __restrict_
             __syncthreads();
             if (acceptorIndex < NUM_ACCEPTORS) {
                 for (int index = 0; index < blockSize; index++) {
-#ifdef USE_EXCLUSIONS
                     int donorIndex = donorStart+index;
+#ifdef USE_EXCLUSIONS
                     if (donorIndex == exclusionIndices.x || donorIndex == exclusionIndices.y || donorIndex == exclusionIndices.z || donorIndex == exclusionIndices.w)
                         continue;
 #endif

platforms/opencl/src/OpenCLKernels.cpp

@@ -4798,12 +4798,12 @@ void OpenCLCalcCustomHbondForceKernel::initialize(const System& system, const Cu
     for (int i = 0; i < (int) donorParams->getBuffers().size(); i++) {
         const OpenCLNonbondedUtilities::ParameterInfo& buffer = donorParams->getBuffers()[i];
         extraArgs << ", __global const "+buffer.getType()+"* restrict donor"+buffer.getName();
-        addDonorAndAcceptorCode(computeDonor, computeAcceptor, buffer.getType()+" donorParams"+cl.intToString(i+1)+" = donor"+buffer.getName()+"[index];\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, buffer.getType()+" donorParams"+cl.intToString(i+1)+" = donor"+buffer.getName()+"[donorIndex];\n");
     }
     for (int i = 0; i < (int) acceptorParams->getBuffers().size(); i++) {
         const OpenCLNonbondedUtilities::ParameterInfo& buffer = acceptorParams->getBuffers()[i];
         extraArgs << ", __global const "+buffer.getType()+"* restrict acceptor"+buffer.getName();
-        addDonorAndAcceptorCode(computeDonor, computeAcceptor, buffer.getType()+" acceptorParams"+cl.intToString(i+1)+" = acceptor"+buffer.getName()+"[index];\n");
+        addDonorAndAcceptorCode(computeDonor, computeAcceptor, buffer.getType()+" acceptorParams"+cl.intToString(i+1)+" = acceptor"+buffer.getName()+"[acceptorIndex];\n");
     }
 
     // Now evaluate the expressions.

platforms/opencl/src/kernels/customHbondForce.cl

@@ -91,8 +91,8 @@ __kernel void computeDonorForces(__global real4* restrict forceBuffers, __global
             barrier(CLK_LOCAL_MEM_FENCE);
             if (donorIndex < NUM_DONORS) {
                 for (int index = 0; index < blockSize; index++) {
-#ifdef USE_EXCLUSIONS
                     int acceptorIndex = acceptorStart+index;
+#ifdef USE_EXCLUSIONS
                     if (acceptorIndex == exclusionIndices.x || acceptorIndex == exclusionIndices.y || acceptorIndex == exclusionIndices.z || acceptorIndex == exclusionIndices.w)
                         continue;
 #endif
@@ -179,8 +179,8 @@ __kernel void computeAcceptorForces(__global real4* restrict forceBuffers, __glo
             barrier(CLK_LOCAL_MEM_FENCE);
             if (acceptorIndex < NUM_ACCEPTORS) {
                 for (int index = 0; index < blockSize; index++) {
-#ifdef USE_EXCLUSIONS
                     int donorIndex = donorStart+index;
+#ifdef USE_EXCLUSIONS
                     if (donorIndex == exclusionIndices.x || donorIndex == exclusionIndices.y || donorIndex == exclusionIndices.z || donorIndex == exclusionIndices.w)
                         continue;
 #endif

tests/TestCustomHbondForce.h

@@ -244,6 +244,33 @@ void testIllegalVariable() {
     ASSERT(threwException);
 }
 
+void testParameters() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomHbondForce* custom = new CustomHbondForce("(2*d+a)*distance(d1,a1)");
+    custom->addPerDonorParameter("d");
+    custom->addPerAcceptorParameter("a");
+    custom->addDonor(1, 0, -1, vector<double>({1.5}));
+    custom->addDonor(2, 0, -1, vector<double>({1.8}));
+    custom->addAcceptor(0, 1, -1, vector<double>({2.1}));
+    system.addForce(custom);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(0, 2, 0);
+    positions[2] = Vec3(2, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    ASSERT_EQUAL_VEC(Vec3((2*1.8+2.1), (2*1.5+2.1), 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, -(2*1.5+2.1), 0), forces[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-(2*1.8+2.1), 0, 0), forces[2], TOL);
+    ASSERT_EQUAL_TOL(2*(2*1.8+2.1)+2*(2*1.5+2.1), state.getPotentialEnergy(), TOL);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -254,6 +281,7 @@ int main(int argc, char* argv[]) {
         testCutoff();
         testCustomFunctions();
         testIllegalVariable();
+        testParameters();
         runPlatformTests();
     }
     catch(const exception& e) {