From 345b4180925667ef47a373a277d53cb2aeaa4fe0 Mon Sep 17 00:00:00 2001 From: Nick Papior Date: Tue, 12 Mar 2024 20:15:30 +0100 Subject: [PATCH] cleaned finder and removed pbc (getting it from Lattice) The neighbor finder got an overhaul to check for PBC and friends. Signed-off-by: Nick Papior --- src/sisl/geom/_neighbors/_finder.py | 41 ++++-------- src/sisl/geom/_neighbors/tests/test_finder.py | 62 +++++++++++++------ 2 files changed, 55 insertions(+), 48 deletions(-) diff --git a/src/sisl/geom/_neighbors/_finder.py b/src/sisl/geom/_neighbors/_finder.py index 15791687b1..2386621f43 100644 --- a/src/sisl/geom/_neighbors/_finder.py +++ b/src/sisl/geom/_neighbors/_finder.py @@ -362,16 +362,16 @@ def _scalar_to_cartesian_index(self, index): third, index = np.divmod(index, self.nbins[0] * self.nbins[1]) second, first = np.divmod(index, self.nbins[0]) - return np.array([first, second, third]).T + return np.column_stack([first, second, third]) def _correct_pairs_R_too_big( self, neighbor_pairs: np.ndarray, # (n_pairs, 5) split_ind: Union[int, np.ndarray], # (n_queried_atoms, ) - pbc: np.ndarray, # (3, ) ): """Correction to atom and supercell indices when the binning has been done on a tiled geometry""" is_sc_neigh = neighbor_pairs[:, 1] >= self.geometry.na + pbc = self.geometry.lattice.pbc invalid = None if not np.any(pbc): @@ -406,7 +406,6 @@ def find_neighbors( atoms: AtomsArgument = None, as_pairs: bool = False, self_interaction: bool = False, - pbc: Union[bool, Tuple[bool, bool, bool]] = (True, True, True), ): """Find neighbors as specified in the finder. @@ -426,9 +425,6 @@ def find_neighbors( the neighbors for each atom is returned. self_interaction: bool, optional whether to consider an atom a neighbor of itself. - pbc: bool or array-like of shape (3, ) - whether periodic conditions should be considered. - If a single bool is passed, all directions use that value. Returns ---------- @@ -443,9 +439,8 @@ def find_neighbors( # Sanitize atoms atoms = self.geometry._sanitize_atoms(atoms) - # Cast R and pbc into arrays of appropiate shape and type. + # Cast R into array of appropiate shape and type. thresholds = np.full(self._bins_geometry.na, self._aux_R, dtype=np.float64) - pbc = np.full(3, pbc, dtype=bool) # Get search indices search_indices, isc = self._get_search_indices( @@ -475,7 +470,7 @@ def find_neighbors( self_interaction, self._bins_geometry.xyz, self._bins_geometry.cell, - pbc, + self.geometry.lattice.pbc, thresholds, self._overlap, init_pairs, @@ -486,7 +481,7 @@ def find_neighbors( # we needed to create an auxiliary supercell. if self._R_too_big: neighbor_pairs, split_ind = self._correct_pairs_R_too_big( - neighbor_pairs, split_ind, pbc + neighbor_pairs, split_ind ) if as_pairs: @@ -496,12 +491,11 @@ def find_neighbors( # Split to get the neighbors of each atom return np.split(neighbor_pairs[:, 1:], split_ind, axis=0)[:-1] - def find_all_unique_pairs( + def find_unique_pairs( self, self_interaction: bool = False, - pbc: Union[bool, Tuple[bool, bool, bool]] = (True, True, True), ): - """Find all unique neighbor pairs within the geometry. + """Find unique neighbor pairs within the geometry. A pair of atoms is considered unique if atoms have the same index and correspond to the same unit cell. As an example, the connection @@ -525,9 +519,6 @@ def find_all_unique_pairs( the neighbors for each atom is returned. self_interaction: bool, optional whether to consider an atom a neighbor of itself. - pbc: bool or array-like of shape (3, ) - whether periodic conditions should be considered. - If a single bool is passed, all directions use that value. Returns ---------- @@ -546,7 +537,7 @@ def find_all_unique_pairs( if self._R_too_big: # Find all neighbors all_neighbors = self.find_neighbors( - as_pairs=True, self_interaction=self_interaction, pbc=pbc + as_pairs=True, self_interaction=self_interaction ) # Find out which of the pairs are uc connections @@ -560,9 +551,8 @@ def find_all_unique_pairs( # Concatenate the uc connections with the rest of the connections. return np.concatenate((uc_neighbors, all_neighbors[~is_uc_neigh])) - # Cast R and pbc into arrays of appropiate shape and type. + # Cast R into array of appropiate shape and type. thresholds = np.full(self.geometry.na, self.R, dtype=np.float64) - pbc = np.full(3, pbc, dtype=bool) # Get search indices search_indices, isc = self._get_search_indices( @@ -591,7 +581,7 @@ def find_all_unique_pairs( self_interaction, self.geometry.xyz, self.geometry.cell, - pbc, + self.geometry.lattice.pbc, thresholds, self._overlap, init_pairs, @@ -604,7 +594,6 @@ def find_close( self, xyz: Sequence, as_pairs: bool = False, - pbc: Union[bool, Tuple[bool, bool, bool]] = (True, True, True), ): """Find all atoms that are close to some coordinates in space. @@ -620,9 +609,6 @@ def find_close( If `True` pairs are returned, where the first item is the index of the point and the second one is the atom index. Otherwise a list containing the neighbors for each point is returned. - pbc: bool or array-like of shape (3, ) - whether periodic conditions should be considered. - If a single bool is passed, all directions use that value. Returns ---------- @@ -634,9 +620,8 @@ def find_close( If `as_pairs=False`: A list containing a numpy array of shape (n_neighs, 4) for each atom. """ - # Cast R and pbc into arrays of appropiate shape and type. + # Cast R into array of appropiate shape and type. thresholds = np.full(self._bins_geometry.na, self._aux_R, dtype=np.float64) - pbc = np.full(3, pbc, dtype=bool) xyz = np.atleast_2d(xyz) # Get search indices @@ -664,7 +649,7 @@ def find_close( self._list, self._bins_geometry.xyz, self._bins_geometry.cell, - pbc, + self.geometry.lattice.pbc, thresholds, init_pairs, self.memory[1], @@ -674,7 +659,7 @@ def find_close( # we needed to create an auxiliary supercell. if self._R_too_big: neighbor_pairs, split_ind = self._correct_pairs_R_too_big( - neighbor_pairs, split_ind, pbc + neighbor_pairs, split_ind ) if as_pairs: diff --git a/src/sisl/geom/_neighbors/tests/test_finder.py b/src/sisl/geom/_neighbors/tests/test_finder.py index 1c244265e1..407eea5d70 100644 --- a/src/sisl/geom/_neighbors/tests/test_finder.py +++ b/src/sisl/geom/_neighbors/tests/test_finder.py @@ -3,10 +3,10 @@ import numpy as np import pytest -from sisl import Geometry +from sisl import Geometry, Lattice from sisl.geom import NeighborFinder -pytestmark = [pytest.mark.neighbor] +pytestmark = [pytest.mark.geometry, pytest.mark.geom, pytest.mark.neighbor] tr_fixture = partial(pytest.fixture, scope="module", params=[True, False]) @@ -23,9 +23,17 @@ def request_param(request): pbc = tr_fixture()(request_param) +def set_pbc(geom, _pbc): + if _pbc: + geom.lattice.set_boundary_condition(Lattice.BC.PERIODIC) + else: + geom.lattice.set_boundary_condition(Lattice.BC.UNKNOWN) + + @pytest.fixture(scope="module") -def neighfinder(sphere_overlap, multiR): +def neighfinder(sphere_overlap, multiR, pbc): geom = Geometry([[0, 0, 0], [1.2, 0, 0], [9, 0, 0]], lattice=[10, 10, 7]) + set_pbc(geom, pbc) R = np.array([1.1, 1.5, 1.2]) if multiR else 1.5 @@ -114,9 +122,9 @@ def test_neighfinder_setup(sphere_overlap, multiR, post_setup): assert finder._counts.sum() == 2 -def test_neighbor_pairs(neighfinder, self_interaction, pbc, expected_neighs): +def test_neighbor_pairs(neighfinder, self_interaction, expected_neighs): neighs = neighfinder.find_neighbors( - as_pairs=True, self_interaction=self_interaction, pbc=pbc + as_pairs=True, self_interaction=self_interaction ) assert isinstance(neighs, np.ndarray) @@ -130,9 +138,9 @@ def test_neighbor_pairs(neighfinder, self_interaction, pbc, expected_neighs): assert np.all(neighs == [*first_at_neighs, *second_at_neighs, *third_at_neighs]) -def test_neighbors_lists(neighfinder, self_interaction, pbc, expected_neighs): +def test_neighbors_lists(neighfinder, self_interaction, expected_neighs): neighs = neighfinder.find_neighbors( - as_pairs=False, self_interaction=self_interaction, pbc=pbc + as_pairs=False, self_interaction=self_interaction ) assert isinstance(neighs, list) @@ -163,17 +171,13 @@ def test_neighbors_lists(neighfinder, self_interaction, pbc, expected_neighs): ), f"Wrong values for neighbors of atom {i}" -def test_all_unique_pairs(neighfinder, self_interaction, pbc, expected_neighs): +def test_all_unique_pairs(neighfinder, self_interaction, expected_neighs): if neighfinder.R.ndim == 1 and not neighfinder._overlap: with pytest.raises(ValueError): - neighfinder.find_all_unique_pairs( - self_interaction=self_interaction, pbc=pbc - ) + neighfinder.find_unique_pairs(self_interaction=self_interaction) return - neighs = neighfinder.find_all_unique_pairs( - self_interaction=self_interaction, pbc=pbc - ) + neighs = neighfinder.find_unique_pairs(self_interaction=self_interaction) first_at_neighs, second_at_neighs, third_at_neighs = expected_neighs @@ -196,7 +200,7 @@ def test_all_unique_pairs(neighfinder, self_interaction, pbc, expected_neighs): def test_close(neighfinder, pbc): - neighs = neighfinder.find_close([0.3, 0, 0], as_pairs=True, pbc=pbc) + neighs = neighfinder.find_close([0.3, 0, 0], as_pairs=True) expected_neighs = [[0, 1, 0, 0, 0], [0, 0, 0, 0, 0]] if pbc and neighfinder.R.ndim == 0: @@ -210,15 +214,16 @@ def test_no_neighbors(pbc): """Test the case where there are no neighbors, to see that it doesn't crash.""" geom = Geometry([[0, 0, 0]]) + set_pbc(geom, pbc) finder = NeighborFinder(geom, R=1.5) - neighs = finder.find_neighbors(as_pairs=True, pbc=pbc) + neighs = finder.find_neighbors(as_pairs=True) assert isinstance(neighs, np.ndarray) assert neighs.shape == (0, 5) - neighs = finder.find_neighbors(as_pairs=False, pbc=pbc) + neighs = finder.find_neighbors(as_pairs=False) assert isinstance(neighs, list) assert len(neighs) == 1 @@ -226,7 +231,7 @@ def test_no_neighbors(pbc): assert isinstance(neighs[0], np.ndarray) assert neighs[0].shape == (0, 4) - neighs = finder.find_all_unique_pairs(pbc=pbc) + neighs = finder.find_unique_pairs() assert isinstance(neighs, np.ndarray) assert neighs.shape == (0, 5) @@ -237,10 +242,11 @@ def test_R_too_big(pbc): than the unit cell.""" geom = Geometry([[0, 0, 0], [1, 0, 0]], lattice=[2, 10, 10]) + set_pbc(geom, pbc) neighfinder = NeighborFinder(geom, R=1.5) - neighs = neighfinder.find_all_unique_pairs(pbc=pbc) + neighs = neighfinder.find_unique_pairs() expected_neighs = [[0, 1, 0, 0, 0]] if pbc: @@ -252,7 +258,7 @@ def test_R_too_big(pbc): neighfinder = NeighborFinder(geom, R=[0.6, 2.2], overlap=True) - neighs = neighfinder.find_close([[0.5, 0, 0]], as_pairs=True, pbc=pbc) + neighs = neighfinder.find_close([[0.5, 0, 0]], as_pairs=True) expected_neighs = [[0, 1, 0, 0, 0], [0, 0, 0, 0, 0]] if pbc: @@ -283,3 +289,19 @@ def test_bin_sizes(): assert n3.nbins[0] == n4.nbins[0] assert n3.nbins[1] > n4.nbins[1] assert n3.nbins[2] > n4.nbins[2] + + +@pytest.mark.xfail(reason="Neighborfinder should locate the correct position") +def test_outside_box(): + # This should be at [0, 0, 0] and [1, 0, 0] + geom = Geometry([[0, 0, 0], [3, 0, 0]], lattice=[2, 10, 10]) + + # We should have fewer bins along the first lattice vector + n = NeighborFinder(geom, R=1.1) + + # Figure out if the neighbor of the first atom, will locate the + # other one, and in which isc + neighs = finder.find_neighbors() + + assert neighs[0] == [1, 0, 0, 0] + assert neighs[1] == [0, 0, 0, 0]