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 <[email protected]>

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:

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,23 +214,24 @@ 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
 
     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]