Fix neighborlist for skewed cells (#976)

* Fix neighborlist for skewed cells

* minor cleanups and fixed regression on test

---------

Co-authored-by: Nick Papior <nickpapior@gmail.com>

Author: pfebrer

changes/orphan.7.fix.rst

@@ -0,0 +1,8 @@
+Fixed skewed cell neighbor search
+
+For the linear neighbor search mechanism
+it needed to account for skewed cells
+when the binsize is very close to the inter-atomic
+distances.
+
+(reported by @ialcon on Discord, thanks!)

src/sisl/geom/_neighbors/_finder.py

@@ -144,7 +144,7 @@ def __init__(
         geometry: Geometry,
         R: Optional[Union[float, np.ndarray]] = None,
         overlap: bool = False,
-        bin_size: Union[float, tuple[float, float, float]] = 2,
+        bin_size: Union[float, tuple[float, float, float]] = 1,
     ):
         self.setup(geometry, R=R, overlap=overlap, bin_size=bin_size)
 
@@ -153,7 +153,7 @@ def setup(
         geometry: Optional[Geometry] = None,
         R: Optional[Union[float, np.ndarray]] = None,
         overlap: bool = None,
-        bin_size: Union[float, tuple[float, float, float]] = 2,
+        bin_size: Union[float, tuple[float, float, float]] = 1,
     ):
         r"""Prepares everything for neighbor finding.
 
@@ -184,14 +184,15 @@ def setup(
             is within the sphere of the first atom. Note that this implies that
             atom :math:`I` might be atom :math:`J`'s neighbor while the opposite is not true.
         bin_size :
-            the factor for the radius to determine how large the bins are,
-            optionally along each lattice vector.
-            It can minimally be 2, meaning that the maximum radius to consider
-            is twice the radius considered. For larger values, more atoms will be in
-            each bin (and thus fewer bins).
-            Hence, this value can be used to fine-tune the memory requirement by
-            decreasing number of bins, at the cost of a bit more run-time searching
-            bins.
+            the neighbor finding algorithm divides space with a uniform grid.
+            The minimal size of each grid bin depends on the radius for neighbor
+            search and the lattice vectors (the more skewed the cell, the bigger
+            the bins need to be).
+            This argument accepts a factor by which to multiply the bin size,
+            with 1 resulting in the minimum bin size.
+            Larger values reduce the amount of memory needed by the algorithm, but
+            will slow down the neighbor finding since there will be more potential
+            neighbors for each atom.
         """
         # Set the geometry. Copy it because we may need to modify the supercell size.
         if geometry is not None:
@@ -235,28 +236,47 @@ def setup(
                 "All R values are 0 or less. Please provide some positive values"
             )
 
+        # Find the minimum length needed in each lattice vector direction
+        # (the more skewed the cell is, the bigger the bins have to be).
+        # The minimum bin size dictated by two lattice vectors (v1, v2) is:
+        #   2 * (max_R / sin[angle(v1, v2)] )
+        # The factor 2 max_R is taken out and applied after.
+        min_bin_sizes = np.ones(3)
+        lattice_norms = self.geometry.length
+        for i in range(2):
+            for j in range(i + 1, 3):
+
+                # Compute angle of vectors i and j
+                scalar_prod = np.dot(self.geometry.cell[i], self.geometry.cell[j])
+                cos_alpha = scalar_prod / (lattice_norms[i] * lattice_norms[j])
+
+                # Required bin size along vector i due to its relation with j
+                # This ensures a positive number
+                required_bin_size = 1 / (1 - cos_alpha**2) ** 0.5
+                min_bin_sizes[i] = max(min_bin_sizes[i], required_bin_size)
+                min_bin_sizes[j] = max(min_bin_sizes[j], required_bin_size)
+
         bin_size = np.asarray(bin_size)
-        if np.any(bin_size < 2):
+        if np.any(bin_size < 1):
             raise ValueError(
-                "The bin_size must be larger than 2 to only search in the "
-                "neighboring bins. Please increase to a value >=2"
+                "The bin_size argument must be larger than 1 so that the search"
+                f"is performed on neighboring bins. Received {bin_size}"
             )
 
-        bin_size = max_R * bin_size
+        # Correct for maximum distance in the bin-size (R on both sides)
+        bin_size = min_bin_sizes * bin_size * 2 * max_R
 
         # We add a small amount to bin_size to avoid ambiguities when
         # a position is exactly at the center of a bin.
         bin_size += 0.001
 
-        lattice_sizes = self.geometry.length
-
-        self._R_too_big = np.any(bin_size > lattice_sizes)
+        self._R_too_big = np.any(bin_size > lattice_norms)
         if self._R_too_big:
             # This means that nsc must be at least 5.
 
             # We round the amount of cells needed in each direction
             # to the closest next odd number.
-            nsc = np.ceil(bin_size / lattice_sizes) // 2 * 2 + 1
+            nsc = np.ceil(bin_size / lattice_norms) // 2 * 2 + 1
             # And then set it as the number of supercells.
             self.geometry.set_nsc(nsc.astype(int))
             if self._aux_R.ndim == 1:
@@ -272,13 +292,13 @@ def setup(
             )
 
             # Recompute lattice sizes
-            lattice_sizes = self._bins_geometry.length
+            lattice_norms = self._bins_geometry.length
 
         else:
             self._bins_geometry = self.geometry
 
         # Get the number of bins along each cell direction.
-        nbins_float = lattice_sizes / bin_size
+        nbins_float = lattice_norms / bin_size
         self.nbins = tuple(np.floor(nbins_float).astype(int))
         self.total_nbins = np.prod(self.nbins)
 

src/sisl/geom/_neighbors/tests/test_finder.py

@@ -325,8 +325,8 @@ def test_bin_sizes():
     geom = Geometry([[0, 0, 0], [1, 0, 0]], lattice=[2, 10, 10])
 
     # We should have fewer bins along the first lattice vector
-    n1 = NeighborFinder(geom, R=1.5, bin_size=2)
-    n2 = NeighborFinder(geom, R=1.5, bin_size=4)
+    n1 = NeighborFinder(geom, R=1.5, bin_size=1)
+    n2 = NeighborFinder(geom, R=1.5, bin_size=2)
 
     assert n1.total_nbins > n2.total_nbins
     # When the bin is bigger than the unit cell, this situation
@@ -336,8 +336,8 @@ def test_bin_sizes():
     assert n1.nbins[2] > n2.nbins[2]
 
     # We should have the same number of bins the 2nd and 3rd lattice vectors
-    n3 = NeighborFinder(geom, R=1.5, bin_size=2)
-    n4 = NeighborFinder(geom, R=1.5, bin_size=(2, 4, 4))
+    n3 = NeighborFinder(geom, R=1.5, bin_size=1)
+    n4 = NeighborFinder(geom, R=1.5, bin_size=(1, 2, 2))
 
     assert n3.nbins[0] == n4.nbins[0]
     assert n3.nbins[1] > n4.nbins[1]
@@ -352,3 +352,33 @@ def test_outside_box():
     # IT SHOULD BE SUPPORTED IN THE FUTURE
     with pytest.raises(ValueError):
         n = NeighborFinder(geom, R=1.1)
+
+
+def test_skewed_cell(pbc):
+    """Skewed cells are more complex for determining bin size,
+    here we check that the finder works for a case where treating
+    the skewed cell as an orthogonal cell will not work.
+    """
+    R = 1.12
+    atom_0_xyz = np.array([R - 0.004, 0, 0])
+    atom_1_xyz = atom_0_xyz + [1, 0.5, 0]  # Atom 1.118 Ang away
+
+    geom = Geometry(
+        [atom_0_xyz, atom_1_xyz],
+        lattice=[
+            [(R + 0.01) * 6, 0, 0],
+            [-3.5, (R + 0.01) * 6, 0],
+            [0, 0, 20],
+        ],
+    )
+    set_pbc(geom, pbc)
+
+    neighfinder = NeighborFinder(geom, R=R, overlap=False)
+
+    neighfinder.assert_consistency()
+
+    # Check that the finder took into account the skewed cell
+    # by dividing space in (2, 3, 8) instead of (3, 3, 8).
+    assert neighfinder.nbins == (2, 3, 8)
+    # Atoms should have one neighbor
+    assert neighfinder.find_neighbors()[0].n_neighbors == 1