Lee-Ping's vdw dispersion correction

added tests to TestCudaAmoebaVdwForce for dispersion correction using box of 216 water molecules
AmoebaVdwForce and AmoebaMultipoleForce now throw execptions if cutoff > 0.5*(box size)

plugins/amoeba/openmmapi/include/openmm/AmoebaVdwForce.h

@@ -129,6 +129,26 @@ public:
      */
     const std::string& getEpsilonCombiningRule( void ) const;
 
+    /**
+     * LPW: Get whether to add a contribution to the energy that approximately represents the effect of VdW
+     * 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 VdW
+     * 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;
+    }
+
     /**
      * Set exclusions for specified particle
      * 
@@ -195,6 +215,7 @@ private:
     int usePBC;
     int useNeighborList;
     double cutoff;
+    bool useDispersionCorrection;
 
     std::string sigmaCombiningRule;
     std::string epsilonCombiningRule;
@@ -234,3 +255,4 @@ public:
 } // namespace OpenMM
 
 #endif /*OPENMM_AMOEBA_VDW_FORCE_H_*/
+

plugins/amoeba/openmmapi/include/openmm/internal/AmoebaVdwForceImpl.h

@@ -41,6 +41,8 @@
 
 namespace OpenMM {
 
+class System;
+
 /**
  * This is the internal implementation of AmoebaVdwForce.
  */
@@ -61,6 +63,11 @@ public:
         return std::map<std::string, double>(); // This force field doesn't define any parameters.
     }
     std::vector<std::string> getKernelNames();
+    /**
+     * 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 AmoebaVdwForce& force);
 private:
     AmoebaVdwForce& owner;
     Kernel kernel;
@@ -69,3 +76,4 @@ private:
 } // namespace OpenMM
 
 #endif /*OPENMM_AMOEBA_VDW_FORCE_IMPL_H_*/
+

plugins/amoeba/openmmapi/src/AmoebaMultipoleForceImpl.cpp

