Implemented long range correction for CustomNonbondedForce in OpenCL

platforms/cuda/src/CudaKernels.h

@@ -9,7 +9,7 @@
  * 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.      *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *

platforms/opencl/src/OpenCLKernels.cpp

@@ -33,6 +33,7 @@
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/internal/CustomCompoundBondForceImpl.h"
 #include "openmm/internal/CustomHbondForceImpl.h"
+#include "openmm/internal/CustomNonbondedForceImpl.h"
 #include "openmm/internal/NonbondedForceImpl.h"
 #include "OpenCLBondedUtilities.h"
 #include "OpenCLExpressionUtilities.h"
@@ -1830,6 +1831,8 @@ OpenCLCalcCustomNonbondedForceKernel::~OpenCLCalcCustomNonbondedForceKernel() {
         delete tabulatedFunctionParams;
     for (int i = 0; i < (int) tabulatedFunctions.size(); i++)
         delete tabulatedFunctions[i];
+    if (forceCopy != NULL)
+        delete forceCopy;
 }
 
 void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, const CustomNonbondedForce& force) {
@@ -1938,6 +1941,17 @@ void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, cons
         cl.getNonbondedUtilities().addArgument(OpenCLNonbondedUtilities::ParameterInfo(prefix+"globals", "float", 1, sizeof(cl_float), globals->getDeviceBuffer()));
     }
     cl.addForce(new OpenCLCustomNonbondedForceInfo(cl.getNonbondedUtilities().getNumForceBuffers(), force));
+    
+    // Record information for the long range correction.
+    
+    if (force.getNonbondedMethod() == CustomNonbondedForce::CutoffPeriodic && force.getUseLongRangeCorrection() && cl.getContextIndex() == 0) {
+        forceCopy = new CustomNonbondedForce(force);
+        hasInitializedLongRangeCorrection = false;
+    }
+    else {
+        longRangeCoefficient = 0.0;
+        hasInitializedLongRangeCorrection = true;
+    }
 }
 
 double OpenCLCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -1949,10 +1963,20 @@ double OpenCLCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool
                 changed = true;
             globalParamValues[i] = value;
         }
-        if (changed)
+        if (changed) {
             globals->upload(globalParamValues);
+            if (forceCopy != NULL) {
+                longRangeCoefficient = CustomNonbondedForceImpl::calcLongRangeCorrection(*forceCopy, context.getOwner());
+                hasInitializedLongRangeCorrection = true;
+            }
+        }
     }
