Properly sort the pair list when there are duplicated i-j pairs

python/vesin/tests/data/Cd2I4O12.xyz

@@ -0,0 +1,20 @@
+18
+Lattice="5.7779061565689 0.0 0.0 0.0 5.9353021159594 0.0 0.0 0.0 32.37048911" Properties=species:S:1:pos:R:3:masses:R:1:forces:R:3 pbc="T T T"
+Cd       1.66027661       1.01057390      14.51270543     112.41400000      -0.03865710      -0.05802875       0.02849028
+Cd       4.54922969       1.95707716      17.85778368     112.41400000      -0.03865710       0.05802875      -0.02849028
+I        5.63439709       3.76166207      14.85737491     126.90447000      -1.47061191       1.06604780       0.08108146
+I        2.74544401       5.14129111      17.51311420     126.90447000      -1.47061191      -1.06604780      -0.08108146
+I        4.31565215       1.38050855      12.64875342     126.90447000       0.74986725       0.79087438      -0.69424655
+I        1.42669907       1.58714251      19.72173569     126.90447000       0.74986725      -0.79087438       0.69424655
+O        2.78364007       0.39473043      12.50000000      15.99900000      -1.11916203      -0.61054277      -0.10327189
+O        5.67259315       2.57292062      19.87048911      15.99900000      -1.11916203       0.61054277       0.10327189
+O        3.64269369       2.38202798      14.05899015      15.99900000      -0.44450476       0.58457571       0.81988628
+O        0.75374061       0.58562308      18.31149896      15.99900000      -0.44450476      -0.58457571      -0.81988628
+O        5.44167386       0.15724415      13.43179857      15.99900000       0.93009003      -0.74236748       0.47707191
+O        2.55272078       2.81040691      18.93869054      15.99900000       0.93009003       0.74236748      -0.47707191
+O        1.34512782       4.79252145      15.30307812      15.99900000       0.80628161       0.77936335       0.23325304
+O        4.23408090       4.11043172      17.06741099      15.99900000       0.80628161      -0.77936335      -0.23325304
+O        0.52359753       2.92010411      13.34319181      15.99900000       0.37932513      -0.74376411      -0.97400744
+O        3.41255061       0.04754695      19.02729730      15.99900000       0.37932513       0.74376411       0.97400744
+O        0.27893312       2.33486569      15.96606058      15.99900000       0.20737178      -1.07627874       0.77025494
+O        3.16788619       0.63278537      16.40442853      15.99900000       0.20737178       1.07627874      -0.77025494

python/vesin/tests/test_neighbors.py

@@ -85,6 +85,15 @@ def test_sorting():
     # check that unsorted is not already sorted by chance
     assert not np.all(sorted_ij == unsorted_ij)
 
+    # https://github.com/Luthaf/vesin/issues/34
+    atoms = ase.io.read(f"{CURRENT_DIR}/data/Cd2I4O12.xyz")
+    calculator = vesin.NeighborList(cutoff=5.0, full_list=True, sorted=True)
+    i, j = calculator.compute(
+        points=atoms.positions, box=atoms.cell[:], periodic=True, quantities="ij"
+    )
+    sorted_ij = np.concatenate((i.reshape(-1, 1), j.reshape(-1, 1)), axis=1)
+    assert np.all(sorted_ij[np.lexsort((j, i))] == sorted_ij)
+
 
 def test_errors():
     points = np.array([[0.0, 0.0, 0.0], [0.5, 0.5, 0.5]])

python/vesin/vesin/_ase.py

@@ -45,8 +45,7 @@ def ase_neighbor_list(quantities, a, cutoff, self_interaction=False, max_nbins=0
         periodic = False
     else:
         raise ValueError(
-            "different periodic boundary conditions on different axis "
-            "are not supported"
+            "different periodic boundary conditions on different axis are not supported"
         )
 
     # sorted=True and full_list=True since that's what ASE does

vesin/src/cpu_cell_list.cpp

@@ -446,6 +446,7 @@ void GrowableNeighborList::sort() {
 
     // step 2: permute all data according to the sorted indices.
     int64_t cur = 0;
+    int64_t is_sorted_up_to = 0;
     // data in `from` should go to `cur`
     auto from = indices[cur];
 
@@ -509,8 +510,10 @@ void GrowableNeighborList::sort() {
         indices[cur] = -1;
 
         // look for the next loop of permutation
+        cur = is_sorted_up_to;
         while (indices[cur] == -1) {
             cur += 1;
+            is_sorted_up_to += 1;
             if (cur == this->length()) {
                 break;
             }