@@ -53,6 +53,16 @@ void AmoebaMultipoleForceImpl::initialize(ContextImpl& context) {
     if (owner.getNumMultipoles() != system.getNumParticles())
         throw OpenMMException("AmoebaMultipoleForce must have exactly as many particles as the System it belongs to.");
 
+    // check cutoff < 0.5*boxSize
+
+    if (owner.getNonbondedMethod() == AmoebaMultipoleForce::PME) {
+        Vec3 boxVectors[3];
+        system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+        double cutoff = owner.getCutoffDistance();
+        if (cutoff > 0.5*boxVectors[0][0] || cutoff > 0.5*boxVectors[1][1] || cutoff > 0.5*boxVectors[2][2])
+            throw OpenMMException("AmoebaMultipoleForce: The cutoff distance cannot be greater than half the periodic box size.");
+    }   
+
     double quadrupoleValidationTolerance = 1.0e-05;
     for( int ii = 0; ii < system.getNumParticles(); ii++ ){
 

plugins/amoeba/openmmapi/src/AmoebaVdwForce.cpp

@@ -38,7 +38,7 @@ using namespace OpenMM;
 using std::string;
 using std::vector;
 
-AmoebaVdwForce::AmoebaVdwForce() : sigmaCombiningRule("CUBIC-MEAN"), epsilonCombiningRule("HHG"), usePBC(0), cutoff(1.0e+10), useNeighborList(0) {
+AmoebaVdwForce::AmoebaVdwForce() : sigmaCombiningRule("CUBIC-MEAN"), epsilonCombiningRule("HHG"), usePBC(0), cutoff(1.0e+10), useNeighborList(0), useDispersionCorrection(true) {
 }
 
 int AmoebaVdwForce::addParticle(int ivIndex, int classIndex, double sigma, double epsilon, double reductionFactor ) {
@@ -131,3 +131,4 @@ int AmoebaVdwForce::getPBC( void ) const {
 ForceImpl* AmoebaVdwForce::createImpl() {
     return new AmoebaVdwForceImpl(*this);
 }
+

plugins/amoeba/openmmapi/src/AmoebaVdwForceImpl.cpp

@@ -32,8 +32,11 @@
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/internal/AmoebaVdwForceImpl.h"
 #include "openmm/amoebaKernels.h"
+#include <map>
+#include <stdio.h>
 
 using namespace OpenMM;
+using namespace std;
 
 using std::pair;
 using std::vector;
@@ -51,6 +54,16 @@ void AmoebaVdwForceImpl::initialize(ContextImpl& context) {
     if (owner.getNumParticles() != system.getNumParticles())
         throw OpenMMException("AmoebaVdwForce must have exactly as many particles as the System it belongs to.");
 
+    // check that cutoff < 0.5*boxSize
+
+    if (owner.getPBC()) {
+        Vec3 boxVectors[3];
+        system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+        double cutoff = owner.getCutoff();
+        if (cutoff > 0.5*boxVectors[0][0] || cutoff > 0.5*boxVectors[1][1] || cutoff > 0.5*boxVectors[2][2])
+            throw OpenMMException("AmoebaVdwForce: The cutoff distance cannot be greater than half the periodic box size.");
+    }   
+
     kernel = context.getPlatform().createKernel(CalcAmoebaVdwForceKernel::Name(), context);
     kernel.getAs<CalcAmoebaVdwForceKernel>().initialize(context.getSystem(), owner);
 }
@@ -61,9 +74,171 @@ double AmoebaVdwForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includ
     return 0.0;
 }
 
+double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const AmoebaVdwForce& force) {
+
+    // Amoeba VdW dispersion correction implemented by LPW
+    // There is no dispersion correction if PBC is off or the cutoff is set to the default value of ten billion (AmoebaVdwForce.cpp)
+    if (force.getPBC() == 0 || force.getUseNeighborList() == 0)
+        return 0.0;
+
+    // Identify all particle classes (defined by sigma and epsilon and reduction), and count the number of
+    // particles in each class.
+
+    map<pair<double, double>, int> classCounts;
+    for (int i = 0; i < force.getNumParticles(); i++) {
+        double sigma, epsilon, reduction;
+        // The variables reduction, ivindex, classindex are not used.
+        int ivindex, classindex;
+        // Get the sigma and epsilon parameters, ignoring everything else.
+        force.getParticleParameters(i, ivindex, classindex, sigma, epsilon, reduction);
+        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++;
+    }
+
+    // Compute the VdW tapering coefficients.  Mostly copied from amoebaCudaGpu.cpp.
+    double cutoff = force.getCutoff();
+    double vdwTaper = 0.90; // vdwTaper is a scaling factor, it is not a distance.
+    double c0 = 0.0;
+    double c1 = 0.0;
+    double c2 = 0.0;
+    double c3 = 0.0;
+    double c4 = 0.0;
+    double c5 = 0.0;
+
+    double vdwCut = cutoff;
+    double vdwTaperCut = vdwTaper*cutoff;
+
+    double vdwCut2 = vdwCut*vdwCut;
+    double vdwCut3 = vdwCut2*vdwCut;
+    double vdwCut4 = vdwCut2*vdwCut2;
+    double vdwCut5 = vdwCut2*vdwCut3;
+    double vdwCut6 = vdwCut3*vdwCut3;
+    double vdwCut7 = vdwCut3*vdwCut4;
+
+    double vdwTaperCut2 = vdwTaperCut*vdwTaperCut;
+    double vdwTaperCut3 = vdwTaperCut2*vdwTaperCut;
+    double vdwTaperCut4 = vdwTaperCut2*vdwTaperCut2;
+    double vdwTaperCut5 = vdwTaperCut2*vdwTaperCut3;
+    double vdwTaperCut6 = vdwTaperCut3*vdwTaperCut3;
+    double vdwTaperCut7 = vdwTaperCut3*vdwTaperCut4;
+
+    // get 5th degree multiplicative switching function coefficients;
+
+    double denom = 1.0 / (vdwCut - vdwTaperCut);
+    double denom2 = denom*denom;
+    denom = denom * denom2*denom2;
+
+    c0 = vdwCut * vdwCut2 * (vdwCut2 - 5.0 * vdwCut * vdwTaperCut + 10.0 * vdwTaperCut2) * denom;
+    c1 = -30.0 * vdwCut2 * vdwTaperCut2*denom;
+    c2 = 30.0 * (vdwCut2 * vdwTaperCut + vdwCut * vdwTaperCut2) * denom;
+    c3 = -10.0 * (vdwCut2 + 4.0 * vdwCut * vdwTaperCut + vdwTaperCut2) * denom;
+    c4 = 15.0 * (vdwCut + vdwTaperCut) * denom;
+    c5 = -6.0 * denom;
+
+    // Loop over all pairs of classes to compute the coefficient.
+    // Copied over from TINKER - numerical integration.
+    double range = 20.0;
+    double cut = vdwTaperCut; // This is where tapering BEGINS
+    double off = vdwCut; // This is where tapering ENDS
+    int nstep = 200;
+    int ndelta = int(double(nstep) * (range - cut));
+    double rdelta = (range - cut) / double(ndelta);
+    double offset = cut - 0.5 * rdelta;
+    double dhal = 0.07; // This magic number also appears in kCalculateAmoebaCudaVdw14_7.cu
+    double ghal = 0.12; // This magic number also appears in kCalculateAmoebaCudaVdw14_7.cu
+    double elrc = 0.0; // This number is incremented and passed out at the end
+    double e = 0.0;
+    double sigma, epsilon; // The pairwise sigma and epsilon parameters.
+    int i = 0, k = 0; // Loop counters.
+
+    // Double loop over different atom types.
+    std::string sigmaCombiningRule = force.getSigmaCombiningRule();
+    std::string epsilonCombiningRule = force.getEpsilonCombiningRule();
+    for (map<pair<double, double>, int>::const_iterator class1 = classCounts.begin(); class1 != classCounts.end(); ++class1) {
+        k = 0;
+        for (map<pair<double, double>, int>::const_iterator class2 = classCounts.begin(); class2 != classCounts.end(); ++class2) { 
+            // AMOEBA combining rules, copied over from the CUDA code.
+            double iSigma = class1->first.first;
+            double jSigma = class2->first.first;
+            double iEpsilon = class1->first.second;
+            double jEpsilon = class2->first.second;
+            // ARITHMETIC = 1
+            // GEOMETRIC  = 2
+            // CUBIC-MEAN = 3
+            if (sigmaCombiningRule == "ARITHMETIC") {
+                sigma = iSigma + jSigma;
+            } else if (sigmaCombiningRule == "GEOMETRIC") {
+                sigma = 2.0f * std::sqrt(iSigma * jSigma);
+            } else {
+                double iSigma2 = iSigma*iSigma;
+                double jSigma2 = jSigma*jSigma;
+                sigma = 2.0f * (iSigma2 * iSigma + jSigma2 * jSigma) / (iSigma2 + jSigma2);
+            }
+            // ARITHMETIC = 1
+            // GEOMETRIC  = 2
+            // HARMONIC   = 3
+            // HHG        = 4
+            if (epsilonCombiningRule == "ARITHMETIC") {
+                epsilon = 0.5f * (iEpsilon + jEpsilon);
+            } else if (epsilonCombiningRule == "GEOMETRIC") {
+                epsilon = std::sqrt(iEpsilon * jEpsilon);
+            } else if (epsilonCombiningRule == "HARMONIC") {
+                epsilon = 2.0f * (iEpsilon * jEpsilon) / (iEpsilon + jEpsilon);
+            } else {
+                double epsilonS = std::sqrt(iEpsilon) + std::sqrt(jEpsilon);
+                epsilon = 4.0f * (iEpsilon * jEpsilon) / (epsilonS * epsilonS);
+            }
+            int count = class1->second * class2->second;
+            // Below is an exact copy of stuff from the previous block.
+            double rv = sigma;
+            double termik = 2.0 * M_PI * count; // termik is equivalent to 2 * pi * count.
+            double rv2 = rv * rv;
+            double rv6 = rv2 * rv2 * rv2;
+            double rv7 = rv6 * rv;
+            double etot = 0.0;
+            double r2 = 0.0;
+            for (int j = 1; j <= ndelta; j++) {
+                double r = offset + double(j) * rdelta;
+                r2 = r*r;
+                double r3 = r2 * r;
+                double r6 = r3 * r3;
+                double r7 = r6 * r;
+                // The following is for buffered 14-7 only.
+                /*
+                double rho = r/rv;
+                double term1 = pow(((dhal + 1.0) / (dhal + rho)),7);
+                double term2 = ((ghal + 1.0) / (ghal + pow(rho,7))) - 2.0;
+                e = epsilon * term1 * term2;
+                */
+                double rho = r7 + ghal*rv7;
+                double tau = (dhal + 1.0) / (r + dhal * rv);
+                double tau7 = pow(tau, 7);
+                e = epsilon * rv7 * tau7 * ((ghal + 1.0) * rv7 / rho - 2.0);
+                double taper = 0.0;
+                if (r < off) {
+                    double r4 = r2 * r2;
+                    double r5 = r2 * r3;
+                    taper = c5 * r5 + c4 * r4 + c3 * r3 + c2 * r2 + c1 * r + c0;
+                    e = e * (1.0 - taper);
+                }
+                etot = etot + e * rdelta * r2;
+            }
+            elrc = elrc + termik * etot;
+            k++;
+        }
+        i++;
+    }
+    return elrc;
+}
+
 std::vector<std::string> AmoebaVdwForceImpl::getKernelNames() {
     std::vector<std::string> names;
     names.push_back(CalcAmoebaVdwForceKernel::Name());
     return names;
 }
 
+

plugins/amoeba/platforms/cuda/src/AmoebaCudaData.h

@@ -197,6 +197,8 @@ public:
 
     CudaPlatform::PlatformData& cudaPlatformData;
 
+    double dispersionCoefficient;
+    
 private:
 
     amoebaGpuContext amoebaGpu;
@@ -217,3 +219,4 @@ private:
 } // namespace OpenMM
 
 #endif /*AMOEBA_CUDA_DATA_H_*/
+

plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp

@@ -29,12 +29,14 @@
 #include "kernels/amoebaGpuTypes.h"
 #include "kernels/cudaKernels.h"
 #include "kernels/amoebaCudaKernels.h"
+#include "openmm/internal/AmoebaVdwForceImpl.h"
 #include "openmm/internal/AmoebaMultipoleForceImpl.h"
 #include "openmm/internal/AmoebaWcaDispersionForceImpl.h"
 #include "openmm/internal/AmoebaTorsionTorsionForceImpl.h"
 #include "openmm/internal/NonbondedForceImpl.h"
 #include "CudaForceInfo.h"
 
+#include <stdio.h>
 #include <cmath>
 #ifdef _MSC_VER
 #include <windows.h>
@@ -1239,11 +1241,21 @@ void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoeba
                         force.getPBC(), force.getUseNeighborList() );
     }
     data.setUseVdwNeighborList( force.getUseNeighborList() );
