CpuNeighborList correctly handles triclinic boxes

platforms/cpu/include/CpuNeighborList.h

@@ -33,6 +33,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "AlignedArray.h"
+#include "RealVec.h"
 #include "windowsExportCpu.h"
 #include "openmm/internal/ThreadPool.h"
 #include <set>
@@ -46,6 +47,8 @@ public:
     class ThreadTask;
     class Voxels;
     CpuNeighborList(int blockSize);
+    void computeNeighborList(int numAtoms, const AlignedArray<float>& atomLocations, const std::vector<std::set<int> >& exclusions,
+            const RealVec* periodicBoxVectors, bool usePeriodic, float maxDistance, ThreadPool& threads);
     void computeNeighborList(int numAtoms, const AlignedArray<float>& atomLocations, const std::vector<std::set<int> >& exclusions,
             const float* periodicBoxSize, bool usePeriodic, float maxDistance, ThreadPool& threads);
     int getNumBlocks() const;
@@ -68,7 +71,7 @@ private:
     Voxels* voxels;
     const std::vector<std::set<int> >* exclusions;
     const float* atomLocations;
-    const float* periodicBoxSize;
+    RealVec periodicBoxVectors[3];
     int numAtoms;
     bool usePeriodic;
     float maxDistance;

platforms/cpu/src/CpuNonbondedForce.cpp

@@ -463,9 +463,9 @@ void CpuNonbondedForce::getDeltaR(const fvec4& posI, const fvec4& posJ, fvec4& d
     deltaR = posJ-posI;
     if (periodic) {
         if (triclinic) {
-            deltaR -= periodicBoxVec4[2]*floor(deltaR[2]*recipBoxSize[2]+0.5f);
-            deltaR -= periodicBoxVec4[1]*floor(deltaR[1]*recipBoxSize[1]+0.5f);
-            deltaR -= periodicBoxVec4[0]*floor(deltaR[0]*recipBoxSize[0]+0.5f);
+            deltaR -= periodicBoxVec4[2]*floorf(deltaR[2]*recipBoxSize[2]+0.5f);
+            deltaR -= periodicBoxVec4[1]*floorf(deltaR[1]*recipBoxSize[1]+0.5f);
+            deltaR -= periodicBoxVec4[0]*floorf(deltaR[0]*recipBoxSize[0]+0.5f);
         }
         else {
             fvec4 base = round(deltaR*invBoxSize)*boxSize;

platforms/cpu/tests/TestCpuNeighborList.cpp

@@ -48,10 +48,21 @@
 using namespace OpenMM;
 using namespace std;
 
-void testNeighborList(bool periodic) {
+void testNeighborList(bool periodic, bool triclinic) {
     const int numParticles = 500;
     const float cutoff = 2.0f;
-    const float boxSize[3] = {20.0f, 15.0f, 22.0f};
+    RealVec boxVectors[3];
+    if (triclinic) {
+        boxVectors[0] = RealVec(20, 0, 0);
+        boxVectors[1] = RealVec(5, 15, 0);
+        boxVectors[2] = RealVec(-3, -7, 22);
+    }
+    else {
+        boxVectors[0] = RealVec(20, 0, 0);
+        boxVectors[1] = RealVec(0, 15, 0);
+        boxVectors[2] = RealVec(0, 0, 22);
+    }
+    const float boxSize[3] = {boxVectors[0][0], boxVectors[1][1], boxVectors[2][2]};
     const int blockSize = 8;
     OpenMM_SFMT::SFMT sfmt;
     init_gen_rand(0, sfmt);
@@ -69,7 +80,7 @@ void testNeighborList(bool periodic) {
     }
     ThreadPool threads;
     CpuNeighborList neighborList(blockSize);
-    neighborList.computeNeighborList(numParticles, positions, exclusions, boxSize, periodic, cutoff, threads);
+    neighborList.computeNeighborList(numParticles, positions, exclusions, boxVectors, periodic, cutoff, threads);
     
     // Convert the neighbor list to a set for faster lookup.
     
@@ -94,15 +105,13 @@ void testNeighborList(bool periodic) {
     for (int i = 0; i < numParticles; i++)
         for (int j = 0; j <= i; j++) {
             bool shouldInclude = (exclusions[i].find(j) == exclusions[i].end());
-            float dx = positions[4*i]-positions[4*j];
-            float dy = positions[4*i+1]-positions[4*j+1];
-            float dz = positions[4*i+2]-positions[4*j+2];
+            Vec3 diff(positions[4*i]-positions[4*j], positions[4*i+1]-positions[4*j+1], positions[4*i+2]-positions[4*j+2]);
             if (periodic) {
-                dx -= floor(dx/boxSize[0]+0.5f)*boxSize[0];
-                dy -= floor(dy/boxSize[1]+0.5f)*boxSize[1];
-                dz -= floor(dz/boxSize[2]+0.5f)*boxSize[2];
+                diff -= boxVectors[2]*floor(diff[2]/boxSize[2]+0.5);
+                diff -= boxVectors[1]*floor(diff[1]/boxSize[1]+0.5);
+                diff -= boxVectors[0]*floor(diff[0]/boxSize[0]+0.5);
             }
-            if (dx*dx + dy*dy + dz*dz > cutoff*cutoff)
+            if (diff.dot(diff) > cutoff*cutoff)
                 shouldInclude = false;
             bool isIncluded = (neighbors.find(make_pair(i, j)) != neighbors.end() || neighbors.find(make_pair(j, i)) != neighbors.end());
             if (shouldInclude)
@@ -116,8 +125,9 @@ int main() {
             cout << "CPU is not supported.  Exiting." << endl;
             return 0;
         }
-        testNeighborList(false);
-        testNeighborList(true);
+        testNeighborList(false, false);
+        testNeighborList(true, false);
+        testNeighborList(true, true);
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/reference/src/SimTKReference/ReferenceNeighborList.cpp

@@ -157,9 +157,9 @@ public:
             int miny = centerVoxelIndex.y-dIndexY;
             int maxy = centerVoxelIndex.y+dIndexY;
             if (usePeriodic) {
-                int yoffset = (int) ceil(boxz*periodicBoxVectors[2][1]/voxelSizeY);
-                miny -= yoffset;
-                maxy -= yoffset-1;
+                double yoffset = boxz*periodicBoxVectors[2][1]/voxelSizeY;
+                miny -= (int) ceil(yoffset);
+                maxy -= (int) floor(yoffset);
                 maxy = min(maxy, miny+ny-1);
             }
             for (int y = miny; y <= maxy; ++y)
@@ -168,9 +168,9 @@ public:
                 int minx = centerVoxelIndex.x-dIndexX;
                 int maxx = centerVoxelIndex.x+dIndexX;
                 if (usePeriodic) {
-                    int xoffset = (int) ceil((boxy*periodicBoxVectors[1][0]+boxz*periodicBoxVectors[2][0])/voxelSizeX);
-                    minx -= xoffset;
-                    maxx -= xoffset-1;
+                    double xoffset = (boxy*periodicBoxVectors[1][0]+boxz*periodicBoxVectors[2][0])/voxelSizeX;
+                    minx -= (int) ceil(xoffset);
+                    maxx -= (int) floor(xoffset);
                     maxx = min(maxx, minx+nx-1);
                 }
                 for (int x = minx; x <= maxx; ++x)