ReferenceNeighborList correctly handles triclinic boxes

platforms/reference/src/SimTKReference/ReferenceNeighborList.cpp

@@ -147,19 +147,33 @@ public:
         int dIndexY = int(maxDistance / voxelSizeY) + 1;
         int dIndexZ = int(maxDistance / voxelSizeZ) + 1;
         VoxelIndex centerVoxelIndex = getVoxelIndex(locationI);
-        int lastx = centerVoxelIndex.x+dIndexX;
-        int lasty = centerVoxelIndex.y+dIndexY;
-        int lastz = centerVoxelIndex.z+dIndexZ;
+        int minz = centerVoxelIndex.z-dIndexZ;
+        int maxz = centerVoxelIndex.z+dIndexZ;
+        if (usePeriodic)
+            maxz = min(maxz, minz+nz-1);
+        for (int z = minz; z <= maxz; ++z)
+        {
+            int boxz = (int) floor((float) z/nz);
+            int miny = centerVoxelIndex.y-dIndexY;
+            int maxy = centerVoxelIndex.y+dIndexY;
             if (usePeriodic) {
-            lastx = min(lastx, centerVoxelIndex.x-dIndexX+nx-1);
-            lasty = min(lasty, centerVoxelIndex.y-dIndexY+ny-1);
-            lastz = min(lastz, centerVoxelIndex.z-dIndexZ+nz-1);
+                int yoffset = (int) ceil(boxz*periodicBoxVectors[2][1]/voxelSizeY);
+                miny -= yoffset;
+                maxy -= yoffset-1;
+                maxy = min(maxy, miny+ny-1);
             }
-        for (int x = centerVoxelIndex.x - dIndexX; x <= lastx; ++x)
-        {
-            for (int y = centerVoxelIndex.y - dIndexY; y <= lasty; ++y)
+            for (int y = miny; y <= maxy; ++y)
             {
-                for (int z = centerVoxelIndex.z - dIndexZ; z <= lastz; ++z)
+                int boxy = (int) floor((float) y/ny);
+                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;
+                    maxx = min(maxx, minx+nx-1);
+                }
+                for (int x = minx; x <= maxx; ++x)
                 {
                     VoxelIndex voxelIndex(x, y, z);
                     if (usePeriodic) {

platforms/reference/tests/TestReferenceNeighborList.cpp

@@ -62,29 +62,24 @@ void testNeighborList()
     assert(neighborList.size() == 0);
 }
 
-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(RealVec& pos1, RealVec& pos2, const RealVec* periodicBoxVectors) {
+    RealVec diff = pos1-pos2;
+    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 diff.dot(diff);
 }
 
-double distance2(RealVec& pos1, RealVec& pos2, const RealVec& 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 numParticles, vector<RealVec>& positions, const RealVec& periodicBoxSize, double cutoff) {
+void verifyNeighborList(NeighborList& list, int numParticles, vector<RealVec>& positions, const RealVec* periodicBoxVectors, double cutoff) {
     for (int i = 0; i < (int) list.size(); i++) {
         int particle1 = list[i].first;
         int particle2 = list[i].second;
-        ASSERT(distance2(positions[particle1], positions[particle2], periodicBoxSize) <= cutoff*cutoff);
+        ASSERT(distance2(positions[particle1], positions[particle2], periodicBoxVectors) <= cutoff*cutoff);
     }
     int count = 0;
     for (int i = 0; i < numParticles; i++)
         for (int j = i+1; j < numParticles; j++)
-            if (distance2(positions[i], positions[j], periodicBoxSize) <= cutoff*cutoff)
+            if (distance2(positions[i], positions[j], periodicBoxVectors) <= cutoff*cutoff)
                 count++;
     ASSERT_EQUAL(count, list.size());
 }
@@ -92,7 +87,6 @@ void verifyNeighborList(NeighborList& list, int numParticles, vector<RealVec>& p
 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);
@@ -102,16 +96,40 @@ void testPeriodic() {
     init_gen_rand(0, sfmt);
 
     for (int i = 0; i <numParticles; i++) {
-        particleList[i][0] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxSize[0]*3);
-        particleList[i][1] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxSize[1]*3);
-        particleList[i][2] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxSize[2]*3);
+        particleList[i][0] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxVectors[0][0]*3);
+        particleList[i][1] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxVectors[1][1]*3);
+        particleList[i][2] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxVectors[2][2]*3);
+    }
+    vector<set<int> > exclusions(numParticles);
+    NeighborList neighborList;
+    computeNeighborListNaive(neighborList, numParticles, particleList, exclusions, periodicBoxVectors, true, cutoff);
+    verifyNeighborList(neighborList, numParticles, particleList, periodicBoxVectors, cutoff);
+    computeNeighborListVoxelHash(neighborList, numParticles, particleList, exclusions, periodicBoxVectors, true, cutoff);
+    verifyNeighborList(neighborList, numParticles, particleList, periodicBoxVectors, cutoff);
+}
+
+void testTriclinic() {
+    const int numParticles = 1000;
+    const double cutoff = 3.0;
+    RealVec periodicBoxVectors[3];
+    periodicBoxVectors[0] = RealVec(20, 0, 0);
+    periodicBoxVectors[1] = RealVec(5, 15, 0);
+    periodicBoxVectors[2] = RealVec(-3, -7, 22);
+    vector<RealVec> particleList(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+
+    for (int i = 0; i <numParticles; i++) {
+        particleList[i][0] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxVectors[0][0]*3);
+        particleList[i][1] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxVectors[1][1]*3);
+        particleList[i][2] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxVectors[2][2]*3);
     }
     vector<set<int> > exclusions(numParticles);
     NeighborList neighborList;
     computeNeighborListNaive(neighborList, numParticles, particleList, exclusions, periodicBoxVectors, true, cutoff);
-    verifyNeighborList(neighborList, numParticles, particleList, periodicBoxSize, cutoff);
+    verifyNeighborList(neighborList, numParticles, particleList, periodicBoxVectors, cutoff);
     computeNeighborListVoxelHash(neighborList, numParticles, particleList, exclusions, periodicBoxVectors, true, cutoff);
-    verifyNeighborList(neighborList, numParticles, particleList, periodicBoxSize, cutoff);
+    verifyNeighborList(neighborList, numParticles, particleList, periodicBoxVectors, cutoff);
 }
 
 int main() 
@@ -119,6 +137,7 @@ int main()
     try {
         testNeighborList();
         testPeriodic();
+        testTriclinic();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;