Began implementing triclinic boxes in reference platform

platforms/reference/include/ReferenceForce.h

@@ -24,6 +24,7 @@
 #ifndef __ReferenceForce_H__
 #define __ReferenceForce_H__
 
+#include "RealVec.h"
 #include "lepton/CompiledExpression.h"
 #include "openmm/internal/windowsExport.h"
 
@@ -97,18 +98,20 @@ class OPENMM_EXPORT  ReferenceForce {
       static void getDeltaRPeriodic( const OpenMM::RealVec& atomCoordinatesI, const OpenMM::RealVec& atomCoordinatesJ,
                                              const RealOpenMM* boxSize, RealOpenMM* deltaR );
 
-    /**---------------------------------------------------------------------------------------
-      
-         Get deltaR between atomI and atomJ (static method): deltaR: j - i
-      
+      /**---------------------------------------------------------------------------------------
+
+         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 deltaR              deltaX, deltaY, deltaZ upon return
-      
+         @param boxVectors          the vectors defining the periodic box
+         @param deltaR              deltaX, deltaY, deltaZ, R2, R upon return
+
          --------------------------------------------------------------------------------------- */
-      
-      static void getDeltaROnly( const RealOpenMM* atomCoordinatesI, const RealOpenMM* atomCoordinatesJ,
-                                RealOpenMM* deltaR );
+
+      static void getDeltaRPeriodic( const OpenMM::RealVec& atomCoordinatesI, const OpenMM::RealVec& atomCoordinatesJ,
+                                             const OpenMM::RealVec* boxVectors, RealOpenMM* deltaR );
 
       /**
        * Get a pointer to the memory for setting a variable in a CompiledExpression.  If the expression

platforms/reference/include/ReferenceLJCoulombIxn.h

@@ -40,7 +40,7 @@ class ReferenceLJCoulombIxn {
       bool ewald;
       bool pme;
       const OpenMM::NeighborList* neighborList;
-      RealOpenMM periodicBoxSize[3];
+      OpenMM::RealVec periodicBoxVectors[3];
       RealOpenMM cutoffDistance, switchingDistance;
       RealOpenMM krf, crf;
       RealOpenMM alphaEwald;
@@ -118,11 +118,11 @@ class ReferenceLJCoulombIxn {
          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
+         @param vectors    the vectors defining the periodic box
       
          --------------------------------------------------------------------------------------- */
       
-      void setPeriodic( OpenMM::RealVec& boxSize );
+      void setPeriodic(OpenMM::RealVec* vectors);
        
       /**---------------------------------------------------------------------------------------
       

platforms/reference/include/ReferenceNeighborList.h

@@ -22,7 +22,7 @@ void OPENMM_EXPORT computeNeighborListNaive(
                               int nAtoms,
                               const AtomLocationList& atomLocations, 
                               const std::vector<std::set<int> >& exclusions,
-                              const RealVec& periodicBoxSize,
+                              const RealVec* periodicBoxVectors,
                               bool usePeriodic,
                               double maxDistance,
                               double minDistance = 0.0,
@@ -35,9 +35,9 @@ void OPENMM_EXPORT computeNeighborListNaive(
 void OPENMM_EXPORT computeNeighborListVoxelHash(
                               NeighborList& neighborList,
                               int nAtoms,
-                              const AtomLocationList& atomLocations, 
+                              const AtomLocationList& atomLocations,
                               const std::vector<std::set<int> >& exclusions,
-                              const RealVec& periodicBoxSize,
+                              const RealVec* periodicBoxVectors,
                               bool usePeriodic,
                               double maxDistance,
                               double minDistance = 0.0,

platforms/reference/include/ReferencePlatform.h

@@ -66,6 +66,7 @@ public:
     void* velocities;
     void* forces;
     void* periodicBoxSize;
+    void* periodicBoxVectors;
     void* constraints;
 };
 } // namespace OpenMM

platforms/reference/src/ReferenceKernels.cpp

@@ -134,6 +134,11 @@ static RealVec& extractBoxSize(ContextImpl& context) {
     return *(RealVec*) data->periodicBoxSize;
 }
 
+static RealVec* extractBoxVectors(ContextImpl& context) {
+    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+    return (RealVec*) data->periodicBoxVectors;
+}
+
 static ReferenceConstraints& extractConstraints(ContextImpl& context) {
     ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
     return *(ReferenceConstraints*) data->constraints;
@@ -256,10 +261,10 @@ void ReferenceUpdateStateDataKernel::getForces(ContextImpl& context, std::vector
 }
 
 void ReferenceUpdateStateDataKernel::getPeriodicBoxVectors(ContextImpl& context, Vec3& a, Vec3& b, Vec3& c) const {
-    RealVec& box = extractBoxSize(context);
-    a = Vec3(box[0], 0, 0);
-    b = Vec3(0, box[1], 0);
-    c = Vec3(0, 0, box[2]);
+    RealVec* vectors = extractBoxVectors(context);
+    a = vectors[0];
+    b = vectors[1];
+    c = vectors[2];
 }
 
 void ReferenceUpdateStateDataKernel::setPeriodicBoxVectors(ContextImpl& context, const Vec3& a, const Vec3& b, const Vec3& c) const {
@@ -267,6 +272,10 @@ void ReferenceUpdateStateDataKernel::setPeriodicBoxVectors(ContextImpl& context,
     box[0] = (RealOpenMM) a[0];
     box[1] = (RealOpenMM) b[1];
     box[2] = (RealOpenMM) c[2];
+    RealVec* vectors = extractBoxVectors(context);
+    vectors[0] = a;
+    vectors[1] = b;
+    vectors[2] = c;
 }
 
 void ReferenceUpdateStateDataKernel::createCheckpoint(ContextImpl& context, ostream& stream) {
@@ -854,15 +863,15 @@ double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool inc
     bool ewald  = (nonbondedMethod == Ewald);
     bool pme  = (nonbondedMethod == PME);
     if (nonbondedMethod != NoCutoff) {
-        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxSize(context), periodic || ewald || pme, nonbondedCutoff, 0.0);
+        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxVectors(context), periodic || ewald || pme, nonbondedCutoff, 0.0);
         clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
     }
     if (periodic || ewald || pme) {
-        RealVec& box = extractBoxSize(context);
+        RealVec* vectors = extractBoxVectors(context);
         double minAllowedSize = 1.999999*nonbondedCutoff;
-        if (box[0] < minAllowedSize || box[1] < minAllowedSize || box[2] < minAllowedSize)
+        if (vectors[0][0] < minAllowedSize || vectors[1][1] < minAllowedSize || vectors[2][2] < minAllowedSize)
             throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
-        clj.setPeriodic(box);
+        clj.setPeriodic(vectors);
     }
     if (ewald)
         clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
@@ -1018,7 +1027,7 @@ double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bo
     ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames);
     bool periodic = (nonbondedMethod == CutoffPeriodic);
     if (nonbondedMethod != NoCutoff) {
-        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxSize(context), periodic, nonbondedCutoff, 0.0);
+        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxVectors(context), periodic, nonbondedCutoff, 0.0);
         ixn.setUseCutoff(nonbondedCutoff, *neighborList);
     }
     if (periodic) {
@@ -1322,7 +1331,7 @@ double ReferenceCalcCustomGBForceKernel::execute(ContextImpl& context, bool incl
     if (periodic)
         ixn.setPeriodic(extractBoxSize(context));
     if (nonbondedMethod != NoCutoff) {
-        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxSize(context), periodic, nonbondedCutoff, 0.0);
+        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxVectors(context), periodic, nonbondedCutoff, 0.0);
         ixn.setUseCutoff(nonbondedCutoff, *neighborList);
     }
     map<string, double> globalParameters;

platforms/reference/src/ReferencePlatform.cpp

@@ -97,6 +97,7 @@ ReferencePlatform::PlatformData::PlatformData(const System& system) : time(0.0),
     velocities = new vector<RealVec>(numParticles);
     forces = new vector<RealVec>(numParticles);
     periodicBoxSize = new RealVec();
+    periodicBoxVectors = new RealVec[3];
     constraints = new ReferenceConstraints(system);
 }
 
@@ -105,5 +106,6 @@ ReferencePlatform::PlatformData::~PlatformData() {
     delete (vector<RealVec>*) velocities;
     delete (vector<RealVec>*) forces;
     delete (RealVec*) periodicBoxSize;
+    delete[] (RealVec*) periodicBoxVectors;
     delete (ReferenceConstraints*) constraints;
 }

platforms/reference/src/SimTKReference/ReferenceForce.cpp

@@ -41,13 +41,6 @@ using OpenMM::RealVec;
    --------------------------------------------------------------------------------------- */
 
 ReferenceForce::ReferenceForce( ){
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceForce::ReferenceForce";
-
-   // ---------------------------------------------------------------------------------------
-
 }
 
 /**---------------------------------------------------------------------------------------
@@ -57,13 +50,6 @@ ReferenceForce::ReferenceForce( ){
    --------------------------------------------------------------------------------------- */
 
 ReferenceForce::~ReferenceForce( ){
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceForce::~ReferenceForce";
-
-   // ---------------------------------------------------------------------------------------
-
 }
 
 /**---------------------------------------------------------------------------------------
@@ -79,27 +65,8 @@ RealOpenMM ReferenceForce::periodicDifference(RealOpenMM val1, RealOpenMM val2,
 
 }
 
-
-/**---------------------------------------------------------------------------------------
-
-   Get deltaR and distance and distance**2 between atomI and atomJ (static method)
-   deltaR: j - i
-
-   @param atomCoordinatesI    atom i coordinates
-   @param atomCoordinatesI    atom j coordinates
-   @param deltaR              deltaX, deltaY, deltaZ, R2, R upon return
-
-   --------------------------------------------------------------------------------------- */
-
 void ReferenceForce::getDeltaR( const RealVec& atomCoordinatesI, const RealVec& atomCoordinatesJ,
                                RealOpenMM* deltaR ){
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const std::string methodName = "\nReferenceForce::getDeltaR";
-
-   // ---------------------------------------------------------------------------------------
-
    deltaR[XIndex]    = atomCoordinatesJ[0] - atomCoordinatesI[0];
    deltaR[YIndex]    = atomCoordinatesJ[1] - atomCoordinatesI[1];
    deltaR[ZIndex]    = atomCoordinatesJ[2] - atomCoordinatesI[2];
@@ -108,27 +75,8 @@ void ReferenceForce::getDeltaR( const RealVec& atomCoordinatesI, const RealVec&
    deltaR[RIndex]    = (RealOpenMM) SQRT( deltaR[R2Index] );
 }
 
-/**---------------------------------------------------------------------------------------
-
-   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
-
-   --------------------------------------------------------------------------------------- */
-
 void ReferenceForce::getDeltaRPeriodic( const RealVec& atomCoordinatesI, const RealVec& 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]);
@@ -137,29 +85,17 @@ void ReferenceForce::getDeltaRPeriodic( const RealVec& atomCoordinatesI, const R
    deltaR[RIndex]    = (RealOpenMM) SQRT( deltaR[R2Index] );
 }
 
-/**---------------------------------------------------------------------------------------
-
-   Get deltaR between atomI and atomJ (static method); deltaR: j - i
-
-   @param atomCoordinatesI    atom i coordinates
-   @param atomCoordinatesI    atom j coordinates
-   @param deltaR              deltaX, deltaY, deltaZ upon return
-
-   --------------------------------------------------------------------------------------- */
-
-void ReferenceForce::getDeltaROnly( const RealOpenMM* atomCoordinatesI,
-                                   const RealOpenMM* atomCoordinatesJ,
-                                   RealOpenMM* deltaR ){
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const std::string methodName = "\nReferenceForce::getDeltaR";
-
-   // ---------------------------------------------------------------------------------------
-
-   deltaR[XIndex]    = atomCoordinatesJ[0] - atomCoordinatesI[0];
-   deltaR[YIndex]    = atomCoordinatesJ[1] - atomCoordinatesI[1];
-   deltaR[ZIndex]    = atomCoordinatesJ[2] - atomCoordinatesI[2];
+void ReferenceForce::getDeltaRPeriodic( const RealVec& atomCoordinatesI, const RealVec& atomCoordinatesJ,
+                               const RealVec* boxVectors, RealOpenMM* deltaR ){
+    RealVec diff = atomCoordinatesJ-atomCoordinatesI;
+    diff -= boxVectors[2]*floor(diff[2]/boxVectors[2][2]+0.5);
+    diff -= boxVectors[1]*floor(diff[1]/boxVectors[1][1]+0.5);
+    diff -= boxVectors[0]*floor(diff[0]/boxVectors[0][0]+0.5);
+    deltaR[XIndex] = diff[0];
+    deltaR[YIndex] = diff[1];
+    deltaR[ZIndex] = diff[2];
+    deltaR[R2Index] = diff.dot(diff);
+    deltaR[RIndex] = SQRT(deltaR[R2Index]);
 }
 
 double* ReferenceForce::getVariablePointer(Lepton::CompiledExpression& expression, const std::string& name) {

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp

@@ -112,20 +112,20 @@ void ReferenceLJCoulombIxn::setUseSwitchingFunction( RealOpenMM distance ) {
      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
+     @param vectors    the vectors defining the periodic box
 
      --------------------------------------------------------------------------------------- */
 
-  void ReferenceLJCoulombIxn::setPeriodic( RealVec& boxSize ) {
+  void ReferenceLJCoulombIxn::setPeriodic(OpenMM::RealVec* vectors) {
 
     assert(cutoff);
-    assert(boxSize[0] >= 2.0*cutoffDistance);
-    assert(boxSize[1] >= 2.0*cutoffDistance);
-    assert(boxSize[2] >= 2.0*cutoffDistance);
+    assert(vectors[0][0] >= 2.0*cutoffDistance);
+    assert(vectors[1][1] >= 2.0*cutoffDistance);
+    assert(vectors[2][2] >= 2.0*cutoffDistance);
     periodic = true;
-    periodicBoxSize[0] = boxSize[0];
-    periodicBoxSize[1] = boxSize[1];
-    periodicBoxSize[2] = boxSize[2];
+    periodicBoxVectors[0] = vectors[0];
+    periodicBoxVectors[1] = vectors[1];
+    periodicBoxVectors[2] = vectors[2];
   }
 
   /**---------------------------------------------------------------------------------------
@@ -197,7 +197,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
     RealOpenMM  factorEwald             = -1 / (4*alphaEwald*alphaEwald);
     RealOpenMM SQRT_PI                  = sqrt(PI_M);
     RealOpenMM TWO_PI                   = 2.0 * PI_M;
-    RealOpenMM recipCoeff               = (RealOpenMM)(ONE_4PI_EPS0*4*PI_M/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon);
+    RealOpenMM recipCoeff               = (RealOpenMM)(ONE_4PI_EPS0*4*PI_M/(periodicBoxVectors[0][0] * periodicBoxVectors[1][1] * periodicBoxVectors[2][2]) /epsilon);
 
     RealOpenMM totalSelfEwaldEnergy     = 0.0;
     RealOpenMM realSpaceEwaldEnergy     = 0.0;
@@ -237,6 +237,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
     vector<RealOpenMM> charges(numberOfAtoms);
     for (int i = 0; i < numberOfAtoms; i++)
         charges[i] = atomParameters[i][QIndex];
+    RealOpenMM periodicBoxSize[] = {periodicBoxVectors[0][0], periodicBoxVectors[1][1], periodicBoxVectors[2][2]};
     pme_exec(pmedata,atomCoordinates,forces,charges,periodicBoxSize,&recipEnergy,virial);
 
     if( totalEnergy )
@@ -255,7 +256,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
 
     // setup reciprocal box
 
-           RealOpenMM recipBoxSize[3] = { TWO_PI / periodicBoxSize[0], TWO_PI / periodicBoxSize[1], TWO_PI / periodicBoxSize[2]};
+         RealOpenMM recipBoxSize[3] = { TWO_PI / periodicBoxVectors[0][0], TWO_PI / periodicBoxVectors[1][1], TWO_PI / periodicBoxVectors[2][2]};
 
 
     // setup K-vectors
@@ -370,7 +371,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
        int jj = pair.second;
 
        RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
-       ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
+       ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR[0] );
        RealOpenMM r         = deltaR[0][ReferenceForce::RIndex];
        RealOpenMM inverseR  = one/(deltaR[0][ReferenceForce::RIndex]);
        RealOpenMM switchValue = 1, switchDeriv = 0;
@@ -556,7 +557,7 @@ void ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& at
     // get deltaR, R2, and R between 2 atoms
 
     if (periodic)
-        ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
+        ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR[0] );
     else
         ReferenceForce::getDeltaR( atomCoordinates[jj], atomCoordinates[ii], deltaR[0] );
 

platforms/reference/src/SimTKReference/ReferenceNeighborList.cpp

@@ -12,26 +12,15 @@ 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
-static double compPairDistanceSquared(const RealVec& pos1, const RealVec& pos2, const RealVec& periodicBoxSize, bool usePeriodic) {
-    double dx, dy, dz;
-    if (!usePeriodic) {
-        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]);
+static double compPairDistanceSquared(const RealVec& pos1, const RealVec& pos2, const RealVec* periodicBoxVectors, bool usePeriodic) {
+    RealVec diff = pos2-pos1;
+    if (usePeriodic) {
+        diff -= periodicBoxVectors[2]*floor(diff[2]/periodicBoxVectors[2][2]+0.5);
+        diff -= periodicBoxVectors[1]*floor(diff[1]/periodicBoxVectors[1][1]+0.5);
+        diff -= periodicBoxVectors[0]*floor(diff[0]/periodicBoxVectors[0][0]+0.5);
     }
-    return dx*dx + dy*dy + dz*dz;
+    return diff.dot(diff);
 }
 
 // Ridiculous O(n^2) version of neighbor list
@@ -41,7 +30,7 @@ void OPENMM_EXPORT computeNeighborListNaive(
                               int nAtoms,
                               const AtomLocationList& atomLocations, 
                               const vector<set<int> >& exclusions,
-                              const RealVec& periodicBoxSize,
+                              const RealVec* periodicBoxVectors,
                               bool usePeriodic,
                               double maxDistance,
                               double minDistance,
@@ -55,7 +44,7 @@ void OPENMM_EXPORT computeNeighborListNaive(
     {
         for (AtomIndex atomJ = atomI + 1; atomJ < (AtomIndex) nAtoms; ++atomJ)
         {
-            double pairDistanceSquared = compPairDistanceSquared(atomLocations[atomI], atomLocations[atomJ], periodicBoxSize, usePeriodic);
+            double pairDistanceSquared = compPairDistanceSquared(atomLocations[atomI], atomLocations[atomJ], periodicBoxVectors, usePeriodic);
             if ( (pairDistanceSquared <= maxDistanceSquared)  && (pairDistanceSquared >= minDistanceSquared))
                 if (exclusions[atomI].find(atomJ) == exclusions[atomI].end())
                 {
@@ -95,15 +84,15 @@ typedef std::vector< VoxelItem > Voxel;
 class VoxelHash
 {
 public:
-    VoxelHash(double vsx, double vsy, double vsz, const RealVec& periodicBoxSize, bool usePeriodic) :
-            voxelSizeX(vsx), voxelSizeY(vsy), voxelSizeZ(vsz), periodicBoxSize(periodicBoxSize), usePeriodic(usePeriodic) {
+    VoxelHash(double vsx, double vsy, double vsz, const RealVec* periodicBoxVectors, bool usePeriodic) :
+            voxelSizeX(vsx), voxelSizeY(vsy), voxelSizeZ(vsz), periodicBoxVectors(periodicBoxVectors), usePeriodic(usePeriodic) {
         if (usePeriodic) {
-            nx = (int) floor(periodicBoxSize[0]/voxelSizeX+0.5);
-            ny = (int) floor(periodicBoxSize[1]/voxelSizeY+0.5);
-            nz = (int) floor(periodicBoxSize[2]/voxelSizeZ+0.5);
-            voxelSizeX = periodicBoxSize[0]/nx;
-            voxelSizeY = periodicBoxSize[1]/ny;
-            voxelSizeZ = periodicBoxSize[2]/nz;
+            nx = (int) floor(periodicBoxVectors[0][0]/voxelSizeX+0.5);
+            ny = (int) floor(periodicBoxVectors[1][1]/voxelSizeY+0.5);
+            nz = (int) floor(periodicBoxVectors[2][2]/voxelSizeZ+0.5);
+            voxelSizeX = periodicBoxVectors[0][0]/nx;
+            voxelSizeY = periodicBoxVectors[1][1]/ny;
+            voxelSizeZ = periodicBoxVectors[2][2]/nz;
         }
     }
 
@@ -118,20 +107,15 @@ public:
 
 
     VoxelIndex getVoxelIndex(const RealVec& location) const {
-        double xperiodic, yperiodic, zperiodic;
-        if (!usePeriodic) {
-            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]);
+        RealVec r = location;
+        if (usePeriodic) {
+            r -= periodicBoxVectors[2]*floor(r[2]/periodicBoxVectors[2][2]);
+            r -= periodicBoxVectors[1]*floor(r[1]/periodicBoxVectors[1][1]);
+            r -= periodicBoxVectors[0]*floor(r[0]/periodicBoxVectors[0][0]);
         }
-        int x = int(floor(xperiodic / voxelSizeX));
-        int y = int(floor(yperiodic / voxelSizeY));
-        int z = int(floor(zperiodic / voxelSizeZ));
+        int x = int(floor(r[0]/voxelSizeX));
+        int y = int(floor(r[1]/voxelSizeY));
+        int z = int(floor(r[2]/voxelSizeZ));
         
         return VoxelIndex(x, y, z);
     }
@@ -194,7 +178,7 @@ public:
                         // Ignore self hits
                         if (atomI == atomJ) continue;
                         
-                        double dSquared = compPairDistanceSquared(locationI, locationJ, periodicBoxSize, usePeriodic);
+                        double dSquared = compPairDistanceSquared(locationI, locationJ, periodicBoxVectors, usePeriodic);
                         if (dSquared > maxDistanceSquared) continue;
                         if (dSquared < minDistanceSquared) continue;
                         
@@ -213,7 +197,7 @@ public:
 private:
     double voxelSizeX, voxelSizeY, voxelSizeZ;
     int nx, ny, nz;
-    const RealVec& periodicBoxSize;
+    const RealVec* periodicBoxVectors;
     const bool usePeriodic;
     std::map<VoxelIndex, Voxel> voxelMap;
 };
@@ -223,9 +207,9 @@ private:
 void OPENMM_EXPORT computeNeighborListVoxelHash(
                               NeighborList& neighborList,
                               int nAtoms,
-                              const AtomLocationList& atomLocations, 
+                              const AtomLocationList& atomLocations,
                               const vector<set<int> >& exclusions,
-                              const RealVec& periodicBoxSize,
+                              const RealVec* periodicBoxVectors,
                               bool usePeriodic,
                               double maxDistance,
                               double minDistance,
@@ -238,11 +222,11 @@ void OPENMM_EXPORT computeNeighborListVoxelHash(
     if (!usePeriodic)
         edgeSizeX = edgeSizeY = edgeSizeZ = maxDistance; // TODO - adjust this as needed
     else {
-        edgeSizeX = 0.5*periodicBoxSize[0]/floor(periodicBoxSize[0]/maxDistance);
-        edgeSizeY = 0.5*periodicBoxSize[1]/floor(periodicBoxSize[1]/maxDistance);
-        edgeSizeZ = 0.5*periodicBoxSize[2]/floor(periodicBoxSize[2]/maxDistance);
+        edgeSizeX = 0.5*periodicBoxVectors[0][0]/floor(periodicBoxVectors[0][0]/maxDistance);
+        edgeSizeY = 0.5*periodicBoxVectors[1][1]/floor(periodicBoxVectors[1][1]/maxDistance);
+        edgeSizeZ = 0.5*periodicBoxVectors[2][2]/floor(periodicBoxVectors[2][2]/maxDistance);
     }
-    VoxelHash voxelHash(edgeSizeX, edgeSizeY, edgeSizeZ, periodicBoxSize, usePeriodic);
+    VoxelHash voxelHash(edgeSizeX, edgeSizeY, edgeSizeZ, periodicBoxVectors, usePeriodic);
     for (AtomIndex atomJ = 0; atomJ < (AtomIndex) nAtoms; ++atomJ) // use "j", because j > i for pairs
     {
         // 1) Find other atoms that are close to this one

platforms/reference/tests/TestReferenceNeighborList.cpp

@@ -48,17 +48,17 @@ void testNeighborList()
     
     NeighborList neighborList;
 
-    RealVec boxSize;
-    computeNeighborListNaive(neighborList, 2, particleList, exclusions, boxSize, false, 13.7, 0.01);
+    RealVec boxVectors[3];
+    computeNeighborListNaive(neighborList, 2, particleList, exclusions, boxVectors, false, 13.7, 0.01);
     assert(neighborList.size() == 1);
     
-    computeNeighborListNaive(neighborList, 2, particleList, exclusions, boxSize, false, 13.5, 0.01);
+    computeNeighborListNaive(neighborList, 2, particleList, exclusions, boxVectors, false, 13.5, 0.01);
     assert(neighborList.size() == 0);
     
-    computeNeighborListVoxelHash(neighborList, 2, particleList, exclusions, boxSize, false, 13.7, 0.01);
+    computeNeighborListVoxelHash(neighborList, 2, particleList, exclusions, boxVectors, false, 13.7, 0.01);
     assert(neighborList.size() == 1);
     
-    computeNeighborListVoxelHash(neighborList, 2, particleList, exclusions, boxSize, false, 13.5, 0.01);
+    computeNeighborListVoxelHash(neighborList, 2, particleList, exclusions, boxVectors, false, 13.5, 0.01);
     assert(neighborList.size() == 0);
 }
 
@@ -93,6 +93,10 @@ void testPeriodic() {
     const int numParticles = 100;
     const double cutoff = 3.0;
     const RealVec periodicBoxSize(20.0, 15.0, 22.0);
+    RealVec periodicBoxVectors[3];
+    periodicBoxVectors[0] = RealVec(20, 0, 0);
+    periodicBoxVectors[1] = RealVec(0, 15, 0);
+    periodicBoxVectors[2] = RealVec(0, 0, 22);
     vector<RealVec> particleList(numParticles);
     OpenMM_SFMT::SFMT sfmt;
     init_gen_rand(0, sfmt);
@@ -104,9 +108,9 @@ void testPeriodic() {
     }
     vector<set<int> > exclusions(numParticles);
     NeighborList neighborList;
-    computeNeighborListNaive(neighborList, numParticles, particleList, exclusions, periodicBoxSize, true, cutoff);
+    computeNeighborListNaive(neighborList, numParticles, particleList, exclusions, periodicBoxVectors, true, cutoff);
     verifyNeighborList(neighborList, numParticles, particleList, periodicBoxSize, cutoff);
-    computeNeighborListVoxelHash(neighborList, numParticles, particleList, exclusions, periodicBoxSize, true, cutoff);
+    computeNeighborListVoxelHash(neighborList, numParticles, particleList, exclusions, periodicBoxVectors, true, cutoff);
     verifyNeighborList(neighborList, numParticles, particleList, periodicBoxSize, cutoff);
 }
 

platforms/reference/tests/TestReferenceNonbondedForce.cpp

@@ -41,6 +41,7 @@
 #include "openmm/VerletIntegrator.h"
 #include "SimTKOpenMMRealType.h"
 #include "openmm/HarmonicBondForce.h"
+#include "sfmt/SFMT.h"
 #include <iostream>
 #include <vector>
 
@@ -351,6 +352,69 @@ void testPeriodic() {
     ASSERT_EQUAL_TOL(2*ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
 }
 
+void testTriclinic() {
+    ReferencePlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    Vec3 a(3.1, 0, 0);
+    Vec3 b(0.4, 3.5, 0);
+    Vec3 c(-0.1, -0.5, 4.0);
+    system.setDefaultPeriodicBoxVectors(a, b, c);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->addParticle(1.0, 1, 0);
+    nonbonded->addParticle(1.0, 1, 0);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    const double cutoff = 1.5;
+    nonbonded->setCutoffDistance(cutoff);
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    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);
+    for (int iteration = 0; iteration < 50; iteration++) {
+        // Generate random positions for the two particles.
+
+        positions[0] = a*genrand_real2(sfmt) + b*genrand_real2(sfmt) + c*genrand_real2(sfmt);
+        positions[1] = a*genrand_real2(sfmt) + b*genrand_real2(sfmt) + c*genrand_real2(sfmt);
+        context.setPositions(positions);
+
+        // Loop over all possible periodic copies and find the nearest one.
+
+        Vec3 delta;
+        double distance2 = 100.0;
+        for (int i = -1; i < 2; i++)
+            for (int j = -1; j < 2; j++)
+                for (int k = -1; k < 2; k++) {
+                    Vec3 d = positions[1]-positions[0]+a*i+b*j+c*k;
+                    if (d.dot(d) < distance2) {
+                        delta = d;
+                        distance2 = d.dot(d);
+                    }
+                }
+        double distance = sqrt(distance2);
+
+        // See if the force and energy are correct.
+
+        State state = context.getState(State::Forces | State::Energy);
+        if (distance >= cutoff) {
+            ASSERT_EQUAL(0.0, state.getPotentialEnergy());
+            ASSERT_EQUAL_VEC(Vec3(0, 0, 0), state.getForces()[0], 0);
+            ASSERT_EQUAL_VEC(Vec3(0, 0, 0), state.getForces()[1], 0);
+        }
+        else {
+            const Vec3 force = delta*ONE_4PI_EPS0*(-1.0/(distance*distance*distance)+2.0*krf);
+            ASSERT_EQUAL_TOL(ONE_4PI_EPS0*(1.0/distance+krf*distance*distance-crf), state.getPotentialEnergy(), TOL);
+            ASSERT_EQUAL_VEC(force, state.getForces()[0], TOL);
+            ASSERT_EQUAL_VEC(-force, state.getForces()[1], TOL);
+        }
+    }
+}
+
 void testDispersionCorrection() {
     // Create a box full of identical particles.
 
@@ -483,6 +547,7 @@ int main() {
         testCutoff();
         testCutoff14();
         testPeriodic();
+        testTriclinic();
         testDispersionCorrection();
         testSwitchingFunction(NonbondedForce::CutoffNonPeriodic);
         testSwitchingFunction(NonbondedForce::PME);

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

@@ -74,6 +74,11 @@ static RealVec& extractBoxSize(ContextImpl& context) {
     return *(RealVec*) data->periodicBoxSize;
 }
 
+static RealVec* extractBoxVectors(ContextImpl& context) {
+    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+    return (RealVec*) data->periodicBoxVectors;
+}
+
 // ***************************************************************************
 
 ReferenceCalcAmoebaBondForceKernel::ReferenceCalcAmoebaBondForceKernel(std::string name, const Platform& platform, const System& system) : 
@@ -968,7 +973,7 @@ double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool inc
     RealOpenMM energy;
     if( useCutoff ){
         vdwForce.setCutoff( cutoff );
-        computeNeighborListVoxelHash( *neighborList, numParticles, posData, allExclusions, extractBoxSize(context), usePBC, cutoff, 0.0);
+        computeNeighborListVoxelHash( *neighborList, numParticles, posData, allExclusions, extractBoxVectors(context), usePBC, cutoff, 0.0);
         if( usePBC ){
             vdwForce.setNonbondedMethod( AmoebaReferenceVdwForce::CutoffPeriodic);
             RealVec& box = extractBoxSize(context);