Merge pull request #634 from zerothi/633-fix

fixed #633 due to weird combination of boundary lookups

CHANGELOG.md

@@ -22,6 +22,9 @@ we hit release version 1.0.0.
   The _info_attributes_ contains a list of attributes that can be
   discovered while reading ascii files see #509
 
+### Fixed
+- fixed cases where `Geometry.close` would not catch all neighbours, #633
+
 ### Changed
 - `Lattice` now holds the boundary conditions (not `Grid`), see #626
 - Some siles exposed certain properties containing basic information

src/sisl/geom/_common.py

@@ -10,7 +10,7 @@ def geometry_define_nsc(geometry, periodic=(True, True, True)):
     """Define the number of supercells for a geometry based on the periodicity """
     if np.all(geometry.maxR(True) > 0.):
         geometry.optimize_nsc()
-        for i, d, per in zip(range(3), "abc", periodic):
+        for d, per in zip("abc", periodic):
             if per:
                 geometry.lattice.set_boundary_condition(**{d: "Periodic"})
             else:

src/sisl/geometry.py

@@ -195,9 +195,9 @@ def __init__(self, xyz: ArrayLike, atoms=None, lattice=None, names=None):
         else:
             self._names = NamedIndex(names)
 
-        self.__init_lattice(lattice)
+        self._init_lattice(lattice)
 
