Implemented cutoffs and periodic boundary conditions for reference platform

olla/include/kernels.h

@@ -45,6 +45,11 @@ namespace OpenMM {
  */
 class CalcStandardMMForceFieldKernel : public KernelImpl {
 public:
+    enum NonbondedMethod {
+        NoCutoff = 0,
+        CutoffNonPeriodic = 1,
+        CutoffPeriodic = 2
+    };
     static std::string Name() {
         return "CalcStandardMMForceField";
     }
@@ -69,13 +74,17 @@ public:
      *                                  nonbonded forces.  Bonded 1-4 pairs are also included in this list, since they should be omitted from
      *                                  the standard nonbonded calculation.
      * @param nonbondedParameters       the nonbonded force parameters (charge, sigma, epsilon) for each atom
+     * @param nonbondedMethod           the method to use for handling long range nonbonded interactions
+     * @param nonbondedCutoff           the cutoff distance for nonbonded interactions (if nonbondedMethod involves a cutoff)
+     * @param periodicBoxSize           the size of the periodic box (if nonbondedMethod involves a periodic boundary conditions)
      */
     virtual void initialize(const std::vector<std::vector<int> >& bondIndices, const std::vector<std::vector<double> >& bondParameters,
             const std::vector<std::vector<int> >& angleIndices, const std::vector<std::vector<double> >& angleParameters,
             const std::vector<std::vector<int> >& periodicTorsionIndices, const std::vector<std::vector<double> >& periodicTorsionParameters,
             const std::vector<std::vector<int> >& rbTorsionIndices, const std::vector<std::vector<double> >& rbTorsionParameters,
             const std::vector<std::vector<int> >& bonded14Indices, double lj14Scale, double coulomb14Scale,
-            const std::vector<std::set<int> >& exclusions, const std::vector<std::vector<double> >& nonbondedParameters) = 0;
+            const std::vector<std::set<int> >& exclusions, const std::vector<std::vector<double> >& nonbondedParameters,
+            NonbondedMethod nonbondedMethod, double nonbondedCutoff, double periodicBoxSize[3]) = 0;
     /**
      * Execute the kernel to calculate the forces.
      * 

openmmapi/include/StandardMMForceField.h

@@ -51,6 +51,27 @@ namespace OpenMM {
 
 class OPENMM_EXPORT StandardMMForceField : public Force {
 public:
+    /**
+     * This is an enumeration of the different methods that may be used for handling long range nonbonded forces.
+     */
+    enum NonbondedMethod {
+        /**
+         * No cutoff is applied to nonbonded interactions.  The full set of N^2 interactions is computed exactly.
+         * This necessarily means that periodic boundary conditions cannot be used.
+         */
+        NoCutoff = 0,
+        /**
+         * Interactions beyond the cutoff distance are ignored.  Coulomb interactions closer than the cutoff distance
+         * are modified based using the reaction field method.
+         */
+        CutoffNonPeriodic = 1,
+        /**
+         * Periodic boundary conditions are used, so that each atom interacts only with the nearest periodic copy of
+         * each other atom.  Interactions beyond the cutoff distance are ignored.  Coulomb interactions closer than the
+         * cutoff distance are modified based using the reaction field method.
+         */
+        CutoffPeriodic = 2
+    };
     /**
      * Create a StandardMMForceField.
      * 
@@ -91,6 +112,42 @@ public:
     int getNumRBTorsions() const {
         return rbTorsions.size();
     }
+    /**
+     * Get the method used for handling long range nonbonded interactions.
+     */
+    NonbondedMethod getNonbondedMethod();
+    /**
+     * Set the method used for handling long range nonbonded interactions.
+     */
+    void setNonbondedMethod(NonbondedMethod method);
+    /**
+     * Get the cutoff distance (in nm) being used for nonbonded interactions.  If the NonbondedMethod in use
+     * does not use cutoffs, this value will have no effect.
+     */
+    double getCutoffDistance();
+    /**
+     * Set the cutoff distance (in nm) being used for nonbonded interactions.  If the NonbondedMethod in use
+     * does not use cutoffs, this value will have no effect.
+     */
+    void setCutoffDistance(double distance);
+    /**
+     * Get the dimensions of the periodic box (in nm).  If the NonbondedMethod in use does not use periodic
+     * boundary conditions, these values will have no effect.
+     *
+     * @param x      on exit, this contains the width of the periodic box along the x axis
+     * @param y      on exit, this contains the width of the periodic box along the y axis
+     * @param z      on exit, this contains the width of the periodic box along the z axis
+     */
+    void getPeriodicBoxSize(double& x, double& y, double& z);
+    /**
+     * Set the dimensions of the periodic box (in nm).  If the NonbondedMethod in use does not use periodic
+     * boundary conditions, these values will have no effect.
+     *
+     * @param x      the width of the periodic box along the x axis
+     * @param y      the width of the periodic box along the y axis
+     * @param z      the width of the periodic box along the z axis
+     */
+    void setPeriodicBoxSize(double x, double y, double z);
     /**
      * Get the nonbonded force parameters for an atom.
      * 
@@ -217,6 +274,9 @@ private:
     class AngleInfo;
     class PeriodicTorsionInfo;
     class RBTorsionInfo;
+    NonbondedMethod nonbondedMethod;
+    double cutoffDistance;
+    double periodicBoxSize[3];
 
 // Retarded visual studio compiler complains about being unable to 
 // export private stl class members.

openmmapi/src/StandardMMForceField.cpp

@@ -37,7 +37,37 @@
 using namespace OpenMM;
 
 StandardMMForceField::StandardMMForceField(int numAtoms, int numBonds, int numAngles, int numPeriodicTorsions, int numRBTorsions) :
-    atoms(numAtoms), bonds(numBonds), angles(numAngles), periodicTorsions(numPeriodicTorsions), rbTorsions(numRBTorsions) {
+        atoms(numAtoms), bonds(numBonds), angles(numAngles), periodicTorsions(numPeriodicTorsions), rbTorsions(numRBTorsions),
+        nonbondedMethod(NoCutoff), cutoffDistance(1.0) {
+    periodicBoxSize[0] = periodicBoxSize[1] = periodicBoxSize[2] = 2.0;
+}
+
+StandardMMForceField::NonbondedMethod StandardMMForceField::getNonbondedMethod() {
+    return nonbondedMethod;
+}
+
+void StandardMMForceField::setNonbondedMethod(NonbondedMethod method) {
+    nonbondedMethod = method;
+}
+
+double StandardMMForceField::getCutoffDistance() {
+    return cutoffDistance;
+}
+
+void StandardMMForceField::setCutoffDistance(double distance) {
+    cutoffDistance = distance;
+}
+
+void StandardMMForceField::getPeriodicBoxSize(double& x, double& y, double& z) {
+    x = periodicBoxSize[0];
+    y = periodicBoxSize[1];
+    z = periodicBoxSize[2];
+}
+
+void StandardMMForceField::setPeriodicBoxSize(double x, double y, double z) {
+    periodicBoxSize[0] = x;
+    periodicBoxSize[1] = y;
+    periodicBoxSize[2] = z;
 }
 
 void StandardMMForceField::getAtomParameters(int index, double& charge, double& radius, double& depth) const {

openmmapi/src/StandardMMForceFieldImpl.cpp

@@ -118,8 +118,11 @@ void StandardMMForceFieldImpl::initialize(OpenMMContextImpl& context) {
         bonded14Indices[index].push_back(iter->first);
         bonded14Indices[index++].push_back(iter->second);
     }
+    double boxSize[3];
+    owner.getPeriodicBoxSize(boxSize[0], boxSize[1], boxSize[2]);
     dynamic_cast<CalcStandardMMForceFieldKernel&>(kernel.getImpl()).initialize(bondIndices, bondParameters, angleIndices, angleParameters,
-            periodicTorsionIndices, periodicTorsionParameters, rbTorsionIndices, rbTorsionParameters, bonded14Indices, 0.5, 1.0/1.2, exclusions, nonbondedParameters);
+            periodicTorsionIndices, periodicTorsionParameters, rbTorsionIndices, rbTorsionParameters, bonded14Indices, 0.5, 1.0/1.2, exclusions,
+            nonbondedParameters, CalcStandardMMForceFieldKernel::NonbondedMethod(owner.getNonbondedMethod()), owner.getCutoffDistance(), boxSize);
 }
 
 void StandardMMForceFieldImpl::calcForces(OpenMMContextImpl& context, Stream& forces) {

platforms/reference/src/ReferenceKernels.cpp

@@ -113,6 +113,8 @@ ReferenceCalcStandardMMForceFieldKernel::~ReferenceCalcStandardMMForceFieldKerne
     disposeIntArray(exclusionArray, numAtoms);
     disposeIntArray(bonded14IndexArray, num14);
     disposeRealArray(bonded14ParamArray, num14);
+    if (neighborList != NULL)
+        delete neighborList;
 }
 
 void ReferenceCalcStandardMMForceFieldKernel::initialize(const vector<vector<int> >& bondIndices, const vector<vector<double> >& bondParameters,
@@ -120,7 +122,8 @@ void ReferenceCalcStandardMMForceFieldKernel::initialize(const vector<vector<int
         const vector<vector<int> >& periodicTorsionIndices, const vector<vector<double> >& periodicTorsionParameters,
         const vector<vector<int> >& rbTorsionIndices, const vector<vector<double> >& rbTorsionParameters,
         const vector<vector<int> >& bonded14Indices, double lj14Scale, double coulomb14Scale,
-        const vector<set<int> >& exclusions, const vector<vector<double> >& nonbondedParameters) {
+        const vector<set<int> >& exclusions, const vector<vector<double> >& nonbondedParameters,
+        NonbondedMethod nonbondedMethod, double nonbondedCutoff, double periodicBoxSize[3]) {
     numAtoms = nonbondedParameters.size();
     numBonds = bondIndices.size();
     numAngles = angleIndices.size();
@@ -142,6 +145,7 @@ void ReferenceCalcStandardMMForceFieldKernel::initialize(const vector<vector<int
         atomParamArray[i][1] = static_cast<RealOpenMM>( 2.0*sqrt(nonbondedParameters[i][2]) );
         atomParamArray[i][2] = static_cast<RealOpenMM>( nonbondedParameters[i][0]*sqrtEps );
     }
+    this->exclusions = exclusions;
     exclusionArray = new int*[numAtoms];
     for (int i = 0; i < numAtoms; ++i) {
         exclusionArray[i] = new int[exclusions[i].size()+1];
@@ -159,6 +163,16 @@ void ReferenceCalcStandardMMForceFieldKernel::initialize(const vector<vector<int
         bonded14ParamArray[i][1] = static_cast<RealOpenMM>( lj14Scale*(atomParamArray[atom1][1]*atomParamArray[atom2][1]) );
         bonded14ParamArray[i][2] = static_cast<RealOpenMM>( coulomb14Scale*(atomParamArray[atom1][2]*atomParamArray[atom2][2]) );
     }
+    this->nonbondedMethod = nonbondedMethod;
+    this->nonbondedCutoff = nonbondedCutoff;
+    this->periodicBoxSize[0] = periodicBoxSize[0];
+    this->periodicBoxSize[1] = periodicBoxSize[1];
+    this->periodicBoxSize[2] = periodicBoxSize[2];
+    if (nonbondedMethod == NoCutoff)
+        neighborList = NULL;
+    else
+        neighborList = new NeighborList();
+        
 }
 
 void ReferenceCalcStandardMMForceFieldKernel::executeForces(const Stream& positions, Stream& forces) {
@@ -174,6 +188,13 @@ void ReferenceCalcStandardMMForceFieldKernel::executeForces(const Stream& positi
     ReferenceRbDihedralBond rbTorsionBond;
     refBondForce.calculateForce(numRBTorsions, rbTorsionIndexArray, posData, rbTorsionParamArray, forceData, 0, 0, 0, rbTorsionBond);
     ReferenceLJCoulombIxn clj;
+    bool periodic = (nonbondedMethod == CutoffPeriodic);
+    if (nonbondedMethod != NoCutoff) {
+        computeNeighborListVoxelHash(*neighborList, numAtoms, posData, exclusions, periodic ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
+        clj.setUseCutoff(nonbondedCutoff, *neighborList, 78.3);
+    }
+    if (periodic)
+        clj.setPeriodic(periodicBoxSize);
     clj.calculatePairIxn(numAtoms, posData, atomParamArray, exclusionArray, 0, forceData, 0, 0);
     ReferenceLJCoulomb14 nonbonded14;
     refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, 0, 0, nonbonded14);
@@ -203,6 +224,13 @@ double ReferenceCalcStandardMMForceFieldKernel::executeEnergy(const Stream& posi
         energyArray[i] = 0;
     refBondForce.calculateForce(numRBTorsions, rbTorsionIndexArray, posData, rbTorsionParamArray, forceData, energyArray, 0, &energy, rbTorsionBond);
     ReferenceLJCoulombIxn clj;
+    bool periodic = (nonbondedMethod == CutoffPeriodic);
+    if (nonbondedMethod != NoCutoff) {
+        computeNeighborListVoxelHash(*neighborList, numAtoms, posData, exclusions, periodic ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
+        clj.setUseCutoff(nonbondedCutoff, *neighborList, 78.3);
+    }
+    if (periodic)
+        clj.setPeriodic(periodicBoxSize);
     clj.calculatePairIxn(numAtoms, posData, atomParamArray, exclusionArray, 0, forceData, 0, &energy);
     ReferenceLJCoulomb14 nonbonded14;
     for (int i = 0; i < arraySize; ++i)

platforms/reference/src/ReferenceKernels.h

@@ -34,6 +34,7 @@
 
 #include "kernels.h"
 #include "SimTKUtilities/SimTKOpenMMRealType.h"
+#include "SimTKReference/ReferenceNeighborList.h"
 
 class CpuObc;
 class ReferenceAndersenThermostat;
@@ -71,13 +72,17 @@ public:
      *                                  nonbonded forces.  Bonded 1-4 pairs are also included in this list, since they should be omitted from
      *                                  the standard nonbonded calculation.
      * @param nonbondedParameters       the nonbonded force parameters (charge, sigma, epsilon) for each atom
+     * @param nonbondedMethod           the method to use for handling long range nonbonded interactions
+     * @param nonbondedCutoff           the cutoff distance for nonbonded interactions (if nonbondedMethod involves a cutoff)
+     * @param periodicBoxSize           the size of the periodic box (if nonbondedMethod involves a periodic boundary conditions)
      */
     void initialize(const std::vector<std::vector<int> >& bondIndices, const std::vector<std::vector<double> >& bondParameters,
             const std::vector<std::vector<int> >& angleIndices, const std::vector<std::vector<double> >& angleParameters,
             const std::vector<std::vector<int> >& periodicTorsionIndices, const std::vector<std::vector<double> >& periodicTorsionParameters,
             const std::vector<std::vector<int> >& rbTorsionIndices, const std::vector<std::vector<double> >& rbTorsionParameters,
             const std::vector<std::vector<int> >& bonded14Indices, double lj14Scale, double coulomb14Scale,
-            const std::vector<std::set<int> >& exclusions, const std::vector<std::vector<double> >& nonbondedParameters);
+            const std::vector<std::set<int> >& exclusions, const std::vector<std::vector<double> >& nonbondedParameters,
+            NonbondedMethod nonbondedMethod, double nonbondedCutoff, double periodicBoxSize[3]);
     /**
      * Execute the kernel to calculate the forces.
      * 
@@ -97,6 +102,10 @@ private:
     int numAtoms, numBonds, numAngles, numPeriodicTorsions, numRBTorsions, num14;
     int **bondIndexArray, **angleIndexArray, **periodicTorsionIndexArray, **rbTorsionIndexArray, **exclusionArray, **bonded14IndexArray;
     RealOpenMM **bondParamArray, **angleParamArray, **periodicTorsionParamArray, **rbTorsionParamArray, **atomParamArray, **bonded14ParamArray;
+    RealOpenMM nonbondedCutoff, periodicBoxSize[3];
+    std::vector<std::set<int> > exclusions;
+    NonbondedMethod nonbondedMethod;
+    NeighborList* neighborList;
 };
 
 /**

platforms/reference/src/SimTKReference/ReferenceForce.cpp

@@ -62,6 +62,20 @@ ReferenceForce::~ReferenceForce( ){
 
 }
 
+/**---------------------------------------------------------------------------------------
+
+   Given two coordinates on a periodic lattice, return the difference between them.
+
+   --------------------------------------------------------------------------------------- */
+
+RealOpenMM ReferenceForce::periodicDifference(RealOpenMM val1, RealOpenMM val2, RealOpenMM period) {
+    RealOpenMM diff = val1-val2;
+    RealOpenMM base = floor(diff/period+0.5)*period;
+    return diff-base;
+
+}
+
+
 /**---------------------------------------------------------------------------------------
 
    Get deltaR and distance and distance**2 between atomI and atomJ (static method)
@@ -94,6 +108,39 @@ int ReferenceForce::getDeltaR( const RealOpenMM* atomCoordinatesI, const RealOpe
    return ReferenceForce::DefaultReturn;
 }
 
+/**---------------------------------------------------------------------------------------
+
+   Get deltaR and distance and distance**2 between atomI and atomJ, assuming periodic
+   boundary conditions (static method); deltaR: j - i
+
+   @param atomCoordinatesI    atom i coordinates
+   @param atomCoordinatesI    atom j coordinates
+   @param boxSize             X, Y, and Z sizes of the periodic box
+   @param deltaR              deltaX, deltaY, deltaZ, R2, R upon return
+
+   @return ReferenceForce::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceForce::getDeltaRPeriodic( const RealOpenMM* atomCoordinatesI, const RealOpenMM* atomCoordinatesJ,
+                               const RealOpenMM* boxSize, RealOpenMM* deltaR ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName = "\nReferenceForce::getDeltaR";
+
+   // ---------------------------------------------------------------------------------------
+
+   deltaR[XIndex]    = periodicDifference(atomCoordinatesJ[0], atomCoordinatesI[0], boxSize[0]);
+   deltaR[YIndex]    = periodicDifference(atomCoordinatesJ[1], atomCoordinatesI[1], boxSize[1]);
+   deltaR[ZIndex]    = periodicDifference(atomCoordinatesJ[2], atomCoordinatesI[2], boxSize[2]);
+
+   deltaR[R2Index]   = DOT3( deltaR, deltaR );
+   deltaR[RIndex]    = (RealOpenMM) SQRT( deltaR[R2Index] );
+
+   return ReferenceForce::DefaultReturn;
+}
+
 /**---------------------------------------------------------------------------------------
 
    Get deltaR between atomI and atomJ (static method); deltaR: j - i

platforms/reference/src/SimTKReference/ReferenceForce.h

@@ -31,6 +31,8 @@ class ReferenceForce {
 
    private:
        
+       static RealOpenMM periodicDifference(RealOpenMM val1, RealOpenMM val2, RealOpenMM period);
+
    public:
 
       /**---------------------------------------------------------------------------------------
@@ -83,6 +85,23 @@ class ReferenceForce {
       static int getDeltaR( const RealOpenMM* atomCoordinatesI, const RealOpenMM* atomCoordinatesJ,
                             RealOpenMM* deltaR );
       
+      /**---------------------------------------------------------------------------------------
+
+         Get deltaR and distance and distance**2 between atomI and atomJ, assuming periodic
+         boundary conditions (static method); deltaR: j - i
+
+         @param atomCoordinatesI    atom i coordinates
+         @param atomCoordinatesI    atom j coordinates
+         @param boxSize             X, Y, and Z sizes of the periodic box
+         @param deltaR              deltaX, deltaY, deltaZ, R2, R upon return
+
+         @return ReferenceForce::DefaultReturn
+
+         --------------------------------------------------------------------------------------- */
+
+      static int ReferenceForce::getDeltaRPeriodic( const RealOpenMM* atomCoordinatesI, const RealOpenMM* atomCoordinatesJ,
+                                             const RealOpenMM* boxSize, RealOpenMM* deltaR );
+
     /**---------------------------------------------------------------------------------------
       
          Get deltaR between atomI and atomJ (static method): deltaR: j - i

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp

@@ -37,7 +37,7 @@
 
    --------------------------------------------------------------------------------------- */
 
