Began implementing long range correction for CustomNonbondedForce

openmmapi/include/openmm/Context.h

@@ -170,7 +170,7 @@ public:
      * 
      * @param name the name of the parameter to get
      */
-    double getParameter(const std::string& name);
+    double getParameter(const std::string& name) const;
     /**
      * Set the value of an adjustable parameter defined by a Force object in the System.
      * 

openmmapi/include/openmm/CustomNonbondedForce.h

@@ -178,6 +178,34 @@ public:
      * @param distance    the cutoff distance, measured in nm
      */
     void setCutoffDistance(double distance);
+    /**
+     * Get whether a switching is applied to the interaction.
+     */
+    bool getUseSwitchingFunction() const;
+    /**
+     * Set whether a switching is applied to the interaction.
+     */
+    void setUseSwitchingFunction(bool use);
+    /**
+     * Get the distance at which the switching function begins to reduce the interaction.  This must be
+     * less than the cutoff distance.
+     */
+    double getSwitchingDistance() const;
+    /**
+     * Set the distance at which the switching function begins to reduce the interaction.  This must be
+     * less than the cutoff distance.
+     */
+    void setSwitchingDistance(double distance);
+    /**
+     * Get whether to add a correction to the energy to compensate for the cutoff and switching function.
+     * This has no effect if periodic boundary conditions are not used.
+     */
+    bool getUseLongRangeCorrection() const;
+    /**
+     * Set whether to add a correction to the energy to compensate for the cutoff and switching function.
+     * This has no effect if periodic boundary conditions are not used.
+     */
+    void setUseLongRangeCorrection(bool use);
     /**
      * Add a new per-particle parameter that the interaction may depend on.
      *
@@ -335,7 +363,8 @@ private:
     class ExclusionInfo;
     class FunctionInfo;
     NonbondedMethod nonbondedMethod;
-    double cutoffDistance;
+    double cutoffDistance, switchingDistance;
+    bool useSwitchingFunction, useLongRangeCorrection;
     std::string energyExpression;
     std::vector<PerParticleParameterInfo> parameters;
     std::vector<GlobalParameterInfo> globalParameters;

openmmapi/include/openmm/internal/CustomNonbondedForceImpl.h

@@ -35,6 +35,7 @@
 #include "ForceImpl.h"
 #include "openmm/CustomNonbondedForce.h"
 #include "openmm/Kernel.h"
+#include "lepton/ExpressionProgram.h"
 #include <utility>
 #include <map>
 #include <string>
@@ -60,7 +61,15 @@ public:
     std::map<std::string, double> getDefaultParameters();
     std::vector<std::string> getKernelNames();
     void updateParametersInContext(ContextImpl& context);
+    /**
+     * Compute the coefficient which, when divided by the periodic box volume, gives the
+     * long range correction to the energy.
+     */
+    static double calcLongRangeCorrection(const CustomNonbondedForce& force, const Context& context);
 private:
+    class TabulatedFunction;
+    static double integrateInteraction(const Lepton::ExpressionProgram& expression, const std::vector<double>& params1, const std::vector<double>& params2,
+            const CustomNonbondedForce& force, const Context& context);
     const CustomNonbondedForce& owner;
     Kernel kernel;
 };

openmmapi/src/Context.cpp

@@ -234,7 +234,7 @@ void Context::setVelocitiesToTemperature(double temperature, int randomSeed) {
     impl->applyVelocityConstraints(1e-5);
 }
 
