Implemented switching function for CustomNonbondedForce

platforms/cuda/src/CudaKernels.cpp

@@ -1929,6 +1929,15 @@ void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const
     compute << cu.getExpressionUtilities().createExpressions(forceExpressions, variables, functionDefinitions, prefix+"temp", prefix+"functionParams");
     map<string, string> replacements;
     replacements["COMPUTE_FORCE"] = compute.str();
+    replacements["USE_SWITCH"] = (useCutoff && force.getUseSwitchingFunction() ? "1" : "0");
+    if (force.getUseSwitchingFunction()) {
+        // Compute the switching coefficients.
+        
+        replacements["SWITCH_CUTOFF"] = cu.doubleToString(force.getSwitchingDistance());
+        replacements["SWITCH_C3"] = cu.doubleToString(10/pow(force.getSwitchingDistance()-force.getCutoffDistance(), 3.0));
+        replacements["SWITCH_C4"] = cu.doubleToString(15/pow(force.getSwitchingDistance()-force.getCutoffDistance(), 4.0));
+        replacements["SWITCH_C5"] = cu.doubleToString(6/pow(force.getSwitchingDistance()-force.getCutoffDistance(), 5.0));
+    }
     string source = cu.replaceStrings(CudaKernelSources::customNonbonded, replacements);
     cu.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source, force.getForceGroup());
     for (int i = 0; i < (int) params->getBuffers().size(); i++) {

platforms/cuda/src/kernels/customNonbonded.cu

@@ -5,5 +5,14 @@ if (!isExcluded) {
 #endif
     real tempForce = 0;
     COMPUTE_FORCE
+#if USE_SWITCH
+    if (r > SWITCH_CUTOFF) {
+        real x = r-SWITCH_CUTOFF;
+        real switchValue = 1+x*x*x*(SWITCH_C3+x*(SWITCH_C4+x*SWITCH_C5));
+        real switchDeriv = x*x*(3*SWITCH_C3+x*(4*SWITCH_C4+x*5*SWITCH_C5));
+        tempForce = tempForce*switchValue - tempEnergy*switchDeriv;
+        tempEnergy *= switchValue;
+    }
+#endif
     dEdR += tempForce*invR;
 }

platforms/cuda/tests/TestCudaCustomNonbondedForce.cpp

@@ -407,6 +407,58 @@ void testParallelComputation() {
         ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
 }
 
+void testSwitchingFunction() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomNonbondedForce* nonbonded = new CustomNonbondedForce("10/r^2");
+    vector<double> params;
+    nonbonded->addParticle(params);
+    nonbonded->addParticle(params);
+    nonbonded->setNonbondedMethod(CustomNonbondedForce::CutoffNonPeriodic);
+    nonbonded->setCutoffDistance(2.0);
+    nonbonded->setUseSwitchingFunction(true);
+    nonbonded->setSwitchingDistance(1.5);
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    
+    // 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 expectedEnergy = 10/(r*r);
+        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);
+    }
+}
+
 void testLongRangeCorrection() {
     // Create a box of particles.
 
@@ -495,6 +547,7 @@ int main(int argc, char* argv[]) {
         testTabulatedFunction();
         testCoulombLennardJones();
         testParallelComputation();
+        testSwitchingFunction();
         testLongRangeCorrection();
     }
     catch(const exception& e) {

platforms/opencl/src/OpenCLKernels.cpp

@@ -1940,6 +1940,15 @@ void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, cons
     compute << cl.getExpressionUtilities().createExpressions(forceExpressions, variables, functionDefinitions, prefix+"temp", prefix+"functionParams");
     map<string, string> replacements;
     replacements["COMPUTE_FORCE"] = compute.str();
+    replacements["USE_SWITCH"] = (useCutoff && force.getUseSwitchingFunction() ? "1" : "0");
+    if (force.getUseSwitchingFunction()) {
+        // Compute the switching coefficients.
+        
+        replacements["SWITCH_CUTOFF"] = cl.doubleToString(force.getSwitchingDistance());
+        replacements["SWITCH_C3"] = cl.doubleToString(10/pow(force.getSwitchingDistance()-force.getCutoffDistance(), 3.0));
+        replacements["SWITCH_C4"] = cl.doubleToString(15/pow(force.getSwitchingDistance()-force.getCutoffDistance(), 4.0));
+        replacements["SWITCH_C5"] = cl.doubleToString(6/pow(force.getSwitchingDistance()-force.getCutoffDistance(), 5.0));
+    }
     string source = cl.replaceStrings(OpenCLKernelSources::customNonbonded, replacements);
     cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source, force.getForceGroup());
     for (int i = 0; i < (int) params->getBuffers().size(); i++) {

platforms/opencl/src/kernels/customNonbonded.cl

@@ -5,5 +5,14 @@ if (!isExcluded) {
 #endif
     real tempForce = 0.0f;
     COMPUTE_FORCE
+#if USE_SWITCH
+    if (r > SWITCH_CUTOFF) {
+        real x = r-SWITCH_CUTOFF;
+        real switchValue = 1+x*x*x*(SWITCH_C3+x*(SWITCH_C4+x*SWITCH_C5));
+        real switchDeriv = x*x*(3*SWITCH_C3+x*(4*SWITCH_C4+x*5*SWITCH_C5));
+        tempForce = tempForce*switchValue - tempEnergy*switchDeriv;
+        tempEnergy *= switchValue;
+    }
+#endif
     dEdR += tempForce*invR;
 }

platforms/opencl/tests/TestOpenCLCustomNonbondedForce.cpp

@@ -407,6 +407,58 @@ void testParallelComputation() {
         ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
 }
 
+void testSwitchingFunction() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomNonbondedForce* nonbonded = new CustomNonbondedForce("10/r^2");
+    vector<double> params;
+    nonbonded->addParticle(params);
+    nonbonded->addParticle(params);
+    nonbonded->setNonbondedMethod(CustomNonbondedForce::CutoffNonPeriodic);
+    nonbonded->setCutoffDistance(2.0);
+    nonbonded->setUseSwitchingFunction(true);
+    nonbonded->setSwitchingDistance(1.5);
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    
+    // 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 expectedEnergy = 10/(r*r);
+        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);
+    }
+}
+
 void testLongRangeCorrection() {
     // Create a box of particles.
 
@@ -495,6 +547,7 @@ int main(int argc, char* argv[]) {
         testTabulatedFunction();
         testCoulombLennardJones();
         testParallelComputation();
+        testSwitchingFunction();
         testLongRangeCorrection();
     }
     catch(const exception& e) {

platforms/reference/src/ReferenceKernels.cpp

@@ -1046,10 +1046,15 @@ void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, c
     }
     nonbondedMethod = CalcCustomNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
     nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
-    if (nonbondedMethod == NoCutoff)
+    if (nonbondedMethod == NoCutoff) {
         neighborList = NULL;
-    else
+        useSwitchingFunction = false;
+    }
+    else {
         neighborList = new NeighborList();
+        useSwitchingFunction = force.getUseSwitchingFunction();
+        switchingDistance = force.getSwitchingDistance();
+    }
 
     // Create custom functions for the tabulated functions.
 
@@ -1115,6 +1120,8 @@ double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bo
             globalParamsChanged = true;
         globalParamValues[globalParameterNames[i]] = value;
     }