+    if (force.getUseDispersionCorrection())
+        data.dispersionCoefficient = AmoebaVdwForceImpl::calcDispersionCorrection(system, force);
+    else
+        data.dispersionCoefficient = 0.0;               
 }
 
 double CudaCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    _gpuContext* gpu = data.cudaPlatformData.gpu;
     computeAmoebaVdwForce( data );
+    if (data.dispersionCoefficient != 0.0) {
+        double Answer = data.dispersionCoefficient/(gpu->sim.periodicBoxSizeX*gpu->sim.periodicBoxSizeY*gpu->sim.periodicBoxSizeZ);
+        return Answer;
+    } else {
         return 0.0;
+    }
 }
 
 /* -------------------------------------------------------------------------- *
@@ -1318,3 +1330,4 @@ double CudaCalcAmoebaWcaDispersionForceKernel::execute(ContextImpl& context, boo
     computeAmoebaWcaDispersionForce( data );
     return 0.0;
 }
+

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp

@@ -618,6 +618,10 @@ void ReferenceCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextI
 
 ReferenceCalcAmoebaVdwForceKernel::ReferenceCalcAmoebaVdwForceKernel(std::string name, const Platform& platform, System& system) :
        CalcAmoebaVdwForceKernel(name, platform), system(system) {
+    useNeighborList = 0;
+    usePBC = 0;
+    cutoff = 1.0e+10;
+
 }
 
 ReferenceCalcAmoebaVdwForceKernel::~ReferenceCalcAmoebaVdwForceKernel() {
@@ -656,13 +660,26 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A
     }   
     sigmaCombiningRule   = force.getSigmaCombiningRule();
     epsilonCombiningRule = force.getEpsilonCombiningRule();
+    useNeighborList      = force.getUseNeighborList();
+    usePBC               = force.getPBC();
+    cutoff               = force.getCutoff();
 }
 
 double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
 
     vector<RealVec>& posData   = extractPositions(context);
     vector<RealVec>& forceData = extractForces(context);
-    AmoebaReferenceVdwForce vdwForce( sigmaCombiningRule, epsilonCombiningRule, AmoebaReferenceVdwForce::NoCutoff );
+    AmoebaReferenceVdwForce vdwForce( sigmaCombiningRule, epsilonCombiningRule );
+    if( useNeighborList ){
+        vdwForce.setCutoff( cutoff );
+        if( usePBC ){
+            vdwForce.setNonbondedMethod( AmoebaReferenceVdwForce::CutoffPeriodic);
+        } else {
+            vdwForce.setNonbondedMethod( AmoebaReferenceVdwForce::CutoffNonPeriodic);
+        }
+    } else {
+        vdwForce.setNonbondedMethod( AmoebaReferenceVdwForce::NoCutoff );
+    }
     RealOpenMM energy = vdwForce.calculateForceAndEnergy( numParticles, posData, indexIVs, indexClasses, sigmas, epsilons, reductions, allExclusions, forceData);
     return static_cast<double>(energy);
 }

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h

@@ -29,6 +29,7 @@
 
 #include "openmm/System.h"
 #include "openmm/amoebaKernels.h"
+#include "SimTKReference/ReferenceNeighborList.h"
 #include "SimTKUtilities/SimTKOpenMMRealType.h"
 
 namespace OpenMM {
@@ -430,33 +431,6 @@ private:
     System& system;
 };
 
-// /**
-//  * This kernel is invoked by AmoebaMultipoleForce to calculate the forces acting on the system and the energy of the system.
-//  */
-// class ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel : public CalcAmoebaGeneralizedKirkwoodForceKernel {
-// public:
-//     ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel(std::string name, const Platform& platform, System& system);
-//     ~ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel();
-//     /**
-//      * Initialize the kernel.
-//      * 
-//      * @param system     the System this kernel will be applied to
-//      * @param force      the AmoebaMultipoleForce this kernel will be used for
-//      */
-//     void initialize(const System& system, const AmoebaGeneralizedKirkwoodForce& force);
-//     /**
-//      * Execute the kernel to calculate the forces and/or energy.
-//      *
-//      * @param context        the context in which to execute this kernel
-//      * @param includeForces  true if forces should be calculated
-//      * @param includeEnergy  true if the energy should be calculated
-//      * @return the potential energy due to the force
-//      */
-//     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
-// private:
-//     System& system;
-// };
-
 /**
  * This kernel is invoked to calculate the vdw forces acting on the system and the energy of the system.
  */