-ReferenceLJCoulombIxn::ReferenceLJCoulombIxn( ){
+ReferenceLJCoulombIxn::ReferenceLJCoulombIxn( ) : cutoff(false), periodic(false) {
 
    // ---------------------------------------------------------------------------------------
 
@@ -63,6 +63,55 @@ ReferenceLJCoulombIxn::~ReferenceLJCoulombIxn( ){
 
 }
 
+  /**---------------------------------------------------------------------------------------
+
+     Set the force to use a cutoff.
+
+     @param distance            the cutoff distance
+     @param neighbors           the neighbor list to use
+     @param solventDielectric   the dielectric constant of the bulk solvent
+
+     @return ReferenceForce::DefaultReturn
+
+     --------------------------------------------------------------------------------------- */
+
+  int ReferenceLJCoulombIxn::setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors, RealOpenMM solventDielectric ) {
+    
+    cutoff = true;
+    cutoffDistance = distance;
+    neighborList = &neighbors;
+    krf = pow(cutoffDistance, -3.0)*(solventDielectric-1.0)/(2.0*solventDielectric+1.0);
+    crf = (1.0/cutoffDistance)*(3.0*solventDielectric)/(2.0*solventDielectric+1.0);
+            
+    return ReferenceForce::DefaultReturn;
+  }
+
+  /**---------------------------------------------------------------------------------------
+
+     Set the force to use periodic boundary conditions.  This requires that a cutoff has
+     also been set, and the smallest side of the periodic box is at least twice the cutoff
+     distance.
+
+     @param boxSize             the X, Y, and Z widths of the periodic box
+
+     @return ReferenceForce::DefaultReturn
+
+     --------------------------------------------------------------------------------------- */
+
+  int ReferenceLJCoulombIxn::setPeriodic( RealOpenMM* boxSize ) {
+
+    assert(cutoff);
+    assert(boxSize[0] >= 2.0*cutoffDistance);
+    assert(boxSize[1] >= 2.0*cutoffDistance);
+    assert(boxSize[2] >= 2.0*cutoffDistance);
+    periodic = true;
+    periodicBoxSize[0] = boxSize[0];
+    periodicBoxSize[1] = boxSize[1];
+    periodicBoxSize[2] = boxSize[2];
+    return ReferenceForce::DefaultReturn;
+
+  }
+
 /**---------------------------------------------------------------------------------------
 
    Calculate parameters for LJ Coulomb ixn
@@ -140,13 +189,69 @@ int ReferenceLJCoulombIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM** ato
                                              RealOpenMM* fixedParameters, RealOpenMM** forces,
                                              RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
 
-   // ---------------------------------------------------------------------------------------
+   if (cutoff) {
+       for (int i = 0; i < neighborList->size(); i++) {
+           OpenMM::AtomPair pair = (*neighborList)[i];
+           calculateOneIxn(pair.first, pair.second, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy);
+       }
+   }
+   else {
+       // allocate and initialize exclusion array
 
-   // static const char* methodName = "\nReferenceLJCoulombIxn::calculatePairIxn";
+       int* exclusionIndices = new int[numberOfAtoms];
+       for( int ii = 0; ii < numberOfAtoms; ii++ ){
+          exclusionIndices[ii] = -1;
+       }
+
+       for( int ii = 0; ii < numberOfAtoms; ii++ ){
+
+          // set exclusions
+
+          for( int jj = 1; jj <= exclusions[ii][0]; jj++ ){
+             exclusionIndices[exclusions[ii][jj]] = ii;
+          }
+
+          // loop over atom pairs
+
+          for( int jj = ii+1; jj < numberOfAtoms; jj++ ){
+
+             if( exclusionIndices[jj] != ii ){
+                 calculateOneIxn(ii, jj, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy);
+             }
+          }
+       }
+
+       delete[] exclusionIndices;
+   }
+
+   return ReferenceForce::DefaultReturn;
+}
+
+  /**---------------------------------------------------------------------------------------
+
+     Calculate LJ Coulomb pair ixn between two atoms
+
+     @param ii               the index of the first atom
+     @param jj               the index of the second atom
+     @param atomCoordinates  atom coordinates
+     @param atomParameters   atom parameters (charges, c6, c12, ...)     atomParameters[atomIndex][paramterIndex]
+     @param forces           force array (forces added)
+     @param energyByAtom     atom energy
+     @param totalEnergy      total energy
+
+     @return ReferenceForce::DefaultReturn
+
+     --------------------------------------------------------------------------------------- */
+
+  int ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, RealOpenMM** atomCoordinates,
+                        RealOpenMM** atomParameters, RealOpenMM** forces,
+                        RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
 
     // ---------------------------------------------------------------------------------------
 
-   static const std::string methodName = "\nReferenceLJCoulombIxn::calculatePairIxn";
+    static const std::string methodName = "\nReferenceLJCoulombIxn::calculateOneIxn";
+
+    // ---------------------------------------------------------------------------------------
 
     // constants -- reduce Visual Studio warnings regarding conversions between float & double
 
@@ -168,33 +273,14 @@ int ReferenceLJCoulombIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM** ato
 
     RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
 
-   // ---------------------------------------------------------------------------------------
-
-   // allocate and initialize exclusion array
-
-   int* exclusionIndices = new int[numberOfAtoms];
-   for( int ii = 0; ii < numberOfAtoms; ii++ ){
-      exclusionIndices[ii] = -1;
-   }
-
-   for( int ii = 0; ii < numberOfAtoms; ii++ ){
-
-      // set exclusions
-
-      for( int jj = 1; jj <= exclusions[ii][0]; jj++ ){
-         exclusionIndices[exclusions[ii][jj]] = ii;
-      }
-     
-      // loop over atom pairs
-
-      for( int jj = ii+1; jj < numberOfAtoms; jj++ ){
-
-         if( exclusionIndices[jj] != ii ){
-
     // get deltaR, R2, and R between 2 atoms
 
+    if (periodic)
+        ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );  
+    else
         ReferenceForce::getDeltaR( atomCoordinates[jj], atomCoordinates[ii], deltaR[0] );  
 
+    RealOpenMM r2        = deltaR[0][ReferenceForce::R2Index];
     RealOpenMM inverseR  = one/(deltaR[0][ReferenceForce::RIndex]);
     RealOpenMM sig       = atomParameters[ii][SigIndex] +  atomParameters[jj][SigIndex];
     RealOpenMM sig2      = inverseR*sig;
@@ -203,6 +289,9 @@ int ReferenceLJCoulombIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM** ato
 
     RealOpenMM eps       = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
     RealOpenMM dEdR      = eps*( twelve*sig6 - six )*sig6;
+               if (cutoff)
+                   dEdR += atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR-2.0*krf*r2);
+               else
                    dEdR += atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
                dEdR     *= inverseR*inverseR;
          
@@ -219,6 +308,9 @@ int ReferenceLJCoulombIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM** ato
     // accumulate energies
 
     if( totalEnergy || energyByAtom ) {
+        if (cutoff)
+            energy = atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR+krf*r2-crf);
+        else
             energy = atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
         energy += eps*(sig6-one)*sig6;
         if( totalEnergy )
@@ -282,10 +374,4 @@ int ReferenceLJCoulombIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM** ato
        SimTKOpenMMLog::printMessage( message );
     }
   }