+    if (useSwitchingFunction)
+        ixn.setUseSwitchingFunction(switchingDistance);
     ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, globalParamValues, forceData, 0, includeEnergy ? &energy : NULL);
     
     // Add in the long range correction.

platforms/reference/src/ReferenceKernels.h

@@ -617,8 +617,8 @@ private:
     int numParticles;
     int **exclusionArray;
     RealOpenMM **particleParamArray;
-    RealOpenMM nonbondedCutoff, periodicBoxSize[3], longRangeCoefficient;
-    bool hasInitializedLongRangeCorrection;
+    RealOpenMM nonbondedCutoff, switchingDistance, periodicBoxSize[3], longRangeCoefficient;
+    bool useSwitchingFunction, hasInitializedLongRangeCorrection;
     CustomNonbondedForce* forceCopy;
     std::map<std::string, double> globalParamValues;
     std::vector<std::set<int> > exclusions;

platforms/reference/src/SimTKReference/ReferenceCustomNonbondedIxn.cpp

@@ -45,7 +45,7 @@ using OpenMM::RealVec;
 
 ReferenceCustomNonbondedIxn::ReferenceCustomNonbondedIxn(const Lepton::ExpressionProgram& energyExpression,
         const Lepton::ExpressionProgram& forceExpression, const vector<string>& parameterNames) :