-    return 0.0;
+    if (!hasInitializedLongRangeCorrection) {
+        longRangeCoefficient = CustomNonbondedForceImpl::calcLongRangeCorrection(*forceCopy, context.getOwner());
+        hasInitializedLongRangeCorrection = true;
+    }
+    mm_double4 boxSize = cl.getPeriodicBoxSizeDouble();
+    return longRangeCoefficient/(boxSize.x*boxSize.y*boxSize.z);
 }
 
 void OpenCLCalcCustomNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force) {
@@ -1972,6 +1996,14 @@ void OpenCLCalcCustomNonbondedForceKernel::copyParametersToContext(ContextImpl&
     }
     params->setParameterValues(paramVector);
     
+    // If necessary, recompute the long range correction.
+    
+    if (forceCopy != NULL) {
+        longRangeCoefficient = CustomNonbondedForceImpl::calcLongRangeCorrection(force, context.getOwner());
+        hasInitializedLongRangeCorrection = true;
+        *forceCopy = force;
+    }
+    
     // Mark that the current reordering may be invalid.
     
     cl.invalidateMolecules();

platforms/opencl/src/OpenCLKernels.h

@@ -9,7 +9,7 @@
  * 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.      *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -634,7 +634,7 @@ private:
 class OpenCLCalcCustomNonbondedForceKernel : public CalcCustomNonbondedForceKernel {
 public:
     OpenCLCalcCustomNonbondedForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomNonbondedForceKernel(name, platform),
-            cl(cl), params(NULL), globals(NULL), tabulatedFunctionParams(NULL), system(system) {
+            cl(cl), params(NULL), globals(NULL), tabulatedFunctionParams(NULL), forceCopy(NULL), system(system) {
     }
     ~OpenCLCalcCustomNonbondedForceKernel();
     /**
@@ -668,6 +668,9 @@ private:
     std::vector<std::string> globalParamNames;
     std::vector<cl_float> globalParamValues;
     std::vector<OpenCLArray*> tabulatedFunctions;
+    double longRangeCoefficient;
+    bool hasInitializedLongRangeCorrection;
+    CustomNonbondedForce* forceCopy;
     const System& system;
 };
 

platforms/opencl/tests/TestOpenCLCustomNonbondedForce.cpp

@@ -7,7 +7,7 @@
  * 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.      *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -407,6 +407,81 @@ void testParallelComputation() {
         ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
 }
 
+void testLongRangeCorrection() {
+    // Create a box of particles.
+
+    int gridSize = 5;
+    int numParticles = gridSize*gridSize*gridSize;
+    double boxSize = gridSize*0.7;
+    double cutoff = boxSize/3;
+    System standardSystem;
+    System customSystem;
+    VerletIntegrator integrator1(0.01);
+    VerletIntegrator integrator2(0.01);
+    NonbondedForce* standardNonbonded = new NonbondedForce();
+    CustomNonbondedForce* customNonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6); sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)");
+    customNonbonded->addPerParticleParameter("sigma");
+    customNonbonded->addPerParticleParameter("eps");
+    vector<Vec3> positions(numParticles);
+    int index = 0;
+    vector<double> params1(2);
+    params1[0] = 1.1;
+    params1[1] = 0.5;
+    vector<double> params2(2);
+    params2[0] = 1;
+    params2[1] = 1;
+    for (int i = 0; i < gridSize; i++)
+        for (int j = 0; j < gridSize; j++)
+            for (int k = 0; k < gridSize; k++) {
+                standardSystem.addParticle(1.0);
+                customSystem.addParticle(1.0);
+                if (index%2 == 0) {
+                    standardNonbonded->addParticle(0, params1[0], params1[1]);
+                    customNonbonded->addParticle(params1);
+                }
+                else {
+                    standardNonbonded->addParticle(0, params2[0], params2[1]);
+                    customNonbonded->addParticle(params2);
+                }
+                positions[index] = Vec3(i*boxSize/gridSize, j*boxSize/gridSize, k*boxSize/gridSize);
+                index++;
+            }
+    standardNonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    customNonbonded->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
+    standardNonbonded->setCutoffDistance(cutoff);
+    customNonbonded->setCutoffDistance(cutoff);
+    standardSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    customSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    standardNonbonded->setUseDispersionCorrection(true);
+    customNonbonded->setUseLongRangeCorrection(true);
+    standardSystem.addForce(standardNonbonded);
+    customSystem.addForce(customNonbonded);
+
+    // Compute the correction for the standard force.
+
+    Context context1(standardSystem, integrator1, platform);
+    context1.setPositions(positions);
+    double standardEnergy1 = context1.getState(State::Energy).getPotentialEnergy();
+    standardNonbonded->setUseDispersionCorrection(false);
+    context1.reinitialize();
+    context1.setPositions(positions);
+    double standardEnergy2 = context1.getState(State::Energy).getPotentialEnergy();
+
+    // Compute the correction for the custom force.
+
+    Context context2(customSystem, integrator2, platform);
+    context2.setPositions(positions);
+    double customEnergy1 = context2.getState(State::Energy).getPotentialEnergy();
+    customNonbonded->setUseLongRangeCorrection(false);
+    context2.reinitialize();
+    context2.setPositions(positions);
+    double customEnergy2 = context2.getState(State::Energy).getPotentialEnergy();
+
+    // See if they agree.
+
+    ASSERT_EQUAL_TOL(standardEnergy1-standardEnergy2, customEnergy1-customEnergy2, 1e-4);
+}
+
 int main(int argc, char* argv[]) {
     try {
         if (argc > 1)
@@ -420,6 +495,7 @@ int main(int argc, char* argv[]) {
         testTabulatedFunction();
         testCoulombLennardJones();
         testParallelComputation();
+        testLongRangeCorrection();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;