Fixed double-counting of periopdic image pairs

rascaline/src/systems/neighbors.rs

@@ -174,7 +174,7 @@ impl CellList {
         let n_cells = self.cells.shape();
         let n_cells = [n_cells[0], n_cells[1], n_cells[2]];
 
-        // find the cell in which this atom should go
+        // find the subcell in which this atom 'should go'
         let cell_index = [
             f64::floor(fractional[0] * n_cells[0] as f64) as i32,
             f64::floor(fractional[1] * n_cells[1] as f64) as i32,
@@ -250,10 +250,28 @@ impl CellList {
                                 let shift = CellShift(cell_shift) + atom_i.shift - atom_j.shift;
                                 let shift_is_zero = shift[0] == 0 && shift[1] == 0 && shift[2] == 0;
 
-                                if atom_i.index == atom_j.index && shift_is_zero {
-                                    // only create pair with the same atom twice
-                                    // if the pair spans more than one unit cell
-                                    continue;
+                                if atom_i.index == atom_j.index {
+                                    if shift_is_zero {
+                                        // only create pair with the same atom twice
+                                        // if the pair spans more than one unit cell
+                                        continue;
+                                    }
+
+                                    if (shift[0] + shift[1] + shift[2] < 0) ||
+                                        (
+                                            (shift[0] + shift[1] + shift[2] == 0) &&
+                                            (shift[2] < 0 || (shift[2] == 0 && shift[1] < 0))
+                                        )
+                                    {
+                                        // When creating pairs between an atom
+                                        // and one of its periodic images, the
+                                        // code generate multiple redundant
+                                        // pairs (e.g. with shifts 0 1 1 and 0
+                                        // -1 -1); and we want to only keep one
+                                        // of these. We keep the pair in the
+                                        // positive half plane of shifts.
+                                        continue;
+                                    }
                                 }
 
                                 if self.unit_cell.is_infinite() && !shift_is_zero {
@@ -426,21 +444,15 @@ mod tests {
         let neighbors = NeighborsList::new(&positions, cell, 3.0);
 
         let expected = [
-            (Vector3D::new(0.0, -1.0, -1.0), [-1, 0, 0]),
             (Vector3D::new(1.0, 0.0, -1.0),  [-1, 0, 1]),
             (Vector3D::new(1.0, -1.0, 0.0),  [-1, 1, 0]),
-            (Vector3D::new(-1.0, 0.0, -1.0), [0, -1, 0]),
             (Vector3D::new(0.0, 1.0, -1.0),  [0, -1, 1]),
-            (Vector3D::new(-1.0, -1.0, 0.0), [0, 0, -1]),
             (Vector3D::new(1.0, 1.0, 0.0),   [0, 0, 1]),
-            (Vector3D::new(0.0, -1.0, 1.0),  [0, 1, -1]),
             (Vector3D::new(1.0, 0.0, 1.0),   [0, 1, 0]),
-            (Vector3D::new(-1.0, 1.0, 0.0),  [1, -1, 0]),
-            (Vector3D::new(-1.0, 0.0, 1.0),  [1, 0, -1]),
             (Vector3D::new(0.0, 1.0, 1.0),   [1, 0, 0]),
         ];
 
-        assert_eq!(neighbors.pairs.len(), 12);
+        assert_eq!(neighbors.pairs.len(), 6);
         for (pair, (vector, shifts)) in neighbors.pairs.iter().zip(&expected) {
             assert_eq!(pair.first, 0);
             assert_eq!(pair.second, 0);
@@ -480,4 +492,26 @@ mod tests {
             assert_ulps_eq!(pair.distance, 2.0);
         }
     }
+
+    #[test]
+    fn small_cell_large_cutoffs() {
+        let cell = UnitCell::cubic(0.5);
+        let positions = [Vector3D::new(0.0, 0.0, 0.0)];
+        let neighbors = NeighborsList::new(&positions, cell, 0.6);
+
+        let expected = [
+            (Vector3D::new(0.0, 0.0, 0.5), [0, 0, 1]),
+            (Vector3D::new(0.0, 0.5, 0.0),  [0, 1, 0]),
+            (Vector3D::new(0.5, 0.0, 0.0),  [1, 0, 0]),
+        ];
+
+        assert_eq!(neighbors.pairs.len(), 3);
+        for (pair, (vector, shifts)) in neighbors.pairs.iter().zip(&expected) {
+            assert_eq!(pair.first, 0);
+            assert_eq!(pair.second, 0);
+            assert_ulps_eq!(pair.distance, 0.5);
+            assert_ulps_eq!(pair.vector, vector);
+            assert_eq!(&pair.cell_shift_indices, shifts);
+        }
+    }
 }