-            cutoff(false), periodic(false), energyExpression(energyExpression), forceExpression(forceExpression), paramNames(parameterNames) {
+            cutoff(false), useSwitch(false), periodic(false), energyExpression(energyExpression), forceExpression(forceExpression), paramNames(parameterNames) {
 
    // ---------------------------------------------------------------------------------------
 
@@ -94,6 +94,19 @@ ReferenceCustomNonbondedIxn::~ReferenceCustomNonbondedIxn( ){
     neighborList = &neighbors;
   }
 
+/**---------------------------------------------------------------------------------------
+
+   Set the force to use a switching function.
+
+   @param distance            the switching distance
+
+   --------------------------------------------------------------------------------------- */
+
+void ReferenceCustomNonbondedIxn::setUseSwitchingFunction( RealOpenMM distance ) {
+    useSwitch = true;
+    switchingDistance = distance;
+}
+
   /**---------------------------------------------------------------------------------------
 
      Set the force to use periodic boundary conditions.  This requires that a cutoff has
@@ -228,13 +241,24 @@ void ReferenceCustomNonbondedIxn::calculateOneIxn( int ii, int jj, vector<RealVe
         ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR );
     else
         ReferenceForce::getDeltaR( atomCoordinates[jj], atomCoordinates[ii], deltaR );
-    if (cutoff && deltaR[ReferenceForce::RIndex] >= cutoffDistance)
+    RealOpenMM r = deltaR[ReferenceForce::RIndex];
+    if (cutoff && r >= cutoffDistance)
         return;
 
     // accumulate forces
 
-    variables["r"] = deltaR[ReferenceForce::RIndex];
+    variables["r"] = r;
     RealOpenMM dEdR = (RealOpenMM) (forceExpression.evaluate(variables)/(deltaR[ReferenceForce::RIndex]));
+    RealOpenMM energy = (RealOpenMM) energyExpression.evaluate(variables);
+    if (useSwitch) {
+        if (r > switchingDistance) {
+            RealOpenMM t = (r-switchingDistance)/(cutoffDistance-switchingDistance);
+            RealOpenMM switchValue = 1+t*t*t*(-10+t*(15-t*6));
+            RealOpenMM switchDeriv = t*t*(-30+t*(60-t*30))/(cutoffDistance-switchingDistance);
+            dEdR = switchValue*dEdR + energy*switchDeriv/r;
+            energy *= switchValue;
+        }
+    }
     for( int kk = 0; kk < 3; kk++ ){
        RealOpenMM force  = -dEdR*deltaR[kk];
        forces[ii][kk]   += force;
@@ -244,7 +268,6 @@ void ReferenceCustomNonbondedIxn::calculateOneIxn( int ii, int jj, vector<RealVe
     // accumulate energies
 
     if( totalEnergy || energyByAtom ) {
-        RealOpenMM energy = (RealOpenMM) energyExpression.evaluate(variables);
         if( totalEnergy )
            *totalEnergy += energy;
         if( energyByAtom ){

platforms/reference/src/SimTKReference/ReferenceCustomNonbondedIxn.h

@@ -38,10 +38,11 @@ class ReferenceCustomNonbondedIxn {
    private:
 
       bool cutoff;
+      bool useSwitch;
       bool periodic;
       const OpenMM::NeighborList* neighborList;
       RealOpenMM periodicBoxSize[3];
-      RealOpenMM cutoffDistance;
+      RealOpenMM cutoffDistance, switchingDistance;
       Lepton::ExpressionProgram energyExpression;
       Lepton::ExpressionProgram forceExpression;
       std::vector<std::string> paramNames;
@@ -96,6 +97,16 @@ class ReferenceCustomNonbondedIxn {
 
       void setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors );
 
+      /**---------------------------------------------------------------------------------------
+      
+         Set the force to use a switching function.
+      
+         @param distance            the switching distance
+      
+         --------------------------------------------------------------------------------------- */
+      
+      void setUseSwitchingFunction( RealOpenMM distance );
+
       /**---------------------------------------------------------------------------------------
 
          Set the force to use periodic boundary conditions.  This requires that a cutoff has

platforms/reference/tests/TestReferenceCustomNonbondedForce.cpp

@@ -337,6 +337,59 @@ void testCoulombLennardJones() {
     }
 }
 
+void testSwitchingFunction() {
+    ReferencePlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomNonbondedForce* nonbonded = new CustomNonbondedForce("10/r^2");
+    vector<double> params;
+    nonbonded->addParticle(params);
+    nonbonded->addParticle(params);
+    nonbonded->setNonbondedMethod(CustomNonbondedForce::CutoffNonPeriodic);
+    nonbonded->setCutoffDistance(2.0);
+    nonbonded->setUseSwitchingFunction(true);
+    nonbonded->setSwitchingDistance(1.5);
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    
+    // 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 expectedEnergy = 10/(r*r);
+        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);
+    }
+}
+
 void testLongRangeCorrection() {
     // Create a box of particles.
 
@@ -422,6 +475,7 @@ int main() {
         testPeriodic();
         testTabulatedFunction();
         testCoulombLennardJones();
+        testSwitchingFunction();
         testLongRangeCorrection();
     }
     catch(const exception& e) {