Skip to content

Commit

Permalink
enh: added (optional) periodic bc
Browse files Browse the repository at this point in the history
  • Loading branch information
@pfebrer
pfebrer committed Nov 8, 2021
1 parent 723dc69 commit 45b7edf
Show file tree
Hide file tree
Showing 2 changed files with 198 additions and 73 deletions.
203 changes: 158 additions & 45 deletions sisl/geom/neighs/fneighs.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
! ---------
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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.
Expand All @@ -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
! ---------
Expand All @@ -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.
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
Loading

0 comments on commit 45b7edf

Please sign in to comment.