-double Context::getParameter(const string& name) {
+double Context::getParameter(const string& name) const {
     return impl->getParameter(name);
 }
 

openmmapi/src/CustomNonbondedForce.cpp

@@ -47,7 +47,8 @@ using std::string;
 using std::stringstream;
 using std::vector;
 
-CustomNonbondedForce::CustomNonbondedForce(const string& energy) : energyExpression(energy), nonbondedMethod(NoCutoff), cutoffDistance(1.0) {
+CustomNonbondedForce::CustomNonbondedForce(const string& energy) : energyExpression(energy), nonbondedMethod(NoCutoff), cutoffDistance(1.0),
+    switchingDistance(-1.0), useSwitchingFunction(false), useLongRangeCorrection(false) {
 }
 
 const string& CustomNonbondedForce::getEnergyFunction() const {
@@ -74,6 +75,30 @@ void CustomNonbondedForce::setCutoffDistance(double distance) {
     cutoffDistance = distance;
 }
 
+bool CustomNonbondedForce::getUseSwitchingFunction() const {
+    return useSwitchingFunction;
+}
+
+void CustomNonbondedForce::setUseSwitchingFunction(bool use) {
+    useSwitchingFunction = use;
+}
+
+double CustomNonbondedForce::getSwitchingDistance() const {
+    return switchingDistance;
+}
+
+void CustomNonbondedForce::setSwitchingDistance(double distance) {
+    switchingDistance = distance;
+}
+
+bool CustomNonbondedForce::getUseLongRangeCorrection() const {
+    return useLongRangeCorrection;
+}
+
+void CustomNonbondedForce::setUseLongRangeCorrection(bool use) {
+    useLongRangeCorrection = use;
+}
+
 int CustomNonbondedForce::addPerParticleParameter(const string& name) {
     parameters.push_back(PerParticleParameterInfo(name));
     return parameters.size()-1;

openmmapi/src/CustomNonbondedForceImpl.cpp

@@ -32,7 +32,12 @@
 #include "openmm/OpenMMException.h"
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/internal/CustomNonbondedForceImpl.h"
+#include "openmm/internal/SplineFitter.h"
 #include "openmm/kernels.h"
+#include "lepton/CustomFunction.h"
+#include "lepton/ParsedExpression.h"
+#include "lepton/Parser.h"
+#include <cmath>
 #include <sstream>
 
 using namespace OpenMM;
@@ -57,6 +62,10 @@ void CustomNonbondedForceImpl::initialize(ContextImpl& context) {
     const System& system = context.getSystem();
     if (owner.getNumParticles() != system.getNumParticles())
         throw OpenMMException("CustomNonbondedForce must have exactly as many particles as the System it belongs to.");
+    if (owner.getUseSwitchingFunction()) {
+        if (owner.getSwitchingDistance() < 0 || owner.getSwitchingDistance() >= owner.getCutoffDistance())
+            throw OpenMMException("CustomNonbondedForce: Switching distance must satisfy 0 <= r_switch < r_cutoff");
+    }
     vector<set<int> > exclusions(owner.getNumParticles());
     vector<double> parameters;
     int numParameters = owner.getNumPerParticleParameters();
@@ -127,3 +136,126 @@ map<string, double> CustomNonbondedForceImpl::getDefaultParameters() {
 void CustomNonbondedForceImpl::updateParametersInContext(ContextImpl& context) {
     kernel.getAs<CalcCustomNonbondedForceKernel>().copyParametersToContext(context, owner);
 }
+
+class CustomNonbondedForceImpl::TabulatedFunction : public Lepton::CustomFunction {
+public:
+    TabulatedFunction(double min, double max, const vector<double>& values) :
+            min(min), max(max), values(values) {
+        int numValues = values.size();
+        x.resize(numValues);
+        for (int i = 0; i < numValues; i++)
+            x[i] = min+i*(max-min)/(numValues-1);
+        SplineFitter::createNaturalSpline(x, values, derivs);
+    }
+    int getNumArguments() const {
+        return 1;
+    }
+    double evaluate(const double* arguments) const {
+        double t = arguments[0];
+        if (t < min || t > max)
+            return 0.0;
+        return SplineFitter::evaluateSpline(x, values, derivs, t);
+    }
+    double evaluateDerivative(const double* arguments, const int* derivOrder) const {
+        double t = arguments[0];
+        if (t < min || t > max)
+            return 0.0;
+        return SplineFitter::evaluateSplineDerivative(x, values, derivs, t);
+    }
+    CustomFunction* clone() const {
+        return new TabulatedFunction(min, max, values);
+    }
+    double min, max;
+    vector<double> x, values, derivs;
+};
+
+double CustomNonbondedForceImpl::calcLongRangeCorrection(const CustomNonbondedForce& force, const Context& context) {
+    if (force.getNonbondedMethod() == CustomNonbondedForce::NoCutoff || force.getNonbondedMethod() == CustomNonbondedForce::CutoffNonPeriodic)
+        return 0.0;
+    
+    // Identify all particle classes (defined by parameters), and count the number of
+    // particles in each class.
+
+    map<vector<double>, int> classCounts;
+    for (int i = 0; i < force.getNumParticles(); i++) {
+        vector<double> parameters;
+        force.getParticleParameters(i, parameters);
+        map<vector<double>, int>::iterator entry = classCounts.find(parameters);
+        if (entry == classCounts.end())
+            classCounts[parameters] = 1;
+        else
+            entry->second++;
+    }
+    
+    // Parse the energy expression.
+    
+    map<string, Lepton::CustomFunction*> functions;
+    for (int i = 0; i < force.getNumFunctions(); i++) {
+        string name;
+        vector<double> values;
+        double min, max;
+        force.getFunctionParameters(i, name, values, min, max);
+        functions[name] = new TabulatedFunction(min, max, values);
+    }
+    Lepton::ExpressionProgram expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize().createProgram();
+
+    // Loop over all pairs of classes to compute the coefficient.
+
+    double sum = 0;
+    for (map<vector<double>, int>::const_iterator entry = classCounts.begin(); entry != classCounts.end(); ++entry) {
+        int count = (entry->second*(entry->second+1))/2;
+        sum += count*integrateInteraction(expression, entry->first, entry->first, force, context);
+    }
+    for (map<vector<double>, int>::const_iterator class1 = classCounts.begin(); class1 != classCounts.end(); ++class1)
+        for (map<vector<double>, int>::const_iterator class2 = classCounts.begin(); class2 != class1; ++class2) {
+            int count = class1->second*class2->second;
+            sum += count*integrateInteraction(expression, class1->first, class2->first, force, context);
+        }
+    int numParticles = force.getNumParticles();
+    int numInteractions = (numParticles*(numParticles+1))/2;
+    sum /= numInteractions;
+    return 2*M_PI*numParticles*numParticles*sum;
+}
+
+double CustomNonbondedForceImpl::integrateInteraction(const Lepton::ExpressionProgram& expression, const vector<double>& params1, const vector<double>& params2,
+        const CustomNonbondedForce& force, const Context& context) {
+    map<std::string, double> variables;
+    for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
+        stringstream name1, name2;
+        name1 << force.getPerParticleParameterName(i) << 1;
+        name2 << force.getPerParticleParameterName(i) << 2;
+        variables[name1.str()] = params1[i];
+        variables[name2.str()] = params2[i];
+    }
+    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
+        const string& name = force.getGlobalParameterName(i);
+        variables[name] = context.getParameter(name);
+    }
+    
+    // To integrate from r_cutoff to infinity, make the change of variables x=r_cutoff/r and integrate from 0 to 1.
+    // This introduces another r^2 into the integral, which along with the r^2 in the formula for the correction
+    // 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);
+    int numPoints = 1;
+    for (int iteration = 0; iteration < 10; iteration++) {
+        double oldSum = sum;
+        double newSum = 0;
+        numPoints *= 3;
+        for (int i = 0; i < numPoints; i++) {
+            if (i%3 == 1)
+                continue;
+            double x = (i+0.5)/numPoints;
+            double r = cutoff/x;
+            variables["r"] = r;
+            double r2 = r*r;
+            newSum += expression.evaluate(variables)*r2*r2;
+        }
+        sum = newSum/numPoints + oldSum/3;
+        if (iteration > 2 && (fabs((sum-oldSum)/sum) < 1e-5 || sum == 0))
+            break;
+    }
+    return sum/cutoff;
+}

platforms/cuda/src/CudaKernels.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 "CudaBondedUtilities.h"
 #include "CudaExpressionUtilities.h"
@@ -1818,6 +1819,8 @@ CudaCalcCustomNonbondedForceKernel::~CudaCalcCustomNonbondedForceKernel() {
         delete tabulatedFunctionParams;
     for (int i = 0; i < (int) tabulatedFunctions.size(); i++)
         delete tabulatedFunctions[i];
+    if (forceCopy != NULL)
+        delete forceCopy;
 }
 
 void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const CustomNonbondedForce& force) {
@@ -1927,6 +1930,17 @@ void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const
         cu.getNonbondedUtilities().addArgument(CudaNonbondedUtilities::ParameterInfo(prefix+"globals", "float", 1, sizeof(float), globals->getDevicePointer()));
     }
     cu.addForce(new CudaCustomNonbondedForceInfo(force));
+    
+    // Record information for the long range correction.
+    
+    if (force.getNonbondedMethod() == CustomNonbondedForce::CutoffPeriodic && force.getUseLongRangeCorrection() && cu.getContextIndex() == 0) {
+        forceCopy = new CustomNonbondedForce(force);
+        hasInitializedLongRangeCorrection = false;
+    }
+    else {
+        longRangeCoefficient = 0.0;
+        hasInitializedLongRangeCorrection = true;
+    }
 }
 
 double CudaCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -1938,10 +1952,20 @@ double CudaCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool in
                 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;
+    }
+    double4 boxSize = cu.getPeriodicBoxSize();
+    return longRangeCoefficient/(boxSize.x*boxSize.y*boxSize.z);
 }
 
 void CudaCalcCustomNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force) {
@@ -1962,6 +1986,14 @@ void CudaCalcCustomNonbondedForceKernel::copyParametersToContext(ContextImpl& co
     }
     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.
     
     cu.invalidateMolecules();

platforms/cuda/src/CudaKernels.h

@@ -633,7 +633,7 @@ private:
 class CudaCalcCustomNonbondedForceKernel : public CalcCustomNonbondedForceKernel {
 public:
     CudaCalcCustomNonbondedForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system) : CalcCustomNonbondedForceKernel(name, platform),
-            cu(cu), params(NULL), globals(NULL), tabulatedFunctionParams(NULL), system(system) {
+            cu(cu), params(NULL), globals(NULL), tabulatedFunctionParams(NULL), forceCopy(NULL), system(system) {
     }
     ~CudaCalcCustomNonbondedForceKernel();
     /**
@@ -667,6 +667,9 @@ private:
     std::vector<std::string> globalParamNames;
     std::vector<float> globalParamValues;
     std::vector<CudaArray*> tabulatedFunctions;
+    double longRangeCoefficient;
+    bool hasInitializedLongRangeCorrection;
+    CustomNonbondedForce* forceCopy;
     const System& system;
 };
 

platforms/cuda/tests/TestCudaCustomNonbondedForce.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;

platforms/reference/src/ReferenceKernels.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:                                                              *
  *                                                                            *
@@ -65,6 +65,7 @@
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/internal/CustomCompoundBondForceImpl.h"
 #include "openmm/internal/CustomHbondForceImpl.h"
+#include "openmm/internal/CustomNonbondedForceImpl.h"
 #include "openmm/internal/CMAPTorsionForceImpl.h"
 #include "openmm/internal/NonbondedForceImpl.h"
 #include "openmm/internal/SplineFitter.h"
@@ -1001,6 +1002,8 @@ ReferenceCalcCustomNonbondedForceKernel::~ReferenceCalcCustomNonbondedForceKerne
     disposeIntArray(exclusionArray, numParticles);
     if (neighborList != NULL)
         delete neighborList;
+    if (forceCopy != NULL)
+        delete forceCopy;
 }
 
 void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, const CustomNonbondedForce& force) {
@@ -1059,18 +1062,32 @@ void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, c
     forceExpression = expression.differentiate("r").optimize().createProgram();
     for (int i = 0; i < numParameters; i++)
         parameterNames.push_back(force.getPerParticleParameterName(i));
-    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
         globalParameterNames.push_back(force.getGlobalParameterName(i));
+        globalParamValues[force.getGlobalParameterName(i)] = force.getGlobalParameterDefaultValue(i);
+    }
 
     // Delete the custom functions.
 
     for (map<string, Lepton::CustomFunction*>::iterator iter = functions.begin(); iter != functions.end(); iter++)
         delete iter->second;
+    
+    // Record information for the long range correction.
+    
+    if (force.getNonbondedMethod() == CustomNonbondedForce::CutoffPeriodic && force.getUseLongRangeCorrection()) {
+        forceCopy = new CustomNonbondedForce(force);
+        hasInitializedLongRangeCorrection = false;
+    }
+    else {
+        longRangeCoefficient = 0.0;
+        hasInitializedLongRangeCorrection = true;
+    }
 }
 
 double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
     vector<RealVec>& posData = extractPositions(context);
     vector<RealVec>& forceData = extractForces(context);
+    RealVec& box = extractBoxSize(context);
     RealOpenMM energy = 0;
     ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames);
     bool periodic = (nonbondedMethod == CutoffPeriodic);
@@ -1079,16 +1096,27 @@ double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bo
         ixn.setUseCutoff(nonbondedCutoff, *neighborList);
     }
     if (periodic) {
-        RealVec& box = extractBoxSize(context);
         double minAllowedSize = 2*nonbondedCutoff;
         if (box[0] < minAllowedSize || box[1] < minAllowedSize || box[2] < minAllowedSize)
             throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
         ixn.setPeriodic(box);
     }
-    map<string, double> globalParameters;
-    for (int i = 0; i < (int) globalParameterNames.size(); i++)
-        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
-    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, globalParameters, forceData, 0, includeEnergy ? &energy : NULL);
+    bool globalParamsChanged = false;
+    for (int i = 0; i < (int) globalParameterNames.size(); i++) {
+        double value = context.getParameter(globalParameterNames[i]);
+        if (globalParamValues[globalParameterNames[i]] != value)
+            globalParamsChanged = true;
+        globalParamValues[globalParameterNames[i]] = value;
+    }
+    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, globalParamValues, forceData, 0, includeEnergy ? &energy : NULL);
+    
+    // Add in the long range correction.
+    
+    if (!hasInitializedLongRangeCorrection || (globalParamsChanged && forceCopy != NULL)) {
+        longRangeCoefficient = CustomNonbondedForceImpl::calcLongRangeCorrection(*forceCopy, context.getOwner());
+        hasInitializedLongRangeCorrection = true;
+    }
+    energy += longRangeCoefficient/(box[0]*box[1]*box[2]);
     return energy;
 }
 
@@ -1106,6 +1134,14 @@ void ReferenceCalcCustomNonbondedForceKernel::copyParametersToContext(ContextImp
         for (int j = 0; j < numParameters; j++)
             particleParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
     }
+    
+    // If necessary, recompute the long range correction.
+    
+    if (forceCopy != NULL) {
+        longRangeCoefficient = CustomNonbondedForceImpl::calcLongRangeCorrection(force, context.getOwner());
+        hasInitializedLongRangeCorrection = true;
+        *forceCopy = force;
+    }
 }
 
 ReferenceCalcGBSAOBCForceKernel::~ReferenceCalcGBSAOBCForceKernel() {

platforms/reference/src/ReferenceKernels.h

@@ -586,7 +586,7 @@ private:
  */
 class ReferenceCalcCustomNonbondedForceKernel : public CalcCustomNonbondedForceKernel {
 public:
-    ReferenceCalcCustomNonbondedForceKernel(std::string name, const Platform& platform) : CalcCustomNonbondedForceKernel(name, platform) {
+    ReferenceCalcCustomNonbondedForceKernel(std::string name, const Platform& platform) : CalcCustomNonbondedForceKernel(name, platform), forceCopy(NULL) {
     }
     ~ReferenceCalcCustomNonbondedForceKernel();
     /**
@@ -616,7 +616,10 @@ private:
     int numParticles;
     int **exclusionArray;
     RealOpenMM **particleParamArray;
-    RealOpenMM nonbondedCutoff, periodicBoxSize[3];
+    RealOpenMM nonbondedCutoff, periodicBoxSize[3], longRangeCoefficient;
+    bool hasInitializedLongRangeCorrection;
+    CustomNonbondedForce* forceCopy;
+    std::map<std::string, double> globalParamValues;
     std::vector<std::set<int> > exclusions;
     Lepton::ExpressionProgram energyExpression, forceExpression;
     std::vector<std::string> parameterNames, globalParameterNames;

platforms/reference/tests/TestReferenceCustomNonbondedForce.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-2009 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -337,6 +337,82 @@ void testCoulombLennardJones() {
     }
 }
 
+void testLongRangeCorrection() {
+    // Create a box of particles.
+
+    int gridSize = 5;
+    int numParticles = gridSize*gridSize*gridSize;
+    double boxSize = gridSize*0.7;
+    double cutoff = boxSize/3;
+    ReferencePlatform platform;
+    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() {
     try {
         testSimpleExpression();
@@ -346,6 +422,7 @@ int main() {
         testPeriodic();
         testTabulatedFunction();
         testCoulombLennardJones();
+        testLongRangeCorrection();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;