-    def __init_lattice(self, lattice):
+    def _init_lattice(self, lattice):
         """ Initializes the supercell by *calculating* the size if not supplied
 
         If the supercell has not been passed we estimate the unit cell size
@@ -892,7 +892,7 @@ def iR(self, na: int=1000, iR: int=20, R: Optional[float]=None) -> int:
 
         # default block iterator
         if R is None:
-            R = self.maxR()
+            R = self.maxR() + 0.001
         if R < 0:
             raise ValueError(f"{self.__class__.__name__}.iR unable to determine a number of atoms within a sphere with negative radius, is maxR() defined?")
 
@@ -910,13 +910,12 @@ def iter_block_rand(self, iR=20, R=None, atoms: Optional[AtomsArgument]=None) ->
         # We implement yields as we can then do nested iterators
         # create a boolean array
         na = len(self)
-        not_passed = np.empty(na, dtype='b')
         if atoms is not None:
+            not_passed = np.zeros(na, dtype=bool)
             # Reverse the values
-            not_passed[:] = False
             not_passed[atoms] = True
         else:
-            not_passed[:] = True
+            not_passed = np.ones(na, dtype=bool)
 
         # Figure out how many we need to loop on
         not_passed_N = np.sum(not_passed)
@@ -925,11 +924,9 @@ def iter_block_rand(self, iR=20, R=None, atoms: Optional[AtomsArgument]=None) ->
             raise SislError(f'{self.__class__.__name__}.iter_block_rand too small iR!')
 
         if R is None:
-            R = self.maxR()
+            R = self.maxR() + 0.001
         # The boundaries (ensure complete overlap)
-        R = np.array([iR - 0.975, iR + .025]) * R
-
-        append = np.append
+        R = np.array([iR - 0.5, iR + 0.501]) * R
 
         # loop until all passed are true
         while not_passed_N > 0:
@@ -942,6 +939,7 @@ def iter_block_rand(self, iR=20, R=None, atoms: Optional[AtomsArgument]=None) ->
             # atoms at any single time.
             # Shuffling will cut down needed iterations.
             np.random.shuffle(all_true)
+            # take one element, after shufling, we can take the first
             idx = all_true[0]
             del all_true
 
@@ -950,15 +948,20 @@ def iter_block_rand(self, iR=20, R=None, atoms: Optional[AtomsArgument]=None) ->
 
             # get all elements within two radii
             all_idx = self.close(idx, R=R)
-            # Get unit-cell atoms
-            all_idx[0] = self.sc2uc(all_idx[0], unique=True)
-            # First extend the search-space (before reducing)
-            all_idx[1] = self.sc2uc(append(all_idx[1], all_idx[0]), unique=True)
+
+            # Get unit-cell atoms, we are drawing a circle, and this
+            # circle only encompasses those already in the unit-cell.
+            all_idx[1] = np.union1d(self.sc2uc(all_idx[0], unique=True),
+                                    self.sc2uc(all_idx[1], unique=True))
+            # If we translated stuff into the unit-cell, we could end up in situations
+            # where the supercell atom is in the circle, but not the UC-equivalent
+            # of that one.
+            all_idx[0] = all_idx[0][all_idx[0] < na]
 
             # Only select those who have not been runned yet
             all_idx[0] = all_idx[0][not_passed[all_idx[0]].nonzero()[0]]
             if len(all_idx[0]) == 0:
-                raise SislError('Internal error, please report to the developers')
+                continue
 
             # Tell the next loop to skip those passed
             not_passed[all_idx[0]] = False
@@ -994,18 +997,18 @@ def iter_block_shape(self, shape=None, iR=20, atoms: Optional[AtomsArgument]=Non
         if iR < 2:
             raise SislError(f'{self.__class__.__name__}.iter_block_shape too small iR!')
 
-        R = self.maxR()
+        R = self.maxR() + 0.001
         if shape is None:
             # we default to the Cube shapes
-            dS = (Cube(R * (iR - 1.975)),
-                  Cube(R * (iR + 0.025)))
+            dS = (Cube((iR - 0.5) * R),
+                  Cube((iR + 1.501) * R))
         else:
             if isinstance(shape, Shape):
                 dS = (shape,)
             else:
                 dS = tuple(shape)
             if len(dS) == 1:
-                dS += (dS[0].expand(R + 0.01), )
+                dS += (dS[0].expand(R), )
         if len(dS) != 2:
             raise ValueError(f'{self.__class__.__name__}.iter_block_shape, number of Shapes *must* be one or two')
 
@@ -1054,7 +1057,6 @@ def iter_block_shape(self, shape=None, iR=20, atoms: Optional[AtomsArgument]=Non
 
         # Shorthand function
         where = np.where
-        append = np.append
 
         # Now we loop in each direction
         for x, y, z in product(range(ixyz[0]), range(ixyz[1]), range(ixyz[2])):
@@ -1070,10 +1072,14 @@ def iter_block_shape(self, shape=None, iR=20, atoms: Optional[AtomsArgument]=Non
             # get all elements within two radii
             all_idx = self.within(dS)
 
-            # Get unit-cell atoms
-            all_idx[0] = self.sc2uc(all_idx[0], unique=True)
-            # First extend the search-space (before reducing)
-            all_idx[1] = self.sc2uc(append(all_idx[1], all_idx[0]), unique=True)
+            # Get unit-cell atoms, we are drawing a circle, and this
+            # circle only encompasses those already in the unit-cell.
+            all_idx[1] = np.union1d(self.sc2uc(all_idx[0], unique=True),
+                                    self.sc2uc(all_idx[1], unique=True))
+            # If we translated stuff into the unit-cell, we could end up in situations
+            # where the supercell atom is in the circle, but not the UC-equivalent
+            # of that one.
+            all_idx[0] = all_idx[0][all_idx[0] < na]
 
             # Only select those who have not been runned yet
             all_idx[0] = all_idx[0][not_passed[all_idx[0]].nonzero()[0]]
@@ -1091,9 +1097,9 @@ def iter_block_shape(self, shape=None, iR=20, atoms: Optional[AtomsArgument]=Non
             yield all_idx[0], all_idx[1]
 
         if np.any(not_passed):
-            print(not_passed.nonzero()[0])
-            print(np.sum(not_passed), len(self))
-            raise SislError(f'{self.__class__.__name__}.iter_block_shape error on iterations. Not all atoms have been visited.')
+            not_passed = not_passed.nonzero()[0]
+            raise SislError(f"{self.__class__.__name__}.iter_block_shape error on iterations. Not all atoms have been visited "
+                            f"{not_passed}")
 
     def iter_block(self, iR: int=20, R: float=None,
                    atoms: Optional[AtomsArgument]=None,
@@ -1118,7 +1124,7 @@ def iter_block(self, iR: int=20, R: float=None,
         iR :
             the number of `R` ranges taken into account when doing the iterator
         R :
-            enables overwriting the local R quantity. Defaults to ``self.maxR()``
+            enables overwriting the local R quantity. Defaults to ``self.maxR() + 0.001``
         atoms :
             enables only effectively looping a subset of the full geometry
         method : {'rand', 'sphere', 'cube'}
@@ -1140,21 +1146,24 @@ def iter_block(self, iR: int=20, R: float=None,
             raise SislError(f'{self.__class__.__name__}.iter_block too small iR!')
 
         method = method.lower()
-        if method == 'rand' or method == 'random':
+        if method in ("rand", "random"):
             yield from self.iter_block_rand(iR, R, atoms)
-        else:
+        elif method in ("sphere", "cube"):
             if R is None:
-                R = self.maxR()
+                R = self.maxR() + 0.001
 
             # Create shapes
             if method == 'sphere':
-                dS = (Sphere(R * (iR - 0.975)),
-                      Sphere(R * (iR + 0.025)))
+                dS = (Sphere((iR - 0.5) * R),
+                      Sphere((iR + 0.501) * R))
             elif method == 'cube':
-                dS = (Cube(R * (2 * iR - 0.975)),
-                      Cube(R * (2 * iR + 0.025)))
+                dS = (Cube((2 * iR - 0.5) * R),
+                      # we need an extra R here since it needs to extend on both sides
+                      Cube((2 * iR + 1.501) * R))
 
             yield from self.iter_block_shape(dS)
+        else:
+            raise ValueError(f"{self.__class__.__name__}.iter_block got unexpected 'method' argument: {method}")
 
     def copy(self) -> Geometry:
         """Create a new object with the same content (a copy)."""
@@ -1737,10 +1746,10 @@ def optimize_nsc(self, axis=None,
             axis = _a.asarrayi(axis).ravel()
 
         if R is None:
-            R = self.maxR()
+            R = self.maxR() + 0.001
         if R < 0:
             R = 0.00001
-            warn(self.__class__.__name__ +
+            warn(f"{self.__class__.__name__}"
                  ".optimize_nsc could not determine the radius from the "
                  "internal atoms (defaulting to zero radius).")
 
@@ -2918,7 +2927,7 @@ def add_vacuum(self,
         """
         new = self.copy()
         new.set_lattice(self.lattice.add_vacuum(vacuum, axis))
