Implemented Lennard-Jones switching function for CUDA and OpenCL

platforms/cuda/src/CudaKernels.cpp

@@ -1438,6 +1438,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
     map<string, string> defines;
     defines["HAS_COULOMB"] = (hasCoulomb ? "1" : "0");
     defines["HAS_LENNARD_JONES"] = (hasLJ ? "1" : "0");
+    defines["USE_LJ_SWITCH"] = (useCutoff && force.getUseSwitchingFunction() ? "1" : "0");
     if (useCutoff) {
         // Compute the reaction field constants.
 
@@ -1445,6 +1446,15 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
         double reactionFieldC = (1.0 / force.getCutoffDistance())*(3.0*force.getReactionFieldDielectric())/(2.0*force.getReactionFieldDielectric()+1.0);
         defines["REACTION_FIELD_K"] = cu.doubleToString(reactionFieldK);
         defines["REACTION_FIELD_C"] = cu.doubleToString(reactionFieldC);
+        
+        // Compute the switching coefficients.
+        
+        if (force.getUseSwitchingFunction()) {
+            defines["LJ_SWITCH_CUTOFF"] = cu.doubleToString(force.getSwitchingDistance());
+            defines["LJ_SWITCH_C3"] = cu.doubleToString(10/pow(force.getSwitchingDistance()-force.getCutoffDistance(), 3.0));
+            defines["LJ_SWITCH_C4"] = cu.doubleToString(15/pow(force.getSwitchingDistance()-force.getCutoffDistance(), 4.0));
+            defines["LJ_SWITCH_C5"] = cu.doubleToString(6/pow(force.getSwitchingDistance()-force.getCutoffDistance(), 5.0));
+        }
     }
     if (force.getUseDispersionCorrection() && cu.getContextIndex() == 0)
         dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(system, force);

platforms/cuda/src/kernels/coulombLennardJones.cu

