Implemented long range dispersion correction

openmmapi/include/openmm/NonbondedForce.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-2009 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2010 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -245,6 +245,24 @@ public:
      *                        multiplied by this factor
      */
     void createExceptionsFromBonds(const std::vector<std::pair<int, int> >& bonds, double coulomb14Scale, double lj14Scale);
+    /**
+     * Get whether to add a contribution to the energy that approximately represents the effect of Lennard-Jones
+     * interactions beyond the cutoff distance.  The energy depends on the volume of the periodic box, and is only
+     * applicable when periodic boundary conditions are used.  When running simulations at constant pressure, adding
+     * this contribution can improve the quality of results.
+     */
+    bool getUseDispersionCorrection() const {
+        return useDispersionCorrection;
+    }
+    /**
+     * Set whether to add a contribution to the energy that approximately represents the effect of Lennard-Jones
+     * interactions beyond the cutoff distance.  The energy depends on the volume of the periodic box, and is only
+     * applicable when periodic boundary conditions are used.  When running simulations at constant pressure, adding
+     * this contribution can improve the quality of results.
+     */
+    void setUseDispersionCorrection(bool useCorrection) {
+        useDispersionCorrection = useCorrection;
+    }
 protected:
     ForceImpl* createImpl();
 private:
@@ -252,6 +270,7 @@ private:
     class ExceptionInfo;
     NonbondedMethod nonbondedMethod;
     double cutoffDistance, rfDielectric, ewaldErrorTol;
+    bool useDispersionCorrection;
     void addExclusionsToSet(const std::vector<std::set<int> >& bonded12, std::set<int>& exclusions, int baseParticle, int fromParticle, int currentLevel) const;
     std::vector<ParticleInfo> particles;
     std::vector<ExceptionInfo> exceptions;

openmmapi/include/openmm/internal/NonbondedForceImpl.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-2009 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2010 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -74,6 +74,11 @@ public:
      * Particle Mesh Ewald.
      */
     static void calcPMEParameters(const System& system, const NonbondedForce& force, double& alpha, int& xsize, int& ysize, int& zsize);
+    /**
+     * Compute the coefficient which, when divided by the periodic box volume, gives the
+     * long range dispersion correction to the energy.
+     */
+    static double calcDispersionCorrection(const System& system, const NonbondedForce& force);
 private:
     class ErrorFunction;
     class EwaldErrorFunction;

openmmapi/src/NonbondedForce.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-2009 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2010 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -46,7 +46,7 @@ using std::string;
 using std::stringstream;
 using std::vector;
 
-NonbondedForce::NonbondedForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), rfDielectric(78.3), ewaldErrorTol(5e-4) {
+NonbondedForce::NonbondedForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), rfDielectric(78.3), ewaldErrorTol(5e-4), useDispersionCorrection(false) {
 }
 
 NonbondedForce::NonbondedMethod NonbondedForce::getNonbondedMethod() const {

openmmapi/src/NonbondedForceImpl.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-2009 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2010 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -34,6 +34,7 @@
 #include "openmm/internal/NonbondedForceImpl.h"
 #include "openmm/kernels.h"
 #include <cmath>
+#include <map>
 #include <sstream>
 
 #ifndef M_PI
@@ -94,6 +95,10 @@ void NonbondedForceImpl::initialize(ContextImpl& context) {
         if (cutoff > 0.5*boxVectors[0][0] || cutoff > 0.5*boxVectors[1][1] || cutoff > 0.5*boxVectors[2][2])
             throw OpenMMException("NonbondedForce: The cutoff distance cannot be greater than half the periodic box size.");
     }
+    else {
+        if (owner.getUseDispersionCorrection())
+            throw OpenMMException("NonbondedForce: Dispersion correction may only be used with periodic boundary conditions.");
+    }
     dynamic_cast<CalcNonbondedForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
 }
 
