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 3099c99dc..41afa4c8f 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: @@ -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) diff --git a/src/sisl/geom/_neighbors/tests/test_finder.py b/src/sisl/geom/_neighbors/tests/test_finder.py index df197d1ef..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] @@ -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