@@ -33,8 +33,19 @@ if (!isExcluded || needCorrection) {
             sig2 *= sig2;
             real sig6 = sig2*sig2*sig2;
             real epssig6 = sig6*(sigmaEpsilon1.y*sigmaEpsilon2.y);
-            tempForce = epssig6*(12.0f*sig6 - 6.0f) + prefactor*(erfcAlphaR+alphaR*expAlphaRSqr*TWO_OVER_SQRT_PI);
-            tempEnergy += epssig6*(sig6 - 1.0f) + prefactor*erfcAlphaR;
+            tempForce = epssig6*(12.0f*sig6 - 6.0f);
+            real ljEnergy = epssig6*(sig6 - 1.0f);
+            #if USE_LJ_SWITCH
+            if (r > LJ_SWITCH_CUTOFF) {
+                real x = r-LJ_SWITCH_CUTOFF;
+                real switchValue = 1+x*x*x*(LJ_SWITCH_C3+x*(LJ_SWITCH_C4+x*LJ_SWITCH_C5));
+                real switchDeriv = x*x*(3*LJ_SWITCH_C3+x*(4*LJ_SWITCH_C4+x*5*LJ_SWITCH_C5));
+                tempForce = tempForce*switchValue - ljEnergy*switchDeriv*r;
+                ljEnergy *= switchValue;
+            }
+            #endif
+            tempForce += prefactor*(erfcAlphaR+alphaR*expAlphaRSqr*TWO_OVER_SQRT_PI);
+            tempEnergy += ljEnergy + prefactor*erfcAlphaR;
 #else
             tempForce = prefactor*(erfcAlphaR+alphaR*expAlphaRSqr*TWO_OVER_SQRT_PI);
             tempEnergy += prefactor*erfcAlphaR;
@@ -58,7 +69,17 @@ if (!isExcluded || needCorrection) {
     real sig6 = sig2*sig2*sig2;
     real epssig6 = sig6*(sigmaEpsilon1.y*sigmaEpsilon2.y);
     tempForce = epssig6*(12.0f*sig6 - 6.0f);
-    tempEnergy += includeInteraction ? epssig6*(sig6 - 1) : 0;
+    real ljEnergy = includeInteraction ? epssig6*(sig6 - 1) : 0;
+    #if USE_LJ_SWITCH
+    if (r > LJ_SWITCH_CUTOFF) {
+        real x = r-LJ_SWITCH_CUTOFF;
+        real switchValue = 1+x*x*x*(LJ_SWITCH_C3+x*(LJ_SWITCH_C4+x*LJ_SWITCH_C5));
+        real switchDeriv = x*x*(3*LJ_SWITCH_C3+x*(4*LJ_SWITCH_C4+x*5*LJ_SWITCH_C5));
+        tempForce = tempForce*switchValue - ljEnergy*switchDeriv*r;
+        ljEnergy *= switchValue;
+    }
+    #endif
+    tempEnergy += ljEnergy;
   #endif
 #if HAS_COULOMB
   #ifdef USE_CUTOFF

platforms/cuda/tests/TestCudaNonbondedForce.cpp

@@ -6,7 +6,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:                                                              *
  *                                                                            *
@@ -813,6 +813,61 @@ void testParallelComputation(bool useCutoff) {
         ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
 }
 
+void testSwitchingFunction(NonbondedForce::NonbondedMethod method) {
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->addParticle(0, 1.2, 1);
+    nonbonded->addParticle(0, 1.4, 2);
+    nonbonded->setNonbondedMethod(method);
+    nonbonded->setCutoffDistance(2.0);
+    nonbonded->setUseSwitchingFunction(true);
+    nonbonded->setSwitchingDistance(1.5);
+    nonbonded->setUseDispersionCorrection(false);
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    double eps = SQRT_TWO;
+    
+    // Compute the interaction at various distances.
+    
+    for (double r = 1.0; r < 2.5; r += 0.1) {
+        positions[1] = Vec3(r, 0, 0);
+        context.setPositions(positions);
+        State state = context.getState(State::Forces | State::Energy);
+        
+        // See if the energy is correct.
+        
+        double x = 1.3/r;
+        double expectedEnergy = 4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0));
+        double switchValue;
+        if (r <= 1.5)
+            switchValue = 1;
+        else if (r >= 2.0)
+            switchValue = 0;
+        else {
+            double t = (r-1.5)/0.5;
+            switchValue = 1+t*t*t*(-10+t*(15-t*6));
+        }
+        ASSERT_EQUAL_TOL(switchValue*expectedEnergy, state.getPotentialEnergy(), TOL);
+        
+        // See if the force is the gradient of the energy.
+        
+        double delta = 1e-3;
+        positions[1] = Vec3(r-delta, 0, 0);
+        context.setPositions(positions);
+        double e1 = context.getState(State::Energy).getPotentialEnergy();
+        positions[1] = Vec3(r+delta, 0, 0);
+        context.setPositions(positions);
+        double e2 = context.getState(State::Energy).getPotentialEnergy();
+        ASSERT_EQUAL_TOL((e2-e1)/(2*delta), state.getForces()[0][0], 1e-3);
+    }
+}
+
 int main(int argc, char* argv[]) {
     try {
         if (argc > 1)
@@ -830,6 +885,8 @@ int main(int argc, char* argv[]) {
         testChangingParameters();
         testParallelComputation(false);
         testParallelComputation(true);
+        testSwitchingFunction(NonbondedForce::CutoffNonPeriodic);
+        testSwitchingFunction(NonbondedForce::PME);
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/opencl/src/OpenCLKernels.cpp

@@ -1411,6 +1411,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
     map<string, string> defines;
     defines["HAS_COULOMB"] = (hasCoulomb ? "1" : "0");
     defines["HAS_LENNARD_JONES"] = (hasLJ ? "1" : "0");
+    defines["USE_LJ_SWITCH"] = (useCutoff && force.getUseSwitchingFunction() ? "1" : "0");
     if (useCutoff) {
         // Compute the reaction field constants.
 
@@ -1418,6 +1419,15 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
         double reactionFieldC = (1.0 / force.getCutoffDistance())*(3.0*force.getReactionFieldDielectric())/(2.0*force.getReactionFieldDielectric()+1.0);
         defines["REACTION_FIELD_K"] = cl.doubleToString(reactionFieldK);
         defines["REACTION_FIELD_C"] = cl.doubleToString(reactionFieldC);
+        
+        // Compute the switching coefficients.
+        
+        if (force.getUseSwitchingFunction()) {
+            defines["LJ_SWITCH_CUTOFF"] = cl.doubleToString(force.getSwitchingDistance());
+            defines["LJ_SWITCH_C3"] = cl.doubleToString(10/pow(force.getSwitchingDistance()-force.getCutoffDistance(), 3.0));
+            defines["LJ_SWITCH_C4"] = cl.doubleToString(15/pow(force.getSwitchingDistance()-force.getCutoffDistance(), 4.0));
+            defines["LJ_SWITCH_C5"] = cl.doubleToString(6/pow(force.getSwitchingDistance()-force.getCutoffDistance(), 5.0));
+        }
     }
     if (force.getUseDispersionCorrection() && cl.getContextIndex() == 0)
         dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(system, force);

platforms/opencl/src/kernels/coulombLennardJones.cl

@@ -33,8 +33,19 @@ if (!isExcluded || needCorrection) {
             sig2 *= sig2;
             real sig6 = sig2*sig2*sig2;
             real epssig6 = sig6*(sigmaEpsilon1.y*sigmaEpsilon2.y);
-            tempForce = epssig6*(12.0f*sig6 - 6.0f) + prefactor*(erfcAlphaR+alphaR*expAlphaRSqr*TWO_OVER_SQRT_PI);
-            tempEnergy += epssig6*(sig6 - 1.0f) + prefactor*erfcAlphaR;
+            tempForce = epssig6*(12.0f*sig6 - 6.0f);
+            real ljEnergy = epssig6*(sig6 - 1.0f);
+            #if USE_LJ_SWITCH
+            if (r > LJ_SWITCH_CUTOFF) {
+                real x = r-LJ_SWITCH_CUTOFF;
+                real switchValue = 1+x*x*x*(LJ_SWITCH_C3+x*(LJ_SWITCH_C4+x*LJ_SWITCH_C5));
+                real switchDeriv = x*x*(3*LJ_SWITCH_C3+x*(4*LJ_SWITCH_C4+x*5*LJ_SWITCH_C5));
+                tempForce = tempForce*switchValue - ljEnergy*switchDeriv*r;
+                ljEnergy *= switchValue;
+            }
+            #endif
+            tempForce += prefactor*(erfcAlphaR+alphaR*expAlphaRSqr*TWO_OVER_SQRT_PI);
+            tempEnergy += ljEnergy + prefactor*erfcAlphaR;
 #else
             tempForce = prefactor*(erfcAlphaR+alphaR*expAlphaRSqr*TWO_OVER_SQRT_PI);
             tempEnergy += prefactor*erfcAlphaR;
@@ -63,7 +74,17 @@ if (!isExcluded || needCorrection) {
     real sig6 = sig2*sig2*sig2;
     real epssig6 = sig6*(sigmaEpsilon1.y*sigmaEpsilon2.y);
     tempForce = epssig6*(12.0f*sig6 - 6.0f);
-    tempEnergy += select((real) 0, epssig6*(sig6-1), includeInteraction);
+    real ljEnergy = select((real) 0, epssig6*(sig6-1), includeInteraction);
+    #if USE_LJ_SWITCH
+    if (r > LJ_SWITCH_CUTOFF) {
+        real x = r-LJ_SWITCH_CUTOFF;
+        real switchValue = 1+x*x*x*(LJ_SWITCH_C3+x*(LJ_SWITCH_C4+x*LJ_SWITCH_C5));
+        real switchDeriv = x*x*(3*LJ_SWITCH_C3+x*(4*LJ_SWITCH_C4+x*5*LJ_SWITCH_C5));
+        tempForce = tempForce*switchValue - ljEnergy*switchDeriv*r;
+        ljEnergy *= switchValue;
+    }
+    #endif
+    tempEnergy += ljEnergy;
   #endif
 #if HAS_COULOMB
   #ifdef USE_CUTOFF

platforms/opencl/tests/TestOpenCLNonbondedForce.cpp

@@ -6,7 +6,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:                                                              *
  *                                                                            *
@@ -816,6 +816,61 @@ void testParallelComputation(bool useCutoff) {
         ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
 }
 
+void testSwitchingFunction(NonbondedForce::NonbondedMethod method) {
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->addParticle(0, 1.2, 1);
+    nonbonded->addParticle(0, 1.4, 2);
+    nonbonded->setNonbondedMethod(method);
+    nonbonded->setCutoffDistance(2.0);
+    nonbonded->setUseSwitchingFunction(true);
+    nonbonded->setSwitchingDistance(1.5);
+    nonbonded->setUseDispersionCorrection(false);
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    double eps = SQRT_TWO;
+    
+    // Compute the interaction at various distances.
+    
+    for (double r = 1.0; r < 2.5; r += 0.1) {
+        positions[1] = Vec3(r, 0, 0);
+        context.setPositions(positions);
+        State state = context.getState(State::Forces | State::Energy);
+        
+        // See if the energy is correct.
+        
+        double x = 1.3/r;
+        double expectedEnergy = 4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0));
+        double switchValue;
+        if (r <= 1.5)
+            switchValue = 1;
+        else if (r >= 2.0)
+            switchValue = 0;
+        else {
+            double t = (r-1.5)/0.5;
+            switchValue = 1+t*t*t*(-10+t*(15-t*6));
+        }
+        ASSERT_EQUAL_TOL(switchValue*expectedEnergy, state.getPotentialEnergy(), TOL);
+        
+        // See if the force is the gradient of the energy.
+        
+        double delta = 1e-3;
+        positions[1] = Vec3(r-delta, 0, 0);
+        context.setPositions(positions);
+        double e1 = context.getState(State::Energy).getPotentialEnergy();
+        positions[1] = Vec3(r+delta, 0, 0);
+        context.setPositions(positions);
+        double e2 = context.getState(State::Energy).getPotentialEnergy();
+        ASSERT_EQUAL_TOL((e2-e1)/(2*delta), state.getForces()[0][0], 1e-3);
+    }
+}
+
 int main(int argc, char* argv[]) {
     try {
         if (argc > 1)
@@ -833,6 +888,8 @@ int main(int argc, char* argv[]) {
         testChangingParameters();
         testParallelComputation(false);
         testParallelComputation(true);
+        testSwitchingFunction(NonbondedForce::CutoffNonPeriodic);
+        testSwitchingFunction(NonbondedForce::PME);
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;