@@ -169,3 +174,49 @@ int NonbondedForceImpl::findZero(const NonbondedForceImpl::ErrorFunction& f, int
         value = f.getValue(++arg);
     return arg;
 }
+
+double NonbondedForceImpl::calcDispersionCorrection(const System& system, const NonbondedForce& force) {
+    // Identify all particle classes (defined by sigma and epsilon), and count the number of
+    // particles in each class.
+
+    map<pair<double, double>, int> classCounts;
+    for (int i = 0; i < force.getNumParticles(); i++) {
+        double charge, sigma, epsilon;
+        force.getParticleParameters(i, charge, sigma, epsilon);
+        pair<double, double> key = make_pair(sigma, epsilon);
+        map<pair<double, double>, int>::iterator entry = classCounts.find(key);
+        if (entry == classCounts.end())
+            classCounts[key] = 1;
+        else
+            entry->second++;
+    }
+
+    // Loop over all pairs of classes to compute the coefficient.
+
+    double sum1 = 0, sum2 = 0;
+    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;
+        int count = (entry->second*(entry->second+1))/2;
+        double sigma2 = sigma*sigma;
+        double sigma6 = sigma2*sigma2*sigma2;
+        sum1 += count*epsilon*sigma6*sigma6;
+        sum2 += count*epsilon*sigma6;
+    }
+    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) {
+            double sigma = 0.5*(class1->first.first+class2->first.first);
+            double epsilon = sqrt(class1->first.second*class2->first.second);
+            int count = class1->second*class2->second;
+            double sigma2 = sigma*sigma;
+            double sigma6 = sigma2*sigma2*sigma2;
+            sum1 += count*epsilon*sigma6*sigma6;
+            sum2 += count*epsilon*sigma6;
+        }
+    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)));
+}

platforms/cuda/include/CudaPlatform.h

@@ -79,7 +79,7 @@ public:
     int nonbondedMethod, customNonbondedMethod;
     int cmMotionFrequency;
     int stepCount, computeForceCount;
-    double time, ewaldSelfEnergy;
+    double time, ewaldSelfEnergy, dispersionCoefficient;
     std::map<std::string, std::string> propertyValues;
 };
 

platforms/cuda/src/CudaKernels.cpp

@@ -102,7 +102,10 @@ double CudaCalcForcesAndEnergyKernel::finishEnergyComputation(ContextImpl& conte
     if (data.hasCustomNonbonded)
         kCalculateCustomNonbondedForces(gpu, data.hasNonbonded);
     kCalculateLocalForces(gpu);
-    return kReduceEnergy(gpu)+data.ewaldSelfEnergy;
+    double energy = kReduceEnergy(gpu)+data.ewaldSelfEnergy;
+    if (data.dispersionCoefficient != 0.0)
+        energy += data.dispersionCoefficient/(gpu->sim.periodicBoxSizeX*gpu->sim.periodicBoxSizeY*gpu->sim.periodicBoxSizeZ);
+    return energy;
 }
 
 void CudaUpdateStateDataKernel::initialize(const System& system) {
@@ -538,6 +541,13 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
                 for (int i = 0; i < numParticles; i++)
                     data.ewaldSelfEnergy -= selfEnergyScale*q[i]*q[i];
         }
+
+        // Compute the long range dispersion correction.
+
+        if (force.getUseDispersionCorrection())
+            data.dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(system, force);
+        else
+            data.dispersionCoefficient = 0.0;
     }
 
     // Initialize 1-4 nonbonded interactions.

platforms/cuda/tests/TestCudaMonteCarloBarostat.cpp

