From 8eae064355503ae3e69cfa204d07d60cf19bb671 Mon Sep 17 00:00:00 2001 From: Pol Febrer Date: Mon, 2 Mar 2026 18:58:12 +0100 Subject: [PATCH 1/2] Fix neighborlist for skewed cells --- src/sisl/geom/_neighbors/_finder.py | 59 ++++++++++++------- src/sisl/geom/_neighbors/tests/test_finder.py | 30 ++++++++++ 2 files changed, 69 insertions(+), 20 deletions(-) diff --git a/src/sisl/geom/_neighbors/_finder.py b/src/sisl/geom/_neighbors/_finder.py index 3099c99dc..cede80718 100644 --- a/src/sisl/geom/_neighbors/_finder.py +++ b/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: @@ -234,29 +235,47 @@ def setup( raise ValueError( "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 dicated by two lattice vectors (v1, v2) is: + # 2 * (max_R / sin[angle(v1, v2)] ) + min_bin_sizes = np.zeros(3) + lattice_norms = self.geometry.length + for i in range(3): + for j in range(3): + if i == j: continue + + # Compute angle of vectors i and j + scalar_prod = np.dot(self.geometry.cell[i], self.geometry.cell[j]) + alpha = np.acos( + scalar_prod / (lattice_norms[i] * lattice_norms[j]) + ) + + # Required bin size along vector i due to its relation with j + required_bin_size = 2 * max_R / np.sin(alpha) + min_bin_sizes[i] = max(min_bin_sizes[i], 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 + bin_size = min_bin_sizes * bin_size # 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 +291,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) diff --git a/src/sisl/geom/_neighbors/tests/test_finder.py b/src/sisl/geom/_neighbors/tests/test_finder.py index df197d1ef..4745aa94b 100644 --- a/src/sisl/geom/_neighbors/tests/test_finder.py +++ b/src/sisl/geom/_neighbors/tests/test_finder.py @@ -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 + From cfd3a14e1df6d875fd4feecda186d23bddf7f2de Mon Sep 17 00:00:00 2001 From: Nick Papior Date: Tue, 3 Mar 2026 14:10:19 +0100 Subject: [PATCH 2/2] minor cleanups and fixed regression on test --- changes/orphan.7.fix.rst | 8 ++++++ src/sisl/geom/_neighbors/_finder.py | 27 ++++++++++--------- src/sisl/geom/_neighbors/tests/test_finder.py | 20 +++++++------- 3 files changed, 32 insertions(+), 23 deletions(-) create mode 100644 changes/orphan.7.fix.rst diff --git a/changes/orphan.7.fix.rst b/changes/orphan.7.fix.rst new file mode 100644 index 000000000..012f7835d --- /dev/null +++ b/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!) diff --git a/src/sisl/geom/_neighbors/_finder.py b/src/sisl/geom/_neighbors/_finder.py index cede80718..41afa4c8f 100644 --- a/src/sisl/geom/_neighbors/_finder.py +++ b/src/sisl/geom/_neighbors/_finder.py @@ -189,7 +189,7 @@ def setup( 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. + 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. @@ -235,26 +235,26 @@ def setup( raise ValueError( "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 dicated by two lattice vectors (v1, v2) is: - # 2 * (max_R / sin[angle(v1, v2)] ) - min_bin_sizes = np.zeros(3) + # 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(3): - for j in range(3): - if i == j: continue + 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]) - alpha = np.acos( - scalar_prod / (lattice_norms[i] * lattice_norms[j]) - ) + cos_alpha = scalar_prod / (lattice_norms[i] * lattice_norms[j]) # Required bin size along vector i due to its relation with j - required_bin_size = 2 * max_R / np.sin(alpha) + # 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 < 1): @@ -263,7 +263,8 @@ def setup( f"is performed on neighboring bins. Received {bin_size}" ) - bin_size = min_bin_sizes * 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. diff --git a/src/sisl/geom/_neighbors/tests/test_finder.py b/src/sisl/geom/_neighbors/tests/test_finder.py index 4745aa94b..9f36579be 100644 --- a/src/sisl/geom/_neighbors/tests/test_finder.py +++ b/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] @@ -353,6 +353,7 @@ def test_outside_box(): 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 @@ -360,15 +361,15 @@ def test_skewed_cell(pbc): """ 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 + atom_1_xyz = atom_0_xyz + [1, 0.5, 0] # Atom 1.118 Ang away geom = Geometry( - [atom_0_xyz, atom_1_xyz], + [atom_0_xyz, atom_1_xyz], lattice=[ - [(R+0.01)*6, 0, 0], - [-3.5, (R+0.01)*6, 0], + [(R + 0.01) * 6, 0, 0], + [-3.5, (R + 0.01) * 6, 0], [0, 0, 20], - ] + ], ) set_pbc(geom, pbc) @@ -381,4 +382,3 @@ def test_skewed_cell(pbc): assert neighfinder.nbins == (2, 3, 8) # Atoms should have one neighbor assert neighfinder.find_neighbors()[0].n_neighbors == 1 -