Long range correction for both NonbondedForce and CustomNonbondedForce takes...

Long range correction for both NonbondedForce and CustomNonbondedForce takes the switching function into account

openmmapi/include/openmm/CustomNonbondedForce.h

@@ -80,6 +80,25 @@ namespace OpenMM {
  * if the labels 1 and 2 are reversed.  In contrast, if it depended on the difference sigma1-sigma2, the results would
  * be undefined, because reversing the labels 1 and 2 would change the energy.
  * 
+ * When using a cutoff, by default the interaction is sharply truncated at the cutoff distance.
+ * Optionally you can instead use a switching function to make the interaction smoothly go to zero over a finite
+ * distance range.  To enable this, call setUseSwitchingFunction().  You must also call setSwitchingDistance()
+ * to specify the distance at which the interaction should begin to decrease.  The switching distance must be
+ * less than the cutoff distance.  Of course, you could also incorporate the switching function directly into your
+ * energy expression, but there are several advantages to keeping it separate.  It makes your energy expression simpler
+ * to write and understand.  It allows you to use the same energy expression with or without a cutoff.  Also, when using
+ * a long range correction (see below), separating out the switching function allows the correction to be calculated
+ * more accurately.
+ * 
+ * Another optional feature of this class is to add a contribution to the energy which approximates the effect of all
+ * interactions beyond the cutoff in a periodic system.  When running a simulation at constant pressure, this can improve
+ * the quality of the result.  Call setUseLongRangeCorrection() to enable it.
+ * 
+ * Computing the long range correction takes negligible work in each time step, but it does require an expensive precomputation
+ * at the start of the simulation.  Furthermore, that precomputation must be repeated every time a global parameter changes
+ * (or when you modify per-particle parameters by calling updateParametersInContext()).  This means that if parameters change
+ * frequently, the long range correction can be very slow.  For this reason, it is disabled by default.
+ * 
  * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
  * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, step, delta.  All trigonometric functions
  * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  delta(x) = 1 if x is 0, 0 otherwise.  The names of per-particle parameters

openmmapi/include/openmm/NonbondedForce.h

@@ -71,6 +71,11 @@ namespace OpenMM {
  * distance range.  To enable this, call setUseSwitchingFunction().  You must also call setSwitchingDistance()
  * to specify the distance at which the interaction should begin to decrease.  The switching distance must be
  * less than the cutoff distance.
+ * 
+ * Another optional feature of this class (enabled by default) is to add a contribution to the energy which approximates
+ * the effect of all Lennard-Jones interactions beyond the cutoff in a periodic system.  When running a simulation
+ * at constant pressure, this can improve the quality of the result.  Call setUseDispersionCorrection() to set whether
+ * this should be used.
  */
 
 class OPENMM_EXPORT NonbondedForce : public Force {

openmmapi/include/openmm/internal/NonbondedForceImpl.h

@@ -83,6 +83,7 @@ private:
     class ErrorFunction;
     class EwaldErrorFunction;
     static int findZero(const ErrorFunction& f, int initialGuess);
+    static double evalIntegral(double r, double rs, double rc, double sigma);
     const NonbondedForce& owner;
     Kernel kernel;
 };

openmmapi/src/CustomNonbondedForceImpl.cpp

@@ -240,13 +240,11 @@ double CustomNonbondedForceImpl::integrateInteraction(const Lepton::ExpressionPr
     // means we multiply the function by r^4.  Use the midpoint method.
 
     double cutoff = force.getCutoffDistance();
-    variables["r"] = 2*cutoff;
-    double sum = expression.evaluate(variables);
+    double sum = 0;
     int numPoints = 1;
-    for (int iteration = 0; iteration < 10; iteration++) {
+    for (int iteration = 0; ; iteration++) {
         double oldSum = sum;
         double newSum = 0;
-        numPoints *= 3;
         for (int i = 0; i < numPoints; i++) {
             if (i%3 == 1)
                 continue;
@@ -259,6 +257,38 @@ double CustomNonbondedForceImpl::integrateInteraction(const Lepton::ExpressionPr
         sum = newSum/numPoints + oldSum/3;
         if (iteration > 2 && (fabs((sum-oldSum)/sum) < 1e-5 || sum == 0))
             break;
+        if (iteration == 8)
+            throw OpenMMException("CustomNonbondedForce: Long range correction did not converge.  Does the energy go to 0 faster than 1/r^2?");
+        numPoints *= 3;
+    }
+    
+    // If a switching function is used, integrate over the switching interval.
+    
+    double sum2 = 0;
+    if (force.getUseSwitchingFunction()) {
+        double rswitch = force.getSwitchingDistance();
+        sum2 = 0;
+        numPoints = 1;
+        for (int iteration = 0; ; iteration++) {
+            double oldSum = sum2;
+            double newSum = 0;
+            for (int i = 0; i < numPoints; i++) {
+                if (i%3 == 1)
+                    continue;
+                double x = (i+0.5)/numPoints;
+                double r = rswitch+x*(cutoff-rswitch);
+                double switchValue = x*x*x*(10+x*(-15+x*6));
+                variables["r"] = r;
+                newSum += switchValue*expression.evaluate(variables)*r*r;
+            }
+            sum2 = newSum/numPoints + oldSum/3;
+            if (iteration > 2 && (fabs((sum2-oldSum)/sum2) < 1e-5 || sum2 == 0))
+                break;
+            if (iteration == 8)
+                throw OpenMMException("CustomNonbondedForce: Long range correction did not converge.  Is the energy finite everywhere in the switching interval?");
+            numPoints *= 3;
+        }
+        sum2 *= cutoff-rswitch;
     }
-    return sum/cutoff;
+    return sum/cutoff+sum2;
 }

openmmapi/src/NonbondedForceImpl.cpp

@@ -174,6 +174,43 @@ int NonbondedForceImpl::findZero(const NonbondedForceImpl::ErrorFunction& f, int
     return arg;
 }
 
+double NonbondedForceImpl::evalIntegral(double r, double rs, double rc, double sigma) {
+    // Compute the indefinite integral of the LJ interaction multiplied by the switching function.
+    // This is a large and somewhat horrifying expression, though it does grow on you if you look
+    // at it long enough.  Perhaps it could be simplified further, but I got tired of working on it.
+    
+    double A = 1/(rc-rs);
+    double A2 = A*A;
+    double A3 = A2*A;
+    double sig2 = sigma*sigma;
+    double sig6 = sig2*sig2*sig2;
+    double rs2 = rs*rs;
+    double rs3 = rs*rs2;
+    double r2 = r*r;
+    double r3 = r*r2;
+    double r4 = r*r3;
+    double r5 = r*r4;
+    double r6 = r*r5;
+    double r9 = r3*r6;
+    return sig6*A3*((
+        sig6*(
+            + rs3*28*(6*rs2*A2 + 15*rs*A + 10)
+            - r*rs2*945*(rs2*A2 + 2*rs*A + 1)
+            + r2*rs*1080*(2*rs2*A2 + 3*rs*A + 1)
+            - r3*420*(6*rs2*A2 + 6*rs*A + 1)
+            + r4*756*(2*rs*A2 + A)
+            - r5*378*A2)
+        -r6*(
+            + rs3*84*(6*rs2*A2 + 15*rs*A + 10)
+            - r*rs2*3780*(rs2*A2 + 2*rs*A + 1)
+            + r2*rs*7560*(2*rs2*A2 + 3*rs*A + 1))
+        )/(252*r9)
+     - log(r)*10*(6*rs2*A2 + 6*rs*A + 1)
+     + r*15*(2*rs*A2 + A)
+     - r2*3*A2
+    );
+}
+
 double NonbondedForceImpl::calcDispersionCorrection(const System& system, const NonbondedForce& force) {
     if (force.getNonbondedMethod() == NonbondedForce::NoCutoff || force.getNonbondedMethod() == NonbondedForce::CutoffNonPeriodic)
         return 0.0;
@@ -195,7 +232,10 @@ double NonbondedForceImpl::calcDispersionCorrection(const System& system, const
 
     // Loop over all pairs of classes to compute the coefficient.
 
-    double sum1 = 0, sum2 = 0;
+    double sum1 = 0, sum2 = 0, sum3 = 0;
+    bool useSwitch = force.getUseSwitchingFunction();
+    double cutoff = force.getCutoffDistance();
+    double switchDist = force.getSwitchingDistance();
     for (map<pair<double, double>, int>::const_iterator entry = classCounts.begin(); entry != classCounts.end(); ++entry) {
         double sigma = entry->first.first;
         double epsilon = entry->first.second;
@@ -204,6 +244,8 @@ double NonbondedForceImpl::calcDispersionCorrection(const System& system, const
         double sigma6 = sigma2*sigma2*sigma2;
         sum1 += count*epsilon*sigma6*sigma6;
         sum2 += count*epsilon*sigma6;
+        if (useSwitch)
+            sum3 += count*epsilon*(evalIntegral(cutoff, switchDist, cutoff, sigma)-evalIntegral(switchDist, switchDist, cutoff, sigma));
     }
     for (map<pair<double, double>, int>::const_iterator class1 = classCounts.begin(); class1 != classCounts.end(); ++class1)
         for (map<pair<double, double>, int>::const_iterator class2 = classCounts.begin(); class2 != class1; ++class2) {
@@ -214,13 +256,15 @@ double NonbondedForceImpl::calcDispersionCorrection(const System& system, const
             double sigma6 = sigma2*sigma2*sigma2;
             sum1 += count*epsilon*sigma6*sigma6;
             sum2 += count*epsilon*sigma6;
+            if (useSwitch)
+                sum3 += count*epsilon*(evalIntegral(cutoff, switchDist, cutoff, sigma)-evalIntegral(switchDist, switchDist, cutoff, sigma));
         }
     int numParticles = system.getNumParticles();
     int numInteractions = (numParticles*(numParticles+1))/2;
     sum1 /= numInteractions;
     sum2 /= numInteractions;
-    double cutoff = force.getCutoffDistance();
-    return 8*numParticles*numParticles*M_PI*(sum1/(9*pow(cutoff, 9))-sum2/(3*pow(cutoff, 3)));
+    sum3 /= numInteractions;
+    return 8*numParticles*numParticles*M_PI*(sum1/(9*pow(cutoff, 9))-sum2/(3*pow(cutoff, 3))+sum3);
 }
 
 void NonbondedForceImpl::updateParametersInContext(ContextImpl& context) {

platforms/cuda/tests/TestCudaCustomNonbondedForce.cpp

@@ -506,6 +506,10 @@ void testLongRangeCorrection() {
     customSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     standardNonbonded->setUseDispersionCorrection(true);
     customNonbonded->setUseLongRangeCorrection(true);
+    standardNonbonded->setUseSwitchingFunction(true);
+    customNonbonded->setUseSwitchingFunction(true);
+    standardNonbonded->setSwitchingDistance(0.8*cutoff);
+    customNonbonded->setSwitchingDistance(0.8*cutoff);
     standardSystem.addForce(standardNonbonded);
     customSystem.addForce(customNonbonded);
 

platforms/opencl/tests/TestOpenCLCustomNonbondedForce.cpp

@@ -506,6 +506,10 @@ void testLongRangeCorrection() {
     customSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     standardNonbonded->setUseDispersionCorrection(true);
     customNonbonded->setUseLongRangeCorrection(true);
+    standardNonbonded->setUseSwitchingFunction(true);
+    customNonbonded->setUseSwitchingFunction(true);
+    standardNonbonded->setSwitchingDistance(0.8*cutoff);
+    customNonbonded->setSwitchingDistance(0.8*cutoff);
     standardSystem.addForce(standardNonbonded);
     customSystem.addForce(customNonbonded);
 

platforms/reference/tests/TestReferenceCustomNonbondedForce.cpp

@@ -438,6 +438,10 @@ void testLongRangeCorrection() {
     customSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     standardNonbonded->setUseDispersionCorrection(true);
     customNonbonded->setUseLongRangeCorrection(true);
+    standardNonbonded->setUseSwitchingFunction(true);
+    customNonbonded->setUseSwitchingFunction(true);
+    standardNonbonded->setSwitchingDistance(0.8*cutoff);
+    customNonbonded->setSwitchingDistance(0.8*cutoff);
     standardSystem.addForce(standardNonbonded);
     customSystem.addForce(customNonbonded);