@@ -204,6 +204,7 @@ void testWater() {
     system.setDefaultPeriodicBoxVectors(Vec3(gridSize*spacing, 0, 0), Vec3(0, gridSize*spacing, 0), Vec3(0, 0, gridSize*spacing));
     NonbondedForce* nonbonded = new NonbondedForce();
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    nonbonded->setUseDispersionCorrection(true);
     vector<Vec3> positions;
     Vec3 offset1(dOH, 0, 0);
     Vec3 offset2(dOH*std::cos(angle), dOH*std::sin(angle), 0);
@@ -249,7 +250,7 @@ void testWater() {
     }
     volume /= steps;
     double density = numMolecules*18/(AVOGADRO*volume*1e-21);
-    ASSERT_USUALLY_EQUAL_TOL(1.0, density, 0.1);
+    ASSERT_USUALLY_EQUAL_TOL(1.0, density, 0.02);
 }
 
 int main() {

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 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2010 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -600,6 +600,76 @@ void testBlockInteractions(bool periodic) {
         }
 }
 
+void testDispersionCorrection() {
+    // Create a box full of identical particles.
+
+    int gridSize = 5;
+    int numParticles = gridSize*gridSize*gridSize;
+    double boxSize = gridSize*0.5;
+    double cutoff = boxSize/3;
+    CudaPlatform platform;
+    System system;
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    vector<Vec3> positions(numParticles);
+    int index = 0;
+    for (int i = 0; i < gridSize; i++)
+        for (int j = 0; j < gridSize; j++)
+            for (int k = 0; k < gridSize; k++) {
+                system.addParticle(1.0);
+                nonbonded->addParticle(0, 1.1, 0.5);
+                positions[index] = Vec3(i*boxSize/gridSize, j*boxSize/gridSize, k*boxSize/gridSize);
+                index++;
+            }
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    nonbonded->setCutoffDistance(cutoff);
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.addForce(nonbonded);
+
+    // See if the correction has the correct value.
+
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    double energy1 = context.getState(State::Energy).getPotentialEnergy();
+    nonbonded->setUseDispersionCorrection(true);
+    context.reinitialize();
+    context.setPositions(positions);
+    double energy2 = context.getState(State::Energy).getPotentialEnergy();
+    double term1 = (0.5*pow(1.1, 12)/pow(cutoff, 9))/9;
+    double term2 = (0.5*pow(1.1, 6)/pow(cutoff, 3))/3;
+    double expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
+    ASSERT_EQUAL_TOL(expected, energy2-energy1, 1e-4);
+
+    // Now modify half the particles to be different, and see if it is still correct.
+
+    int numType2 = 0;
+    for (int i = 0; i < numParticles; i += 2) {
+        nonbonded->setParticleParameters(i, 0, 1, 1);
+        numType2++;
+    }
+    int numType1 = numParticles-numType2;
+    nonbonded->setUseDispersionCorrection(false);
+    context.reinitialize();
+    context.setPositions(positions);
+    energy1 = context.getState(State::Energy).getPotentialEnergy();
+    nonbonded->setUseDispersionCorrection(true);
+    context.reinitialize();
+    context.setPositions(positions);
+    energy2 = context.getState(State::Energy).getPotentialEnergy();
+    term1 = ((numType1*(numType1+1))/2)*(0.5*pow(1.1, 12)/pow(cutoff, 9))/9;
+    term2 = ((numType1*(numType1+1))/2)*(0.5*pow(1.1, 6)/pow(cutoff, 3))/3;
+    term1 += ((numType2*(numType2+1))/2)*(1*pow(1.0, 12)/pow(cutoff, 9))/9;
+    term2 += ((numType2*(numType2+1))/2)*(1*pow(1.0, 6)/pow(cutoff, 3))/3;
+    double combinedSigma = 0.5*(1+1.1);
+    double combinedEpsilon = sqrt(1*0.5);
+    term1 += (numType1*numType2)*(combinedEpsilon*pow(combinedSigma, 12)/pow(cutoff, 9))/9;
+    term2 += (numType1*numType2)*(combinedEpsilon*pow(combinedSigma, 6)/pow(cutoff, 3))/3;
+    term1 /= (numParticles*(numParticles+1))/2;
+    term2 /= (numParticles*(numParticles+1))/2;
+    expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
+    ASSERT_EQUAL_TOL(expected, energy2-energy1, 1e-4);
+}
+
 int main() {
     try {
         testCoulomb();
@@ -611,6 +681,7 @@ int main() {
         testLargeSystem();
         testBlockInteractions(false);
         testBlockInteractions(true);
+        testDispersionCorrection();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/opencl/src/OpenCLKernels.cpp

@@ -1098,6 +1098,10 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
         defines["REACTION_FIELD_K"] = doubleToString(reactionFieldK);
         defines["REACTION_FIELD_C"] = doubleToString(reactionFieldC);
     }
+    if (force.getUseDispersionCorrection())
+        dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(system, force);
+    else
+        dispersionCoefficient = 0.0;
     double alpha = 0;
     if (force.getNonbondedMethod() == NonbondedForce::Ewald) {
         // Compute the Ewald parameters.
@@ -1358,7 +1362,12 @@ void OpenCLCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
 
 double OpenCLCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
     executeForces(context);
-    return ewaldSelfEnergy;
+    double energy = ewaldSelfEnergy;
+    if (dispersionCoefficient != 0.0) {
+        mm_float4 boxSize = cl.getPeriodicBoxSize();
+        energy += dispersionCoefficient/(boxSize.x*boxSize.y*boxSize.z);
+    }
+    return energy;
 }
 
 class OpenCLCustomNonbondedForceInfo : public OpenCLForceInfo {

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-2009 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2010 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -505,7 +505,7 @@ private:
     cl::Kernel pmeConvolutionKernel;
     cl::Kernel pmeInterpolateForceKernel;
     std::map<std::string, std::string> pmeDefines;
-    double ewaldSelfEnergy;
+    double ewaldSelfEnergy, dispersionCoefficient;
     static const int PmeOrder = 5;
 };
 

platforms/opencl/tests/TestOpenCLMonteCarloBarostat.cpp

@@ -204,6 +204,7 @@ void testWater() {
     system.setDefaultPeriodicBoxVectors(Vec3(gridSize*spacing, 0, 0), Vec3(0, gridSize*spacing, 0), Vec3(0, 0, gridSize*spacing));
     NonbondedForce* nonbonded = new NonbondedForce();
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    nonbonded->setUseDispersionCorrection(true);
     vector<Vec3> positions;
     Vec3 offset1(dOH, 0, 0);
     Vec3 offset2(dOH*std::cos(angle), dOH*std::sin(angle), 0);
@@ -249,7 +250,7 @@ void testWater() {
     }
     volume /= steps;
     double density = numMolecules*18/(AVOGADRO*volume*1e-21);
-    ASSERT_USUALLY_EQUAL_TOL(1.0, density, 0.1);
+    ASSERT_USUALLY_EQUAL_TOL(1.0, density, 0.02);
 }
 
 int main() {

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-2009 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2010 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -605,6 +605,76 @@ void testBlockInteractions(bool periodic) {
         }
 }
 
+void testDispersionCorrection() {
+    // Create a box full of identical particles.
+
+    int gridSize = 5;
+    int numParticles = gridSize*gridSize*gridSize;
+    double boxSize = gridSize*0.5;
+    double cutoff = boxSize/3;
+    OpenCLPlatform platform;
+    System system;
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    vector<Vec3> positions(numParticles);
+    int index = 0;
+    for (int i = 0; i < gridSize; i++)
+        for (int j = 0; j < gridSize; j++)
+            for (int k = 0; k < gridSize; k++) {
+                system.addParticle(1.0);
+                nonbonded->addParticle(0, 1.1, 0.5);
+                positions[index] = Vec3(i*boxSize/gridSize, j*boxSize/gridSize, k*boxSize/gridSize);
+                index++;
+            }
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    nonbonded->setCutoffDistance(cutoff);
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.addForce(nonbonded);
+
+    // See if the correction has the correct value.
+
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    double energy1 = context.getState(State::Energy).getPotentialEnergy();
+    nonbonded->setUseDispersionCorrection(true);
+    context.reinitialize();
+    context.setPositions(positions);
+    double energy2 = context.getState(State::Energy).getPotentialEnergy();
+    double term1 = (0.5*pow(1.1, 12)/pow(cutoff, 9))/9;
+    double term2 = (0.5*pow(1.1, 6)/pow(cutoff, 3))/3;
+    double expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
+    ASSERT_EQUAL_TOL(expected, energy2-energy1, 1e-4);
+
+    // Now modify half the particles to be different, and see if it is still correct.
+
+    int numType2 = 0;
+    for (int i = 0; i < numParticles; i += 2) {
+        nonbonded->setParticleParameters(i, 0, 1, 1);
+        numType2++;
+    }
+    int numType1 = numParticles-numType2;
+    nonbonded->setUseDispersionCorrection(false);
+    context.reinitialize();
+    context.setPositions(positions);
+    energy1 = context.getState(State::Energy).getPotentialEnergy();
+    nonbonded->setUseDispersionCorrection(true);
+    context.reinitialize();
+    context.setPositions(positions);
+    energy2 = context.getState(State::Energy).getPotentialEnergy();
+    term1 = ((numType1*(numType1+1))/2)*(0.5*pow(1.1, 12)/pow(cutoff, 9))/9;
+    term2 = ((numType1*(numType1+1))/2)*(0.5*pow(1.1, 6)/pow(cutoff, 3))/3;
+    term1 += ((numType2*(numType2+1))/2)*(1*pow(1.0, 12)/pow(cutoff, 9))/9;
+    term2 += ((numType2*(numType2+1))/2)*(1*pow(1.0, 6)/pow(cutoff, 3))/3;
+    double combinedSigma = 0.5*(1+1.1);
+    double combinedEpsilon = sqrt(1*0.5);
+    term1 += (numType1*numType2)*(combinedEpsilon*pow(combinedSigma, 12)/pow(cutoff, 9))/9;
+    term2 += (numType1*numType2)*(combinedEpsilon*pow(combinedSigma, 6)/pow(cutoff, 3))/3;
+    term1 /= (numParticles*(numParticles+1))/2;
+    term2 /= (numParticles*(numParticles+1))/2;
+    expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
+    ASSERT_EQUAL_TOL(expected, energy2-energy1, 1e-4);
+}
+
 int main() {
     try {
         testCoulomb();
@@ -616,6 +686,7 @@ int main() {
         testLargeSystem();
         testBlockInteractions(false);
         testBlockInteractions(true);
+        testDispersionCorrection();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/reference/src/ReferenceKernels.cpp

@@ -697,6 +697,10 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
         ewaldAlpha = (RealOpenMM) alpha;
     }
     rfDielectric = (RealOpenMM)force.getReactionFieldDielectric();
+    if (force.getUseDispersionCorrection())
+        dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(system, force);
+    else
+        dispersionCoefficient = 0.0;
 }
 
 void ReferenceCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
@@ -749,6 +753,10 @@ double ReferenceCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
     refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, energyArray, 0, &energy, nonbonded14);
     disposeRealArray(forceData, numParticles);
     delete[] energyArray;
+    if (periodic || ewald || pme) {
+        RealOpenMM* boxSize = extractBoxSize(context);
+        energy += dispersionCoefficient/(boxSize[0]*boxSize[1]*boxSize[2]);
+    }
     return energy;
 }
 

platforms/reference/src/ReferenceKernels.h

@@ -452,7 +452,7 @@ private:
     int numParticles, num14;
     int **exclusionArray, **bonded14IndexArray;
     RealOpenMM **particleParamArray, **bonded14ParamArray;
-    RealOpenMM nonbondedCutoff, rfDielectric, ewaldAlpha;
+    RealOpenMM nonbondedCutoff, rfDielectric, ewaldAlpha, dispersionCoefficient;
     int kmax[3], gridSize[3];
     std::vector<std::set<int> > exclusions;
     NonbondedMethod nonbondedMethod;

platforms/reference/tests/TestReferenceNonbondedForce.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 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2010 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -353,6 +353,76 @@ void testPeriodic() {
     ASSERT_EQUAL_TOL(2*ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
 }
 
+void testDispersionCorrection() {
+    // Create a box full of identical particles.
+
+    int gridSize = 5;
+    int numParticles = gridSize*gridSize*gridSize;
+    double boxSize = gridSize*0.5;
+    double cutoff = boxSize/3;
+    ReferencePlatform platform;
+    System system;
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    vector<Vec3> positions(numParticles);
+    int index = 0;
+    for (int i = 0; i < gridSize; i++)
+        for (int j = 0; j < gridSize; j++)
+            for (int k = 0; k < gridSize; k++) {
+                system.addParticle(1.0);
+                nonbonded->addParticle(0, 1.1, 0.5);
+                positions[index] = Vec3(i*boxSize/gridSize, j*boxSize/gridSize, k*boxSize/gridSize);
+                index++;
+            }
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    nonbonded->setCutoffDistance(cutoff);
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.addForce(nonbonded);
+
+    // See if the correction has the correct value.
+
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    double energy1 = context.getState(State::Energy).getPotentialEnergy();
+    nonbonded->setUseDispersionCorrection(true);
+    context.reinitialize();
+    context.setPositions(positions);
+    double energy2 = context.getState(State::Energy).getPotentialEnergy();
+    double term1 = (0.5*pow(1.1, 12)/pow(cutoff, 9))/9;
+    double term2 = (0.5*pow(1.1, 6)/pow(cutoff, 3))/3;
+    double expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
+    ASSERT_EQUAL_TOL(expected, energy2-energy1, 1e-4);
+
+    // Now modify half the particles to be different, and see if it is still correct.
+
+    int numType2 = 0;
+    for (int i = 0; i < numParticles; i += 2) {
+        nonbonded->setParticleParameters(i, 0, 1, 1);
+        numType2++;
+    }
+    int numType1 = numParticles-numType2;
+    nonbonded->setUseDispersionCorrection(false);
+    context.reinitialize();
+    context.setPositions(positions);
+    energy1 = context.getState(State::Energy).getPotentialEnergy();
+    nonbonded->setUseDispersionCorrection(true);
+    context.reinitialize();
+    context.setPositions(positions);
+    energy2 = context.getState(State::Energy).getPotentialEnergy();
+    term1 = ((numType1*(numType1+1))/2)*(0.5*pow(1.1, 12)/pow(cutoff, 9))/9;
+    term2 = ((numType1*(numType1+1))/2)*(0.5*pow(1.1, 6)/pow(cutoff, 3))/3;
+    term1 += ((numType2*(numType2+1))/2)*(1*pow(1.0, 12)/pow(cutoff, 9))/9;
+    term2 += ((numType2*(numType2+1))/2)*(1*pow(1.0, 6)/pow(cutoff, 3))/3;
+    double combinedSigma = 0.5*(1+1.1);
+    double combinedEpsilon = sqrt(1*0.5);
+    term1 += (numType1*numType2)*(combinedEpsilon*pow(combinedSigma, 12)/pow(cutoff, 9))/9;
+    term2 += (numType1*numType2)*(combinedEpsilon*pow(combinedSigma, 6)/pow(cutoff, 3))/3;
+    term1 /= (numParticles*(numParticles+1))/2;
+    term2 /= (numParticles*(numParticles+1))/2;
+    expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
+    ASSERT_EQUAL_TOL(expected, energy2-energy1, 1e-4);
+}
+
 int main() {
     try {
         testCoulomb();
@@ -361,6 +431,7 @@ int main() {
         testCutoff();
         testCutoff14();
         testPeriodic();
+        testDispersionCorrection();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;