Skip to content

Commit

Permalink
prevent double-counting of atom pairs in the same block
Browse files Browse the repository at this point in the history
  • Loading branch information
aoterodelaroza committed Mar 26, 2024
1 parent 73cdde7 commit 39a6b8b
Showing 1 changed file with 9 additions and 1 deletion.
10 changes: 9 additions & 1 deletion src/[email protected]
Original file line number Diff line number Diff line change
Expand Up @@ -757,6 +757,7 @@ module subroutine find_asterisms_covalent(c)
integer :: iblock_stride
integer :: ix, iy, iz, jx, jy, jz, iid, jid, iidx(3), jidx(3)
integer :: iat, jat, iaux(3), lvecx(3)
logical :: sameblock

! return if there are no atoms
if (c%ncel == 0) return
Expand Down Expand Up @@ -834,16 +835,23 @@ module subroutine find_asterisms_covalent(c)
else
iaux = (/jx,jy,jz/)
jidx = (/modulo(jx,c%nblock(1)), modulo(jy,c%nblock(2)), modulo(jz,c%nblock(3))/)
jid = (jidx(1) * iblock_stride + jidx(2)) * iblock_stride + jidx(3)
if (jid < iid) cycle
xdelta = real((iaux - jidx) / c%nblock,8)
xdelta = matmul(c%m_xr2c,xdelta)
end if
if (c%iblock0(jidx(1),jidx(2),jidx(3)) == 0) cycle
sameblock = (ix == jx) .and. (iy == jy) .and. (iz == jz)

jat = c%iblock0(jidx(1),jidx(2),jidx(3))
do while (jat /= 0)
xj = c%atcel(jat)%rxc + xdelta

iat = c%iblock0(iidx(1),iidx(2),iidx(3))
if (sameblock) then
iat = c%atcel(jat)%inext
else
iat = c%iblock0(iidx(1),iidx(2),iidx(3))
end if
do while (iat /= 0)
xi = c%atcel(iat)%rxc

Expand Down

0 comments on commit 39a6b8b

Please sign in to comment.