@@ -482,6 +456,9 @@ public:
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 private:
     int numParticles;
+    int useNeighborList;
+    int usePBC;
+    double cutoff;
     std::vector<int> indexIVs;
     std::vector<int> indexClasses;
     std::vector< std::vector<int> > allExclusions;
@@ -491,6 +468,7 @@ private:
     std::string sigmaCombiningRule;
     std::string epsilonCombiningRule;
     System& system;
+    NeighborList* neighborList;
 };
 
 /**

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.cpp

@@ -32,12 +32,12 @@ using OpenMM::RealVec;
 
 AmoebaReferenceVdwForce::AmoebaReferenceVdwForce( ) : _nonbondedMethod(NoCutoff) {
 
+    _cutoff  = 1.0e+10;
     setSigmaCombiningRule( "ARITHMETIC" );
     setEpsilonCombiningRule( "GEOMETRIC" );
 }
 
-AmoebaReferenceVdwForce::AmoebaReferenceVdwForce( const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule, NonbondedMethod nonbondedMethod ) :
-                                                  _nonbondedMethod(nonbondedMethod) {
+AmoebaReferenceVdwForce::AmoebaReferenceVdwForce( const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule ) : _nonbondedMethod(NoCutoff) {
 
     setSigmaCombiningRule( sigmaCombiningRule );
     setEpsilonCombiningRule( epsilonCombiningRule );
@@ -51,6 +51,14 @@ void AmoebaReferenceVdwForce::setNonbondedMethod( AmoebaReferenceVdwForce::Nonbo
     _nonbondedMethod = nonbondedMethod;
 }
 
+void AmoebaReferenceVdwForce::setCutoff( double cutoff ){
+    _cutoff = cutoff;
+}
+
+double AmoebaReferenceVdwForce::getCutoff( void ) const {
+    return _cutoff;
+}
+
 void AmoebaReferenceVdwForce::setSigmaCombiningRule( const std::string& sigmaCombiningRule ){
 
     _sigmaCombiningRule = sigmaCombiningRule;
@@ -116,7 +124,7 @@ RealOpenMM AmoebaReferenceVdwForce::arithmeticEpsilonCombiningRule( RealOpenMM e
 }
 
 RealOpenMM AmoebaReferenceVdwForce::geometricEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const {
-    return 2.0*SQRT(epsilonI*epsilonJ);
+    return SQRT(epsilonI*epsilonJ);
 }
 
 RealOpenMM AmoebaReferenceVdwForce::harmonicEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const {
@@ -199,7 +207,7 @@ RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn( RealOpenMM combindedSigma,
 
 }
 
-RealOpenMM AmoebaReferenceVdwForce::calculateNoCutoffForceAndEnergy( int numParticles,
+RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergyNoCutoff( int numParticles,
                                                                      const vector<RealVec>& particlePositions,
                                                                      const std::vector<int>& indexIVs, 
                                                                      const std::vector<int>& indexClasses, 
@@ -211,7 +219,88 @@ RealOpenMM AmoebaReferenceVdwForce::calculateNoCutoffForceAndEnergy( int numPart
 
     // ---------------------------------------------------------------------------------------
 
-    //static const std::string methodName = "AmoebaReferenceVdwForce::calculateNoCutoffForceAndEnergy";
+    static const RealOpenMM zero          = 0.0;
+    static const RealOpenMM one           = 1.0;
+    static const RealOpenMM two           = 2.0;
+
+    // ---------------------------------------------------------------------------------------
+
+    // set reduced coordinates
+
+    std::vector<Vec3> reducedPositions;
+    reducedPositions.resize(numParticles);
+    for( unsigned int ii = 0; ii <  static_cast<unsigned int>(numParticles); ii++ ){
+        if( reductions[ii] != zero ){
+            int reductionIndex     = indexIVs[ii];
+            reducedPositions[ii]   = Vec3( reductions[ii]*( particlePositions[ii][0] - particlePositions[reductionIndex][0] ) + particlePositions[reductionIndex][0], 
+                                           reductions[ii]*( particlePositions[ii][1] - particlePositions[reductionIndex][1] ) + particlePositions[reductionIndex][1], 
+                                           reductions[ii]*( particlePositions[ii][2] - particlePositions[reductionIndex][2] ) + particlePositions[reductionIndex][2] ); 
+        } else {
+            reducedPositions[ii] = Vec3( particlePositions[ii][0], particlePositions[ii][1], particlePositions[ii][2] ); 
+        }
+    }
+ 
+    // loop over all ixns
+    //    (1) initialize exclusion vector
+
+    RealOpenMM energy = zero;
+    std::vector<unsigned int> exclusions(numParticles, 0);
+    for( unsigned int ii = 0; ii < static_cast<unsigned int>(numParticles); ii++ ){
+ 
+        RealOpenMM sigmaI      = sigmas[ii];
+        RealOpenMM epsilonI    = epsilons[ii];
+        for( unsigned int jj = 0; jj < allExclusions[ii].size(); jj++ ){
+            exclusions[allExclusions[ii][jj]] = 1;
+        }
+
+        for( unsigned int jj = ii+1; jj < static_cast<unsigned int>(numParticles); jj++ ){
+            if( exclusions[jj] == 0 ){
+
+                RealOpenMM combindedSigma   = (this->*_combineSigmas)(sigmaI, sigmas[jj] );
+                RealOpenMM combindedEpsilon = (this->*_combineEpsilons)(epsilonI, epsilons[jj] );
+
+                Vec3 force;
+                energy                     += calculatePairIxn( combindedSigma, combindedEpsilon,
+                                                                reducedPositions[ii], reducedPositions[jj],
+                                                                force );
+                
+                if( indexIVs[ii] == ii ){
+                    forces[ii][0] -= force[0];
+                    forces[ii][1] -= force[1];
+                    forces[ii][2] -= force[2];
+                } else {
+                    addReducedForce( ii, indexIVs[ii], reductions[ii], -one, force, forces );
+                }
+                if( indexIVs[jj] == jj ){
+                    forces[jj][0] += force[0];
+                    forces[jj][1] += force[1];
+                    forces[jj][2] += force[2];
+                } else {
+                    addReducedForce( jj, indexIVs[jj], reductions[jj], one, force, forces );
+                }
+
+            }
+        }
+
+        for( unsigned int jj = 0; jj < allExclusions[ii].size(); jj++ ){
+            exclusions[allExclusions[ii][jj]] = 0;
+        }
+    }
+
+    return energy;
+}
+
+RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergyApplyCutoff( int numParticles,
+                                                                        const vector<RealVec>& particlePositions,
+                                                                        const std::vector<int>& indexIVs, 
+                                                                        const std::vector<int>& indexClasses, 
+                                                                        const std::vector<RealOpenMM>& sigmas,
+                                                                        const std::vector<RealOpenMM>& epsilons,
+                                                                        const std::vector<RealOpenMM>& reductions,
+                                                                        const std::vector< std::vector<int> >& allExclusions,
+                                                                        vector<RealVec>& forces ) const {
+
+    // ---------------------------------------------------------------------------------------
 
     static const RealOpenMM zero          = 0.0;
     static const RealOpenMM one           = 1.0;
@@ -295,11 +384,24 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy( int numParticles,
                                                              vector<RealVec>& forces ) const {
 
     
-    if( getNonbondedMethod() == NoCutoff || 1 ){
-        return calculateNoCutoffForceAndEnergy( numParticles, particlePositions, 
+    if( getNonbondedMethod() == NoCutoff ){
+        return calculateForceAndEnergyNoCutoff( numParticles, particlePositions, 
+                                                indexIVs, indexClasses, sigmas, epsilons,
+                                                reductions, allExclusions, forces );
+    } else {
+        return calculateForceAndEnergyApplyCutoff( numParticles, particlePositions, 
                                                    indexIVs, indexClasses, sigmas, epsilons,
                                                    reductions, allExclusions, forces );
     }
 
 }    
-
+/*
+    double cutoff = force.getCutoff();
+    double taperCutoff = cutoff*0.9;
+    replacements["CUTOFF_DISTANCE"] = cu.doubleToString(force.getCutoff());
+    replacements["TAPER_CUTOFF"] = cu.doubleToString(taperCutoff);
+    replacements["TAPER_C3"] = cu.doubleToString(10/pow(taperCutoff-cutoff, 3.0));
+    replacements["TAPER_C4"] = cu.doubleToString(15/pow(taperCutoff-cutoff, 4.0));
+    replacements["TAPER_C5"] = cu.doubleToString(6/pow(taperCutoff-cutoff, 5.0));
+
+*/

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.h

