Skip to content

Commit

Permalink
final correction that passes tests
Browse files Browse the repository at this point in the history
Signed-off-by: Nick Papior <[email protected]>
  • Loading branch information
zerothi committed Oct 24, 2024
1 parent a332378 commit e774c5a
Show file tree
Hide file tree
Showing 2 changed files with 31 additions and 10 deletions.
25 changes: 15 additions & 10 deletions src/sisl/_core/geometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -3627,14 +3627,16 @@ def within_inf(
# 1. Number of times each lattice vector must be expanded to fit
# inside the "possibly" larger `lattice`.
idx = lattice.cell @ self.icell.T
tile_min = floor(idx.min(0))
tile_min = floor(idx.min(0)).astype(dtype=int32)
tile_max = ceil(idx.max(0)).astype(dtype=int32)

# Intrinsic offset (when atomic coordinates are outside primary unit-cell)
fxyz = np.round(self.fxyz, decimals=5)
fxyz_floor = floor(fxyz).astype(dtype=int32)
tile_min = np.minimum(tile_min, fxyz_floor.min(0)).astype(dtype=int32)
tile_max = np.maximum(tile_max, ceil(idx.max(0))).astype(dtype=int32)
# We don't collapse this as it is necessary for correcting isc further below
fxyz_ifloor = floor(fxyz).astype(dtype=int32)
fxyz_iceil = ceil(fxyz).max(0).astype(dtype=int32)
tile_min = np.minimum(tile_min, fxyz_ifloor.min(0))
tile_max = np.maximum(tile_max, fxyz_iceil)
del idx, fxyz

# 1a) correct for origin displacement
Expand Down Expand Up @@ -3662,7 +3664,14 @@ def within_inf(

# Make sure that full_geom doesn't return coordinates outside the unit cell
# for non periodic directions
nsc = full_geom.nsc.copy()
nsc = full_geom.nsc.copy() // 2

# If we have atoms outside the primary unit-cell in the original
# cell, then we should consider an nsc large enough to encompass this
nsc = np.maximum(nsc, fxyz_iceil)
nsc = np.maximum(nsc, -fxyz_ifloor.min(0))
nsc = nsc * 2 + 1

nsc[non_periodic] = 1
full_geom.set_nsc(nsc)

Expand Down Expand Up @@ -3695,7 +3704,7 @@ def within_inf(
# Convert indices to unit-cell indices and also return coordinates and
# infinite supercell indices
ia = self.asc2uc(idx)
return ia, xyz, isc - fxyz_floor[ia]
return ia, xyz, isc - fxyz_ifloor[ia]

def _orbital_values(
self, grid_shape: tuple[int, int, int], truncate_with_nsc: bool = False
Expand Down Expand Up @@ -3736,10 +3745,6 @@ def _orbital_values(
IA, XYZ, ISC = self.within_inf(lattice, periodic=self.pbc)
XYZ -= self.lattice.origin.reshape(1, 3)

# within_inf translates atoms to the unit cell to compute
# supercell indices. Here we revert that
ISC -= np.floor(self.fxyz[IA]).astype(int32)

# Don't consider atoms that are outside of the geometry's auxiliary cell.
if truncate_with_nsc:
mask = (abs(ISC) <= self.nsc // 2).all(axis=1)
Expand Down
16 changes: 16 additions & 0 deletions src/sisl/_core/tests/test_geometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -982,6 +982,22 @@ def test_within_inf_duplicates(self, setup):
lattice_3x3 = g.lattice.tile(3, 0).tile(3, 1)
assert len(g.within_inf(lattice_3x3)[0]) == 25

def test_within_inf_gh649(self):
# see https://github.com/zerothi/sisl/issues/649

# Create a geometry with an atom outside of the unit cell
geometry = Geometry([-0.5, 0, 0], lattice=np.diag([2, 10, 10]))

search = Lattice(np.diag([3, 10, 10]))
ia, xyz, isc = geometry.within_inf(search, periodic=True)
assert np.allclose(ia, 0)
assert np.allclose(isc, [1, 0, 0])

search = Lattice(np.diag([2, 10, 10]))
ia, xyz, isc = geometry.within_inf(search, periodic=True)
assert np.allclose(ia, 0)
assert np.allclose(isc, [1, 0, 0])

def test_close_sizes(self, setup):
point = 0

Expand Down

0 comments on commit e774c5a

Please sign in to comment.