diff --git a/sisl/geom/neighs/fneighs.f90 b/sisl/geom/neighs/fneighs.f90 index 281d2dff8..ce5b5b06d 100644 --- a/sisl/geom/neighs/fneighs.f90 +++ b/sisl/geom/neighs/fneighs.f90 @@ -27,6 +27,7 @@ subroutine build_table(nbins, bin_indices, Nat, list, heads, counts) ! first atom that is contained in it. ! counts: ! For each bin, the number of atoms that it contains. + implicit none integer, intent(in) :: nbins integer, intent(in) :: bin_indices(Nat) @@ -58,11 +59,10 @@ subroutine build_table(nbins, bin_indices, Nat, list, heads, counts) end subroutine build_table -subroutine get_pairs(at_indices, indices, heads, list, max_npairs, & - self_interaction, xyz, thresholds, sphere_overlap, N_ind, N, & +subroutine get_pairs(at_indices, indices, iscs, heads, list, max_npairs, & + self_interaction, xyz, cell, pbc, thresholds, sphere_overlap, N_ind, N, & neighs, split_indices) - ! Gets (possibly duplicated) pairs of atoms that can potentially - ! be neighbours. + ! Gets (possibly duplicated) pairs of neighbour atoms. ! ! INPUTS ! --------- @@ -84,16 +84,36 @@ subroutine get_pairs(at_indices, indices, heads, list, max_npairs, & ! Indices are 0-indexed (python indices). If an item is ! -1, that there are no more atoms in the same bin. ! This array is constructed by `build_table`. - ! npairs: - ! The number of pairs that are to be found in total. This - ! is computed in python with the help of the `counts` array - ! constructed by `build_table`. + ! max_npairs: + ! The number of maximum pairs that can be found. + ! It is used to allocate the `neighs` array. This is computed + ! in python with the help of the `counts` array constructed + ! by `build_table`. ! self_interaction: bool, optional ! whether to consider an atom a neighbour of itself. + ! xyz: + ! the cartesian coordinates of all atoms in the geometry. + ! cell: + ! the lattice vectors of the geometry, where cell(i, :) is the + ! ith lattice vector. + ! pbc: + ! for each direction, whether periodic boundary conditions should + ! be considered. + ! thresholds: + ! the threshold radius for each atom in the geometry. + ! sphere_overlap: + ! If true, two atoms are considered neighbours if their spheres + ! overlap. + ! If false, two atoms are considered neighbours if the second atom + ! is within the sphere of the first atom. Note that this implies that + ! atom `i` might be atom `j`'s neighbour while the opposite is not true. ! N_ind: ! The number of atoms for which we will look for neighbours. ! No need to pass it when calling from python, it will be ! inferred from the shape of `at_indices` and `indices`. + ! N: + ! total number of atoms. No need to pass it when calling from + ! python. It will be inferred from the shapes of the arrays. ! ! ! OUTPUTS @@ -106,20 +126,26 @@ subroutine get_pairs(at_indices, indices, heads, list, max_npairs, & ! breakpoints delimit the position where one can find pairs ! for each 8 bin search. It can be used later to quickly ! split `neighs` on the multiple individual searches. + implicit none + integer, intent(in) :: at_indices(N_ind) - integer, intent(in) :: heads(:), list(:), indices(N_ind, 8) + integer, intent(in) :: indices(N_ind, 8), iscs(N_ind, 8, 3) + integer, intent(in) :: heads(:), list(N) integer, intent(in) :: max_npairs logical, intent(in) :: self_interaction - real(8), intent(in) :: xyz(N, 3), thresholds(N) + real(8), intent(in) :: xyz(N, 3), cell(3,3), thresholds(N) + logical, intent(in) :: pbc(3) logical, intent(in) :: sphere_overlap integer, intent(in) :: N, N_ind - integer, intent(out) :: neighs(max_npairs, 2) + integer, intent(out) :: neighs(max_npairs, 5) integer, intent(out) :: split_indices(N_ind) ! Auxiliar variables integer :: search_index integer :: bin_index, at, neigh_at, i_pair, j - real(8) :: dist, threshold + real(8) :: dist, threshold, ref_xyz(3) + integer :: isc(3) + logical :: is_neigh_cell ! Initialize the pair index, which will tell us the location ! to store each pair that we find. Each time we store a new pair @@ -133,27 +159,61 @@ subroutine get_pairs(at_indices, indices, heads, list, max_npairs, & do j=1, 8 ! Find the bin index. Sum 1 to get the fortran index. bin_index = indices(search_index, j) + 1 + ! Get the first atom index in this bin neigh_at = heads(bin_index) + ! If there are no atoms in this bin, do not even bother + ! checking supercell indices. + if (neigh_at == -1) cycle + + ! Find the supercell indices for this bin + isc(:) = iscs(search_index, j, :) + ! And check if this bin corresponds to a neighbor cell. + is_neigh_cell = count(isc /= 0) > 0 + + ref_xyz = xyz(at + 1, :) + if (is_neigh_cell) then + ! If we are looking at a neighbouring cell in a direction + ! where there are no periodic boundary conditions, go to + ! next bin. + if (any(isc /= 0 .and. .not. pbc)) cycle + ! Otherwise, move the atom to the neighbor cell. We do this + ! instead of moving potential neighbors to the unit cell + ! because in this way we reduce the number of operations. + ref_xyz = ref_xyz - matmul(isc, cell) + end if + ! Loop through all atoms that are in this bin. ! If neigh_at == -1, this means no more atoms are in this bin. do while (neigh_at >= 0) - dist = norm2(xyz(at + 1, :) - xyz(neigh_at + 1, :)) + ! If this is a self interaction and the user didn't want them, + ! go to next atom. + if (.not. self_interaction .and. neigh_at == at) then + ! Get the next atom in this bin. Sum 1 to get fortran index. + neigh_at = list(neigh_at + 1) + cycle + end if + ! Calculate the distance between the atom and the potential + ! neighbor. + dist = norm2(ref_xyz - xyz(neigh_at + 1, :)) + + ! Get the threshold for this pair of atoms threshold = thresholds(at + 1) if (sphere_overlap) then + ! If the user wants to check for sphere overlaps, we have + ! to sum the radius of the neighbor to the threshold threshold = threshold + thresholds(neigh_at + 1) endif if (dist < threshold) then - if (.not. (.not. self_interaction .and. neigh_at == at)) then - ! Store the pair of potential neighbours. - neighs(i_pair, 1) = at - neighs(i_pair, 2) = neigh_at - ! Increment the pair index - i_pair = i_pair + 1 - end if + ! Store the pair of neighbours. + neighs(i_pair, 1) = at + neighs(i_pair, 2) = neigh_at + neighs(i_pair, 3:) = isc + ! Increment the pair index + i_pair = i_pair + 1 end if ! Get the next atom in this bin. Sum 1 to get fortran index. @@ -167,13 +227,10 @@ subroutine get_pairs(at_indices, indices, heads, list, max_npairs, & end subroutine get_pairs -subroutine get_all_unique_pairs(indices, heads, list, & - max_npairs, self_interaction, xyz, thresholds, sphere_overlap, & +subroutine get_all_unique_pairs(indices, iscs, heads, list, & + max_npairs, self_interaction, xyz, cell, pbc, thresholds, sphere_overlap, & N, neighs, npairs) - ! Gets all pairs of atoms that are neighbours according to a given - ! threshold. - ! - ! This routine removes duplicates. + ! Gets all unique pairs of atoms that are neighbours. ! ! INPUTS ! --------- @@ -182,6 +239,9 @@ subroutine get_all_unique_pairs(indices, heads, list, & ! 8 bins that contain potential neighbours. Bin indices are ! 0-indexed (python indices). The first dimension should be ! as large as the total number of atoms. + ! iscs: + ! Supercell indices of each bin index in `indices`. Supercell + ! indices are 0-indexed (python indices). ! heads: ! For each bin, the index of `list` where we can find the ! first atom that is contained in it. @@ -198,8 +258,24 @@ subroutine get_all_unique_pairs(indices, heads, list, & ! It is used to allocate the `neighs` array. This is computed ! in python with the help of the `counts` array constructed ! by `build_table`. - ! self_interaction: bool, optional + ! self_interaction: ! whether to consider an atom a neighbour of itself. + ! xyz: + ! the cartesian coordinates of all atoms in the geometry. + ! cell: + ! the lattice vectors of the geometry, where cell(i, :) is the + ! ith lattice vector. + ! pbc: + ! for each direction, whether periodic boundary conditions should + ! be considered. + ! thresholds: + ! the threshold radius for each atom in the geometry. + ! sphere_overlap: + ! If true, two atoms are considered neighbours if their spheres + ! overlap. + ! If false, two atoms are considered neighbours if the second atom + ! is within the sphere of the first atom. Note that this implies that + ! atom `i` might be atom `j`'s neighbour while the opposite is not true. ! N: ! total number of atoms. No need to pass it when calling from ! python. It will be inferred from the shapes of `indices` and @@ -208,31 +284,36 @@ subroutine get_all_unique_pairs(indices, heads, list, & ! OUTPUTS ! -------- ! neighs: - ! contains the atom indices of all unique pairs that are - ! potential neighbours. First column of atoms is sorted - ! in increasing order. Each pair is also sorted in - ! increasing order. - ! Atom indices are 0-indexed (python indices). + ! contains the atom indices of all unique neighbor pairs + ! First column of atoms is sorted in increasing order. + ! Columns 3 to 5 contain the supercell indices of the neighbour + ! atom (the one in column 2, column 1 atoms are always in the + ! unit cell). + ! All indices are 0-indexed (python indices). ! npairs: ! The number of unique pairs that have been found, which ! will be less than `max_npairs`. This should presumably ! be used to slice the `neighs` array once in python. - implicit none + integer, parameter :: dp = selected_real_kind(p=15) - integer, intent(in) :: indices(N, 8), heads(:), list(N) + integer, intent(in) :: indices(N, 8), iscs(N, 8, 3) + integer, intent(in) :: heads(:), list(N) integer, intent(in) :: max_npairs logical, intent(in) :: self_interaction - real(8), intent(in) :: xyz(N, 3), thresholds(N) + real(8), intent(in) :: xyz(N, 3), cell(3,3), thresholds(N) + logical, intent(in) :: pbc(3) logical, intent(in) :: sphere_overlap integer, intent(in) :: N - integer, intent(out) :: neighs(max_npairs, 2), npairs + integer, intent(out) :: neighs(max_npairs, 5), npairs ! Auxiliar variables. integer :: at, neigh_at integer :: bin_index, i_pair, j - real(8) :: dist, threshold + real(8) :: dist, threshold, ref_xyz(3) + integer :: isc(3) + logical :: is_neigh_cell ! Initialize the pair index, which will tell us the location ! to store each pair that we find. Each time we store a new pair @@ -248,23 +329,54 @@ subroutine get_all_unique_pairs(indices, heads, list, & ! Get the first atom index in this bin. neigh_at = heads(bin_index) - ! Loop through all atoms that are in this bin. - ! If neigh_at is smaller than at, we already stored - ! this pair when performing the search for neigh_at. - ! The following atoms will have even lower indices - ! So we can just move to the next bin. - do while (neigh_at > at) - dist = norm2(xyz(at + 1, :) - xyz(neigh_at + 1, :)) + ! If there are no atoms in this bin, do not even bother + ! checking supercell indices. + if (neigh_at == -1) cycle + + ! Find the supercell indices for this bin + isc(:) = iscs(at + 1, j, :) + ! And check if this bin corresponds to a neighbor cell. + is_neigh_cell = any(isc /= 0) + + ref_xyz = xyz(at + 1, :) + if (is_neigh_cell) then + ! If we are looking at a neighbouring cell in a direction + ! where there are no periodic boundary conditions, go to + ! next bin. + if (any(isc /= 0 .and. .not. pbc)) cycle + ! Otherwise, move the atom to the neighbor cell. We do this + ! instead of moving potential neighbors to the unit cell + ! because in this way we reduce the number of operations. + ref_xyz = ref_xyz - matmul(isc, cell) + end if + ! Loop through all atoms that are in this bin. + do while (neigh_at >= 0) + ! If neigh_at is smaller than at, we already stored + ! this pair when performing the search for neigh_at. + ! The following atoms will have even lower indices + ! So we can just move to the next bin. However, if + ! we are checking a neighboring cell, this connection + ! will always be unique. + if (.not. is_neigh_cell .and. neigh_at <= at) exit + + ! Calculate the distance between the atom and the potential + ! neighbor. + dist = norm2(ref_xyz - xyz(neigh_at + 1, :)) + + ! Get the threshold for this pair of atoms threshold = thresholds(at + 1) if (sphere_overlap) then + ! If the user wants to check for sphere overlaps, we have + ! to sum the radius of the neighbor to the threshold threshold = threshold + thresholds(neigh_at + 1) endif if (dist < threshold) then - ! Store the pair of potential neighbours. + ! Store the pair of neighbours. neighs(i_pair, 1) = at neighs(i_pair, 2) = neigh_at + neighs(i_pair, 3:) = isc ! Increment the pair index i_pair = i_pair + 1 end if @@ -276,6 +388,7 @@ subroutine get_all_unique_pairs(indices, heads, list, & if (j==1 .and. self_interaction) then ! Add the self interaction neighs(i_pair, :) = at + neighs(i_pair, 3:) = 0 ! Increment the pair index i_pair = i_pair + 1 end if diff --git a/sisl/geom/neighs/neighfinder.py b/sisl/geom/neighs/neighfinder.py index 952071bd1..d9c1096e0 100644 --- a/sisl/geom/neighs/neighfinder.py +++ b/sisl/geom/neighs/neighfinder.py @@ -56,6 +56,8 @@ def setup_finder(self, geometry=None, R=None, sphere_overlap=None): If not provided, an array is constructed, where the radius for each atom is its maxR. + + Note that negative R values are allowed sphere_overlap: bool, optional If `True`, two atoms are considered neighbours if their spheres overlap. @@ -80,10 +82,6 @@ def setup_finder(self, geometry=None, R=None, sphere_overlap=None): # Set the radius self._R = R - - # Check that the radius values make sense. - if np.min(self._R) <= 0: - raise ValueError("R contains values that are not greater than 0") # If sphere overlap was not provided, set it to False if R is a single float # and True otherwise. @@ -94,6 +92,9 @@ def setup_finder(self, geometry=None, R=None, sphere_overlap=None): # Determine the bin_size as the maximum DIAMETER to ensure that we ALWAYS # only need to look one bin away for neighbors. bin_size = np.max(self._R) * 2 + # Check that the bin size is positive. + if bin_size <= 0: + raise ValueError("All R values are 0 or less. Please provide some positive values") if sphere_overlap: # In the case when we want to check for sphere overlap, the size should # be twice as big. @@ -260,7 +261,7 @@ def _scalar_to_cartesian_index(self, index): second, first = np.divmod(index, self.nbins[0]) return np.array([first, second, third]).T - def find_neighbours(self, atoms=None, as_pairs=False, self_interaction=False): + def find_neighbours(self, atoms=None, as_pairs=False, self_interaction=False, pbc=(True, True, True)): """Find neighbours as specified in the finder. This routine only executes the action of finding neighbours, @@ -279,18 +280,26 @@ def find_neighbours(self, atoms=None, as_pairs=False, self_interaction=False): the neighbours for each atom is returned. self_interaction: bool, optional whether to consider an atom a neighbour 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 ---------- np.ndarray or list: If `as_pairs=True`: - A `np.ndarray` of shape (n_pairs, 2) is returned. Each pair `ij` - means that `j` is a neighbour of `i`. + A `np.ndarray` of shape (n_pairs, 5) is returned. + Each pair `ij` means that `j` is a neighbour of `i`. + The three extra columns are the supercell indices of atom `j`. If `as_pairs=False`: A list containing a numpy array of shape (n_neighs) for each atom. """ # Sanitize atoms atoms = self.geometry._sanitize_atoms(atoms) + + # Cast R and pbc into arrays of appropiate shape and type. + thresholds = np.full(self.geometry.na, self._R, dtype=np.float64) + pbc = np.full(3, pbc, dtype=np.bool) # Get search indices search_indices, isc = self._get_search_indices(self.geometry.fxyz[atoms], cartesian=False) @@ -301,25 +310,27 @@ def find_neighbours(self, atoms=None, as_pairs=False, self_interaction=False): total_pairs = at_counts.sum() if not self_interaction: total_pairs -= search_indices.shape[0] - - thresholds = np.empty(self.geometry.na ,dtype=np.float64) - thresholds[:] = self._R # Find the candidate pairs neighbour_pairs, split_ind = _fneighs.get_pairs( - atoms, search_indices, self._heads, self._list, total_pairs, self_interaction, - self.geometry.xyz, thresholds, self._sphere_overlap + atoms, search_indices, isc, self._heads, self._list, total_pairs, self_interaction, + self.geometry.xyz, self.geometry.cell, pbc, thresholds, self._sphere_overlap ) if as_pairs: - # Just returned the filtered pairs - return neighbour_pairs[split_ind[-1]] + # Just return the neighbour pairs + return neighbour_pairs[:split_ind[-1]] else: # Split to get the neighbours of each atom - return np.split(neighbour_pairs[:, 1], split_ind)[:-1] + return np.split(neighbour_pairs[:, 1:], split_ind, axis=0)[:-1] - def find_all_unique_pairs(self, self_interaction=False): + def find_all_unique_pairs(self, self_interaction=False, pbc=(True, True, True)): """Find all unique neighbour 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 + atom 0 (unit cell) to atom 5 (1, 0, 0) is not the same as the + connection atom 5 (unit cell) to atom 0 (-1, 0, 0). Note that this routine can not be called if `sphere_overlap` is set to `False` and the radius is not a single float. In that case, @@ -338,23 +349,27 @@ def find_all_unique_pairs(self, self_interaction=False): the neighbours for each atom is returned. self_interaction: bool, optional whether to consider an atom a neighbour 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 ---------- - np.ndarray or list: - If `as_pairs=True`: - A `np.ndarray` of shape (n_pairs, 2) is returned. Each pair `ij` - means that `j` is a neighbour of `i`. - If `as_pairs=False`: - A list containing a numpy array of shape (n_neighs) for each atom. + np.ndarray of shape (n_pairs, 5): + Each pair `ij` means that `j` is a neighbour of `i`. + The three extra columns are the supercell indices of atom `j`. """ if not self._sphere_overlap and not isinstance(self._R, Real): raise ValueError("Unique atom pairs do not make sense if we are not looking for sphere overlaps." " Please setup the finder again setting `sphere_overlap` to `True` if you wish so.") + + # Cast R and pbc into arrays of appropiate shape and type. + thresholds = np.full(self.geometry.na, self._R, dtype=np.float64) + pbc = np.full(3, pbc, dtype=np.bool) # Get search indices search_indices, isc = self._get_search_indices(self.geometry.fxyz, cartesian=False) - + # Get atom counts at_counts = self._get_search_atom_counts(search_indices) @@ -362,13 +377,10 @@ def find_all_unique_pairs(self, self_interaction=False): if not self_interaction: max_pairs -= search_indices.shape[0] - thresholds = np.empty(self.geometry.na ,dtype=np.float64) - thresholds[:] = self._R - # Find the candidate pairs candidate_pairs, n_pairs = _fneighs.get_all_unique_pairs( - search_indices, self._heads, self._list, max_pairs, self_interaction, - self.geometry.xyz, thresholds, self._sphere_overlap + search_indices, isc, self._heads, self._list, max_pairs, self_interaction, + self.geometry.xyz, self.geometry.cell, pbc, thresholds, self._sphere_overlap ) return candidate_pairs[:n_pairs]