@@ -79,8 +79,7 @@ public:
        --------------------------------------------------------------------------------------- */
  
     AmoebaReferenceVdwForce( const std::string& sigmaCombiningRule,
-                             const std::string& epsilonCombiningRule,
-                             NonbondedMethod nonbondedMethod );
+                             const std::string& epsilonCombiningRule );
  
     /**---------------------------------------------------------------------------------------
        
@@ -110,6 +109,26 @@ public:
     
     void setNonbondedMethod( NonbondedMethod nonbondedMethod );
 
+    /**---------------------------------------------------------------------------------------
+    
+       Get cutoff
+    
+       @return cutoff
+    
+       --------------------------------------------------------------------------------------- */
+    
+    double getCutoff( void ) const;
+
+    /**---------------------------------------------------------------------------------------
+    
+       Set cutof
+    
+       @param cutoff
+    
+       --------------------------------------------------------------------------------------- */
+    
+    void setCutoff( double cutoff );
+
     /**---------------------------------------------------------------------------------------
     
        Set sigma combining rule
@@ -181,6 +200,7 @@ private:
     std::string _sigmaCombiningRule;
     std::string _epsilonCombiningRule;
     NonbondedMethod _nonbondedMethod;
+    double _cutoff;
 
     CombiningFunction _combineSigmas;
     RealOpenMM arithmeticSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const;
@@ -246,7 +266,33 @@ private:
     
        --------------------------------------------------------------------------------------- */
     