-        if vacuum > self.maxR():
+        if vacuum > self.maxR() + 0.001:
             # only overwrite along axis
             nsc = [None for _ in range(3)]
             nsc[axis] = 1
@@ -3243,7 +3252,7 @@ def axyz(self, atoms: Optional[AtomsArgument]=None, isc=None) -> ndarray:
             # get offsets from atomic indices (note that this will be per atom)
             isc = self.a2isc(atoms)
             offset = self.lattice.offset(isc)
-            return self.xyz[self.sc2uc(atoms), :] + offset
+            return self.xyz[self.sc2uc(atoms)] + offset
 
         # Neither of atoms, or isc are `None`, we add the offset to all coordinates
         return self.axyz(atoms) + self.lattice.offset(isc)
@@ -3500,7 +3509,7 @@ def close_sc(self,
             If a single float it will return ``x <= R``.
         atoms :
             List of atoms that will be considered. This can
-            be used to only take out a certain atoms.
+            be used to only take out a certain atom.
         atoms_xyz : array_like of float, optional
             The atomic coordinates of the equivalent `atoms` variable (`atoms` must also be passed)
         ret_xyz :
@@ -3519,13 +3528,19 @@ def close_sc(self,
         rij
             distance of the indexed atoms to the center coordinate (only for true `ret_rij`)
         """
+        maxR = self.maxR() + 0.001
         if R is None:
-            R = np.array([self.maxR()], np.float64)
+            R = np.array([maxR], np.float64)
         elif not isndarray(R):
             R = _a.asarrayd(R).ravel()
 
         # Maximum distance queried
         max_R = R[-1]
+        if atoms is not None and max_R > maxR:
+            warn(f"{self.__class__.__name__}.close_sc has been passed an 'atoms' argument "
+                 "together with an R value larger than the orbital ranges. "
+                 "If used together with 'sparse-matrix.construct' this can result in wrong couplings.",
+                 register=True)
 
         # Convert to actual array
         if atoms is not None:
@@ -3535,25 +3550,25 @@ def close_sc(self,
             atoms_xyz = None
 
         if isinstance(xyz_ia, Integral):
-            off = self.xyz[xyz_ia, :]
+            off = self.xyz[xyz_ia]
         elif not isndarray(xyz_ia):
             off = _a.asarrayd(xyz_ia)
         elif xyz_ia.ndim == 0:
-            off = self.xyz[xyz_ia, :]
+            off = self.xyz[xyz_ia]
         else:
             off = xyz_ia
 
         # Calculate the complete offset
-        foff = self.lattice.offset(isc)[:] - off[:]
+        foff = self.lattice.offset(isc) - off
 
-        # Get atomic coordinate in principal cell
+        # Get distances between `xyz_ia` and `atoms`
         if atoms_xyz is None:
-            dxa = self.axyz(atoms) + foff.reshape(1, 3)
+            dxa = self.axyz(atoms) + foff
         else:
             # For extremely large systems re-using the
             # atoms_xyz is faster than indexing
             # a very large array
-            dxa = atoms_xyz + foff.reshape(1, 3)
+            dxa = atoms_xyz + foff
 
         # Immediately downscale by easy checking
         # This will reduce the computation of the vector-norm
@@ -3564,20 +3579,20 @@ def close_sc(self,
         # method..
         if atoms is None:
             atoms, d = indices_in_sphere_with_dist(dxa, max_R)
-            dxa = dxa[atoms, :].reshape(-1, 3)
+            dxa = dxa[atoms].reshape(-1, 3)
         else:
             ix, d = indices_in_sphere_with_dist(dxa, max_R)
             atoms = atoms[ix]
-            dxa = dxa[ix, :].reshape(-1, 3)
+            dxa = dxa[ix].reshape(-1, 3)
             del ix
 
         if len(atoms) == 0:
             # Create default return
-            ret = [[_a.emptyi([0])] * len(R)]
+            ret = [[_a.emptyi([0]) for _ in R]]
             if ret_xyz:
-                ret.append([_a.emptyd([0, 3])] * len(R))
+                ret.append([_a.emptyd([0, 3]) for _ in R])
             if ret_rij:
-                ret.append([_a.emptyd([0])] * len(R))
+                ret.append([_a.emptyd([0]) for _ in R])
 
             # Quick return if there are
             # no entries...
@@ -3592,7 +3607,7 @@ def close_sc(self,
             return ret[0]
 
         if ret_xyz:
-            xa = dxa[:, :] + off[None, :]
+            xa = dxa + off
         del dxa  # just because this array could be very big...
 
         # Check whether we only have one range to check.
@@ -3889,13 +3904,13 @@ def close(self,
             integer lattice offsets for the couplings (related to `rij` without atomic coordinates)
         """
         if R is None:
-            R = self.maxR()
+            R = self.maxR() + 0.001
         R = _a.asarrayd(R).ravel()
         nR = R.size
 
         # Convert index coordinate to point
         if isinstance(xyz_ia, Integral):
-            xyz_ia = self.xyz[xyz_ia, :]
+            xyz_ia = self.xyz[xyz_ia]
         elif not isndarray(xyz_ia):
             xyz_ia = _a.asarrayd(xyz_ia)
 
@@ -3924,7 +3939,7 @@ def isc_tile(isc, n):
         for s in range(self.n_s):
 
             na = self.na * s
-            isc = self.lattice.sc_off[s, :]
+            isc = self.lattice.sc_off[s]
             sret = self.close_sc(xyz_ia, isc, R=R,
                                  atoms=atoms, atoms_xyz=atoms_xyz,
                                  ret_xyz=ret_xyz, ret_rij=ret_rij)
@@ -4398,7 +4413,7 @@ def sparserij(self, dtype=np.float64, na_iR: int=1000, method='rand'):
         rij = SparseAtom(self, nnzpr=20, dtype=dtype)
 
         # Get R
-        R = (0.1, self.maxR())
+        R = (0.1, self.maxR() + 0.001)
         iR = self.iR(na_iR)
 
         # Do the loop
@@ -4495,9 +4510,9 @@ def distance(self,
                 off = self.lattice.offset(sc)
 
                 for ii, jj, kk in product([0, 1], [0, 1], [0, 1]):
-                    o = self.cell[0, :] * ii + \
-                        self.cell[1, :] * jj + \
-                        self.cell[2, :] * kk
+                    o = self.cell[0] * ii + \
+                        self.cell[1] * jj + \
+                        self.cell[2] * kk
                     maxR = max(maxR, fnorm(off + o))
 
             if R > maxR: