diff --git a/src/sisl/_core/geometry.py b/src/sisl/_core/geometry.py
index b0ecf710e6..38a0b6675a 100644
--- a/src/sisl/_core/geometry.py
+++ b/src/sisl/_core/geometry.py
@@ -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
@@ -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)
 
@@ -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
@@ -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)
diff --git a/src/sisl/_core/tests/test_geometry.py b/src/sisl/_core/tests/test_geometry.py
index 257b9583b8..081cba296a 100644
--- a/src/sisl/_core/tests/test_geometry.py
+++ b/src/sisl/_core/tests/test_geometry.py
@@ -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