-    RealOpenMM calculateNoCutoffForceAndEnergy( int numParticles, const std::vector<OpenMM::RealVec>& particlePositions, 
+    RealOpenMM calculateForceAndEnergyNoCutoff( int numParticles, const std::vector<OpenMM::RealVec>& particlePositions, 
+                                                const std::vector<int>& indexIVs, 
+                                                const std::vector<int>& indexClasses, 
+                                                const std::vector<RealOpenMM>& sigmas, const std::vector<RealOpenMM>& epsilons,
+                                                const std::vector<RealOpenMM>& reductions,
+                                                const std::vector< std::vector<int> >& vdwExclusions,
+                                                std::vector<OpenMM::RealVec>& forces ) const;
+         
+    /**---------------------------------------------------------------------------------------
+    
+       Calculate Vdw ixns w/ cutoff
+    
+       @param numParticles            number of particles
+       @param particlePositions       Cartesian coordinates of particles
+       @param indexIVs                position index for associated reducing particle
+       @param indexClasses            class index for combining sigmas/epsilons (not currently used)
+       @param sigmas                  particle sigmas 
+       @param epsilons                particle epsilons
+       @param reductions              particle reduction factors
+       @param vdwExclusions           particle exclusions
+       @param forces                  add forces to this vector
+    
+       @return energy
+    
+       --------------------------------------------------------------------------------------- */
+    
+    RealOpenMM calculateForceAndEnergyApplyCutoff( int numParticles, const std::vector<OpenMM::RealVec>& particlePositions, 
                                                    const std::vector<int>& indexIVs, 
                                                    const std::vector<int>& indexClasses, 
                                                    const std::vector<RealOpenMM>& sigmas, const std::vector<RealOpenMM>& epsilons,

wrappers/python/simtk/openmm/app/forcefield.py

@@ -2793,6 +2793,11 @@ class AmoebaVdwGenerator:
             else:
                 force.setCutoff(nonbondedCutoff)
 
+            # dispersion correction
+
+            if ('useDispersionCorrection' in args):
+                force.setUseDispersionCorrection(int(args['useDispersionCorrection']))
+
             if (nonbondedMethod == PME):
                 force.setPBC(1)
                 force.setUseNeighborList(1)