-      }
-   }
-
-   delete[] exclusionIndices;
 
-   return ReferenceForce::DefaultReturn;
-}

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.h

@@ -26,6 +26,7 @@
 #define __ReferenceLJCoulombIxn_H__
 
 #include "ReferencePairIxn.h"
+#include "ReferenceNeighborList.h"
 
 // ---------------------------------------------------------------------------------------
 
@@ -33,12 +34,40 @@ class ReferenceLJCoulombIxn : public ReferencePairIxn {
 
    private:
        
+      bool cutoff;
+      bool periodic;
+      const OpenMM::NeighborList* neighborList;
+      RealOpenMM periodicBoxSize[3];
+      RealOpenMM cutoffDistance;
+      RealOpenMM krf, crf;
+
       // parameter indices
 
       static const int SigIndex = 0;
       static const int EpsIndex = 1;
       static const int   QIndex = 2;
             
+      /**---------------------------------------------------------------------------------------
+      
+         Calculate LJ Coulomb pair ixn between two atoms
+      
+         @param atom1            the index of the first atom
+         @param atom2            the index of the second atom
+         @param atomCoordinates  atom coordinates
+         @param atomParameters   atom parameters (charges, c6, c12, ...)     atomParameters[atomIndex][paramterIndex]
+         @param forces           force array (forces added)
+         @param energyByAtom     atom energy
+         @param totalEnergy      total energy
+      
+         @return ReferenceForce::DefaultReturn
+            
+         --------------------------------------------------------------------------------------- */
+          
+      int calculateOneIxn( int atom1, int atom2, RealOpenMM** atomCoordinates,
+                            RealOpenMM** atomParameters, RealOpenMM** forces,
+                            RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
+
+
    public:
 
       /**---------------------------------------------------------------------------------------
@@ -57,6 +86,34 @@ class ReferenceLJCoulombIxn : public ReferencePairIxn {
 
        ~ReferenceLJCoulombIxn( );
 
+      /**---------------------------------------------------------------------------------------
+      
+         Set the force to use a cutoff.
+      
+         @param distance            the cutoff distance
+         @param neighbors           the neighbor list to use
+         @param solventDielectric   the dielectric constant of the bulk solvent
+      
+         @return ReferenceForce::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors, RealOpenMM solventDielectric );
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Set the force to use periodic boundary conditions.  This requires that a cutoff has
+         already been set, and the smallest side of the periodic box is at least twice the cutoff
+         distance.
+      
+         @param boxSize             the X, Y, and Z widths of the periodic box
+      
+         @return ReferenceForce::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int setPeriodic( RealOpenMM* boxSize );      
+
       /**---------------------------------------------------------------------------------------
       
          Calculate parameters for LJ 1-4 ixn

platforms/reference/src/SimTKReference/NeighborList.cpp → platforms/reference/src/SimTKReference/ReferenceNeighborList.cpp

-#include "NeighborList.h"
+#include "ReferenceNeighborList.h"
 #include <set>
 #include <map>
 #include <cmath>
@@ -10,11 +10,25 @@ namespace OpenMM {
 
 typedef std::vector<AtomIndex> AtomList;
 
+static double periodicDifference(double val1, double val2, double period) {
+    double diff = val1-val2;
+    double base = floor(diff/period+0.5)*period;
+    return diff-base;
+}
+
 // squared distance between two points
-double compPairDistanceSquared(const Vec3& pos1, const Vec3& pos2) {
-    double dx = pos2[0] - pos1[0];
-    double dy = pos2[1] - pos1[1];
-    double dz = pos2[2] - pos1[2];
+static double compPairDistanceSquared(const RealOpenMM* pos1, const RealOpenMM* pos2, const RealOpenMM* periodicBoxSize) {
+    double dx, dy, dz;
+    if (periodicBoxSize == NULL) {
+        dx = pos2[0] - pos1[0];
+        dy = pos2[1] - pos1[1];
+        dz = pos2[2] - pos1[2];
+    }
+    else {
+        dx = periodicDifference(pos2[0], pos1[0], periodicBoxSize[0]);
+        dy = periodicDifference(pos2[1], pos1[1], periodicBoxSize[1]);
+        dz = periodicDifference(pos2[2], pos1[2], periodicBoxSize[2]);
+    }
     return dx*dx + dy*dy + dz*dz;
 }
 
@@ -22,7 +36,10 @@ double compPairDistanceSquared(const Vec3& pos1, const Vec3& pos2) {
 // for pedagogical purposes and simplicity
 void computeNeighborListNaive(
                               NeighborList& neighborList,
+                              int nAtoms,
                               const AtomLocationList& atomLocations, 
+                              const vector<set<int> >& exclusions,
+                              const RealOpenMM* periodicBoxSize,
                               double maxDistance,
                               double minDistance,
                               bool reportSymmetricPairs
@@ -30,8 +47,6 @@ void computeNeighborListNaive(
 {
     neighborList.clear();
     
-    int nAtoms = atomLocations.size();
-
     double maxDistanceSquared = maxDistance * maxDistance;
     double minDistanceSquared = minDistance * minDistance;
 
@@ -39,8 +54,9 @@ void computeNeighborListNaive(
     {
         for (AtomIndex atomJ = atomI + 1; atomJ < nAtoms; ++atomJ)
         {
-            double pairDistanceSquared = compPairDistanceSquared(atomLocations[atomI], atomLocations[atomJ]);
-            if ( (pairDistanceSquared <= maxDistanceSquared)  && (pairDistanceSquared >= minDistanceSquared) )
+            double pairDistanceSquared = compPairDistanceSquared(atomLocations[atomI], atomLocations[atomJ], periodicBoxSize);
+            if ( (pairDistanceSquared <= maxDistanceSquared)  && (pairDistanceSquared >= minDistanceSquared))
+                if (exclusions[atomI].find(atomJ) == exclusions[atomI].end())
                 {
                     neighborList.push_back( AtomPair(atomI, atomJ) );
                     if (reportSymmetricPairs)
@@ -72,15 +88,22 @@ public:
 };
 
 
-typedef std::pair<Vec3, AtomIndex> VoxelItem;
+typedef std::pair<const RealOpenMM*, AtomIndex> VoxelItem;
 typedef std::vector< VoxelItem > Voxel;
 
 class VoxelHash
 {
 public:
-    VoxelHash(double vs) : voxelSize(vs) {}
+    VoxelHash(double vsx, double vsy, double vsz, const RealOpenMM* periodicBoxSize) :
+            voxelSizeX(vsx), voxelSizeY(vsy), voxelSizeZ(vsz), periodicBoxSize(periodicBoxSize) {
+        if (periodicBoxSize != NULL) {
+            nx = floor(periodicBoxSize[0]/voxelSizeX+0.5);
+            ny = floor(periodicBoxSize[1]/voxelSizeY+0.5);
+            nz = floor(periodicBoxSize[2]/voxelSizeZ+0.5);
+        }
+    }
 
-    void insert(const AtomIndex& item, const Vec3& location)
+    void insert(const AtomIndex& item, const RealOpenMM* location)
     {
         VoxelIndex voxelIndex = getVoxelIndex(location);
         if ( voxelMap.find(voxelIndex) == voxelMap.end() ) voxelMap[voxelIndex] = Voxel(); 
@@ -89,10 +112,21 @@ public:
     }
 
 
-    VoxelIndex getVoxelIndex(const Vec3& location) const {
-        int x = int(floor(location[0] / voxelSize));
-        int y = int(floor(location[1] / voxelSize));
-        int z = int(floor(location[2] / voxelSize));
+    VoxelIndex getVoxelIndex(const RealOpenMM* location) const {
+        double xperiodic, yperiodic, zperiodic;
+        if (periodicBoxSize == NULL) {
+            xperiodic = location[0];
+            yperiodic = location[1];
+            zperiodic = location[2];
+        }
+        else {
+            xperiodic = location[0]-periodicBoxSize[0]*floor(location[0]/periodicBoxSize[0]);
+            yperiodic = location[1]-periodicBoxSize[1]*floor(location[1]/periodicBoxSize[1]);
+            zperiodic = location[2]-periodicBoxSize[2]*floor(location[2]/periodicBoxSize[2]);
+        }
+        int x = int(floor(xperiodic / voxelSizeX));
+        int y = int(floor(yperiodic / voxelSizeY));
+        int z = int(floor(zperiodic / voxelSizeZ));
         
         return VoxelIndex(x, y, z);
     }
@@ -100,6 +134,7 @@ public:
     void getNeighbors(
             NeighborList& neighbors, 
             const VoxelItem& referencePoint, 
+            const vector<set<int> >& exclusions,
             bool reportSymmetricPairs,
             double maxDistance, 
             double minDistance) const 
@@ -109,35 +144,54 @@ public:
         // TODO use more clever selection of neighboring voxels
         assert(maxDistance > 0);
         assert(minDistance >= 0);
-        assert(voxelSize > 0);
+        assert(voxelSizeX > 0);
+        assert(voxelSizeY > 0);
+        assert(voxelSizeZ > 0);
 
         const AtomIndex atomI = referencePoint.second;
-        const Vec3& locationI = referencePoint.first;
+        const RealOpenMM* locationI = referencePoint.first;
         
         double maxDistanceSquared = maxDistance * maxDistance;
         double minDistanceSquared = minDistance * minDistance;
 
-        int dIndex = int(maxDistance / voxelSize) + 1; // How may voxels away do we have to look?
+        int dIndexX = int(maxDistance / voxelSizeX) + 1; // How may voxels away do we have to look?
+        int dIndexY = int(maxDistance / voxelSizeY) + 1;
+        int dIndexZ = int(maxDistance / voxelSizeZ) + 1;
         VoxelIndex centerVoxelIndex = getVoxelIndex(locationI);
-        for (int x = centerVoxelIndex.x - dIndex; x <= centerVoxelIndex.x + dIndex; ++x)
+        int lastx = centerVoxelIndex.x+dIndexX;
+        int lasty = centerVoxelIndex.y+dIndexY;
+        int lastz = centerVoxelIndex.z+dIndexZ;
+        if (periodicBoxSize != NULL) {
+            lastx = min(lastx, centerVoxelIndex.x-dIndexX+nx-1);
+            lasty = min(lasty, centerVoxelIndex.y-dIndexY+ny-1);
+            lastz = min(lastz, centerVoxelIndex.z-dIndexZ+nz-1);
+        }
+        for (int x = centerVoxelIndex.x - dIndexX; x <= lastx; ++x)
         {
-            for (int y = centerVoxelIndex.y - dIndex; y <= centerVoxelIndex.y + dIndex; ++y)
+            for (int y = centerVoxelIndex.y - dIndexY; y <= lasty; ++y)
             {
-                for (int z = centerVoxelIndex.z - dIndex; z <= centerVoxelIndex.z + dIndex; ++z)
+                for (int z = centerVoxelIndex.z - dIndexZ; z <= lastz; ++z)
                 {
                     VoxelIndex voxelIndex(x, y, z);
+                    if (periodicBoxSize != NULL) {
+                        voxelIndex.x = (x+nx)%nx;
+                        voxelIndex.y = (y+ny)%ny;
+                        voxelIndex.z = (z+nz)%nz;
+                    }
                     if (voxelMap.find(voxelIndex) == voxelMap.end()) continue; // no such voxel; skip
                     const Voxel& voxel = voxelMap.find(voxelIndex)->second;
                     for (Voxel::const_iterator itemIter = voxel.begin(); itemIter != voxel.end(); ++itemIter)
                     {
                         const AtomIndex atomJ = itemIter->second;
-                        const Vec3& locationJ = itemIter->first;
+                        const RealOpenMM* locationJ = itemIter->first;
                         
                         // Ignore self hits
                         if (atomI == atomJ) continue;
                         
-                        Vec3 dv = locationI - locationJ;
-                        double dSquared = dv.dot(dv);
+                        // Ignore exclusions.
+                        if (exclusions[atomI].find(atomJ) != exclusions[atomI].end()) continue;
+                        
+                        double dSquared = compPairDistanceSquared(locationI, locationJ, periodicBoxSize);
                         if (dSquared > maxDistanceSquared) continue;
                         if (dSquared < minDistanceSquared) continue;
                         
@@ -151,7 +205,9 @@ public:
     }
 
 private:
-    double voxelSize;
+    double voxelSizeX, voxelSizeY, voxelSizeZ;
+    int nx, ny, nz;
+    const RealOpenMM* periodicBoxSize;
     std::map<VoxelIndex, Voxel> voxelMap;
 };
 
@@ -159,7 +215,10 @@ private:
 // O(n) neighbor list method using voxel hash data structure
 void computeNeighborListVoxelHash(
                               NeighborList& neighborList,
+                              int nAtoms,
                               const AtomLocationList& atomLocations, 
+                              const vector<set<int> >& exclusions,
+                              const RealOpenMM* periodicBoxSize,
                               double maxDistance,
                               double minDistance,
                               bool reportSymmetricPairs
@@ -167,17 +226,23 @@ void computeNeighborListVoxelHash(
 {
     neighborList.clear();
 
-    const int nAtoms = atomLocations.size();
-
-    const double edgeSize = maxDistance; // TODO - adjust this as needed
-    VoxelHash voxelHash(edgeSize);
+    double edgeSizeX, edgeSizeY, edgeSizeZ;
+    if (periodicBoxSize == NULL)
+        edgeSizeX = edgeSizeY = edgeSizeZ = maxDistance; // TODO - adjust this as needed
+    else {
+        edgeSizeX = periodicBoxSize[0]/floor(periodicBoxSize[0]/maxDistance);
+        edgeSizeY = periodicBoxSize[1]/floor(periodicBoxSize[1]/maxDistance);
+        edgeSizeZ = periodicBoxSize[2]/floor(periodicBoxSize[2]/maxDistance);
+    }
+    VoxelHash voxelHash(edgeSizeX, edgeSizeY, edgeSizeZ, periodicBoxSize);
     for (AtomIndex atomJ = 0; atomJ < nAtoms; ++atomJ) // use "j", because j > i for pairs
     {
         // 1) Find other atoms that are close to this one
-        const Vec3& location = atomLocations[atomJ];
+        const RealOpenMM* location = atomLocations[atomJ];
         voxelHash.getNeighbors( 
             neighborList, 
             VoxelItem(location, atomJ),
+            exclusions,
             reportSymmetricPairs, 
             maxDistance, 
             minDistance);

platforms/reference/src/SimTKReference/NeighborList.h → platforms/reference/src/SimTKReference/ReferenceNeighborList.h

 #ifndef OPENMM_REFERENCE_NEIGHBORLIST_H_
 #define OPENMM_REFERENCE_NEIGHBORLIST_H_
 
-#include "Vec3.h"
+#include "../SimTKUtilities/SimTKOpenMMRealType.h"
+#include <set>
 #include <vector>
 
 namespace OpenMM {
 
-typedef std::vector<Vec3> AtomLocationList;
+typedef RealOpenMM** AtomLocationList;
 typedef unsigned int AtomIndex;
 typedef std::pair<AtomIndex, AtomIndex> AtomPair;
 typedef std::vector<AtomPair>  NeighborList;
@@ -17,7 +18,10 @@ typedef std::vector<AtomPair>  NeighborList;
 // neighbors are added
 void computeNeighborListNaive(
                               NeighborList& neighborList,
+                              int nAtoms,
                               const AtomLocationList& atomLocations, 
+                              const std::vector<std::set<int> >& exclusions,
+                              const RealOpenMM* periodicBoxSize,
                               double maxDistance,
                               double minDistance = 0.0,
                               bool reportSymmetricPairs = false
@@ -28,9 +32,12 @@ void computeNeighborListNaive(
 // neighbors are added
 void computeNeighborListVoxelHash(
                               NeighborList& neighborList,
+                              int nAtoms,
                               const AtomLocationList& atomLocations, 
+                              const std::vector<std::set<int> >& exclusions,
+                              const RealOpenMM* periodicBoxSize,
                               double maxDistance,
-                              double minDistance,
+                              double minDistance = 0.0,
                               bool reportSymmetricPairs = false
                              );
 

platforms/reference/tests/TestNeighborList.cpp

-#include "../src/SimTKReference/NeighborList.h"
-#include <cassert>
-#include <iostream>
-
-using namespace std;
-using namespace OpenMM;
-
-void testNeighborList()
-{
-    AtomLocationList atomList;
-    atomList.push_back(Vec3(13.6, 0, 0));
-    atomList.push_back(Vec3(0, 0, 0));
-    
-    NeighborList neighborList;
-    
-    computeNeighborListNaive(neighborList, atomList, 13.7, 0.01);
-    assert(neighborList.size() == 1);
-    
-    computeNeighborListNaive(neighborList, atomList, 13.5, 0.01);
-    assert(neighborList.size() == 0);
-    
-    computeNeighborListVoxelHash(neighborList, atomList, 13.7, 0.01);
-    assert(neighborList.size() == 1);
-    
-    computeNeighborListVoxelHash(neighborList, atomList, 13.5, 0.01);
-    assert(neighborList.size() == 0);
-}
-
-int main() 
-{
-try {
-    testNeighborList();
-    
-    cout << "Test Passed" << endl;
-    return 0;
-}
-catch (...) {
-    cerr << "*** ERROR: Test Failed ***" << endl;
-    return 1;
-}
-}
-

platforms/reference/tests/TestReferenceNeighborList.cpp

+#include "../../../tests/AssertionUtilities.h"
+#include "../src/SimTKReference/ReferenceNeighborList.h"
+#include "../src/sfmt/SFMT.h"
+#include <cassert>
+#include <iostream>
+
+using namespace std;
+using namespace OpenMM;
+
+void testNeighborList()
+{
+    RealOpenMM* atomList[2];
+    atomList[0] = new RealOpenMM[3];
+    atomList[1] = new RealOpenMM[3];
+    atomList[2] = new RealOpenMM[3];
+    atomList[0][0] = 13.6;
+    atomList[0][1] = 0;
+    atomList[0][2] = 0;
+    atomList[1][0] = 0;
+    atomList[1][1] = 0;
+    atomList[1][2] = 0;
+    vector<set<int> > exclusions(2);
+    
+    NeighborList neighborList;
+    
+    computeNeighborListNaive(neighborList, 2, atomList, exclusions, NULL, 13.7, 0.01);
+    assert(neighborList.size() == 1);
+    
+    computeNeighborListNaive(neighborList, 2, atomList, exclusions, NULL, 13.5, 0.01);
+    assert(neighborList.size() == 0);
+    
+    computeNeighborListVoxelHash(neighborList, 2, atomList, exclusions, NULL, 13.7, 0.01);
+    assert(neighborList.size() == 1);
+    
+    computeNeighborListVoxelHash(neighborList, 2, atomList, exclusions, NULL, 13.5, 0.01);
+    assert(neighborList.size() == 0);
+    
+    delete[] atomList[0];
+    delete[] atomList[1];
+    delete[] atomList[2];
+}
+
+double periodicDifference(double val1, double val2, double period) {
+    double diff = val1-val2;
+    double base = floor(diff/period+0.5)*period;
+    return diff-base;
+}
+
+double distance2(RealOpenMM* pos1, RealOpenMM* pos2, const RealOpenMM* periodicBoxSize) {
+    double dx = periodicDifference(pos1[0], pos2[0], periodicBoxSize[0]);
+    double dy = periodicDifference(pos1[1], pos2[1], periodicBoxSize[1]);
+    double dz = periodicDifference(pos1[2], pos2[2], periodicBoxSize[2]);
+    return dx*dx+dy*dy+dz*dz;
+}
+
+void verifyNeighborList(NeighborList& list, int numAtoms, RealOpenMM** positions, const RealOpenMM* periodicBoxSize, double cutoff) {
+    for (int i = 0; i < list.size(); i++) {
+        int atom1 = list[i].first;
+        int atom2 = list[i].second;
+        ASSERT(distance2(positions[atom1], positions[atom2], periodicBoxSize) <= cutoff*cutoff);
+    }
+    int count = 0;
+    for (int i = 0; i < numAtoms; i++)
+        for (int j = i+1; j < numAtoms; j++)
+            if (distance2(positions[i], positions[j], periodicBoxSize) <= cutoff*cutoff)
+                count++;
+    ASSERT(count == list.size());
+}
+
+void testPeriodic() {
+    const int numAtoms = 100;
+    const double cutoff = 3.0;
+    const RealOpenMM periodicBoxSize[3] = {20.0, 15.0, 22.0};
+    RealOpenMM* atomList[numAtoms];
+    init_gen_rand(0);
+    for (int i = 0; i <numAtoms; i++) {
+        atomList[i] = new RealOpenMM[3];
+        atomList[i][0] = genrand_real2()*periodicBoxSize[0]*3;
+        atomList[i][1] = genrand_real2()*periodicBoxSize[1]*3;
+        atomList[i][2] = genrand_real2()*periodicBoxSize[2]*3;
+    }
+    vector<set<int> > exclusions(numAtoms);
+    NeighborList neighborList;
+    computeNeighborListNaive(neighborList, numAtoms, atomList, exclusions, periodicBoxSize, cutoff);
+    verifyNeighborList(neighborList, numAtoms, atomList, periodicBoxSize, cutoff);
+    computeNeighborListVoxelHash(neighborList, numAtoms, atomList, exclusions, periodicBoxSize, cutoff);
+    verifyNeighborList(neighborList, numAtoms, atomList, periodicBoxSize, cutoff);
+    for (int i = 0; i <numAtoms; i++)
+        delete[] atomList[i];
+}
+
+int main() 
+{
+try {
+    testNeighborList();
+    testPeriodic();
+    
+    cout << "Test Passed" << endl;
+    return 0;
+}
+catch (...) {
+    cerr << "*** ERROR: Test Failed ***" << endl;
+    return 1;
+}
+}
+

platforms/reference/tests/TestReferenceStandardMMForceField.cpp

@@ -265,6 +265,71 @@ void testExclusionsAnd14() {
     }
 }
 
+void testCutoff() {
+    ReferencePlatform platform;
+    System system(3, 0);
+    VerletIntegrator integrator(0.01);
+    StandardMMForceField* forceField = new StandardMMForceField(3, 0, 0, 0, 0);
+    forceField->setAtomParameters(0, 1.0, 1, 0);
+    forceField->setAtomParameters(1, 1.0, 1, 0);
+    forceField->setAtomParameters(2, 1.0, 1, 0);
+    forceField->setNonbondedMethod(StandardMMForceField::CutoffNonPeriodic);
+    const double cutoff = 2.9;
+    forceField->setCutoffDistance(cutoff);
+    system.addForce(forceField);
+    OpenMMContext context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(0, 2, 0);
+    positions[2] = Vec3(0, 3, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    const double eps = 78.3;
+    const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
+    const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
+    const double force1 = 138.935485*(1.0)*(0.25-2.0*krf*2.0);
+    const double force2 = 138.935485*(1.0)*(1.0-2.0*krf*1.0);
+    ASSERT_EQUAL_VEC(Vec3(0, -force1, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, force1-force2, 0), forces[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, force2, 0), forces[2], TOL);
+    const double energy1 = 138.935485*(1.0)*(0.5+krf*4.0-crf);
+    const double energy2 = 138.935485*(1.0)*(1.0+krf*1.0-crf);
+    ASSERT_EQUAL_TOL(energy1+energy2, state.getPotentialEnergy(), TOL);
+}
+
+void testPeriodic() {
+    ReferencePlatform platform;
+    System system(3, 0);
+    VerletIntegrator integrator(0.01);
+    StandardMMForceField* forceField = new StandardMMForceField(3, 1, 0, 0, 0);
+    forceField->setAtomParameters(0, 1.0, 1, 0);
+    forceField->setAtomParameters(1, 1.0, 1, 0);
+    forceField->setAtomParameters(2, 1.0, 1, 0);
+    forceField->setBondParameters(0, 0, 1, 1.0, 0.0);
+    forceField->setNonbondedMethod(StandardMMForceField::CutoffPeriodic);
+    const double cutoff = 2.0;
+    forceField->setCutoffDistance(cutoff);
+    forceField->setPeriodicBoxSize(4.0, 4.0, 4.0);
+    system.addForce(forceField);
+    OpenMMContext context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(2, 0, 0);
+    positions[2] = Vec3(3, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    const double eps = 78.3;
+    const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
+    const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
+    const double force = 138.935485*(1.0)*(1.0-2.0*krf*1.0);
+    ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL);
+    ASSERT_EQUAL_TOL(2*138.935485*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
+}
+
 int main() {
     try {
         testBonds();
@@ -274,6 +339,8 @@ int main() {
         testCoulomb();
         testLJ();
         testExclusionsAnd14();
+        testCutoff();
+        testPeriodic();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

tests/TestFindExclusions.cpp

@@ -145,7 +145,8 @@ public:
             const vector<vector<int> >& periodicTorsionIndices, const vector<vector<double> >& periodicTorsionParameters,
             const vector<vector<int> >& rbTorsionIndices, const vector<vector<double> >& rbTorsionParameters,
             const vector<vector<int> >& bonded14Indices, double lj14Scale, double coulomb14Scale,
-            const vector<set<int> >& exclusions, const vector<vector<double> >& nonbondedParameters) {
+            const vector<set<int> >& exclusions, const vector<vector<double> >& nonbondedParameters,
+            NonbondedMethod nonbondedMethod, double nonbondedCutoff, double periodicBoxSize[3]) {
         verifyExclusions(exclusions);
         verify14(bonded14Indices);
     }