diff --git a/CHANGELOG.md b/CHANGELOG.md index 8b6d39572b..136fac2e4e 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/src/sisl/geom/_common.py b/src/sisl/geom/_common.py index 83aa1da4c8..dd726a7704 100644 --- a/src/sisl/geom/_common.py +++ b/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: diff --git a/src/sisl/geometry.py b/src/sisl/geometry.py index e06e6ceabe..04c8eca3d6 100644 --- a/src/sisl/geometry.py +++ b/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: diff --git a/src/sisl/sparse_geometry.py b/src/sisl/sparse_geometry.py index bd414a104f..c5b0364d01 100644 --- a/src/sisl/sparse_geometry.py +++ b/src/sisl/sparse_geometry.py @@ -617,7 +617,7 @@ def iter_nnz(self): __iter__ = iter_nnz - def create_construct(self, R, param): + def create_construct(self, R, params): """ Create a simple function for passing to the `construct` function. This is simply to leviate the creation of simplistic @@ -627,7 +627,7 @@ def create_construct(self, R, param): >>> def func(self, ia, atoms, atoms_xyz=None): ... idx = self.geometry.close(ia, R=R, atoms=atoms, atoms_xyz=atoms_xyz) - ... for ix, p in zip(idx, param): + ... for ix, p in zip(idx, params): ... self[ia, ix] = p Notes @@ -640,24 +640,26 @@ def create_construct(self, R, param): ---------- R : array_like radii parameters for different shells. - Must have same length as `param` or one less. + Must have same length as `params` or one less. If one less it will be extended with ``R[0]/100`` - param : array_like + params : array_like coupling constants corresponding to the `R` - ranges. ``param[0,:]`` are the elements + ranges. ``params[0, :]`` are the elements for the all atoms within ``R[0]`` of each atom. See Also -------- construct : routine to create the sparse matrix from a generic function (as returned from `create_construct`) """ - if len(R) != len(param): + if len(R) != len(params): raise ValueError(f"{self.__class__.__name__}.create_construct got different lengths of `R` and `param`") def func(self, ia, atoms, atoms_xyz=None): idx = self.geometry.close(ia, R=R, atoms=atoms, atoms_xyz=atoms_xyz) - for ix, p in zip(idx, param): + for ix, p in zip(idx, params): self[ia, ix] = p + func.R = R + func.params = params return func @@ -705,7 +707,6 @@ def construct(self, func, na_iR=1000, method='rand', eta=None): tile : tiling *after* construct is much faster for very large systems repeat : repeating *after* construct is much faster for very large systems """ - if not callable(func): if not isinstance(func, (tuple, list)): raise ValueError('Passed `func` which is not a function, nor tuple/list of `R, param`') @@ -718,17 +719,29 @@ def construct(self, func, na_iR=1000, method='rand', eta=None): # Convert to a proper function func = self.create_construct(func[0], func[1]) - iR = self.geometry.iR(na_iR) + try: + # if the function was created through `create_construct`, then + # we have access to the radii used. + R = func.R + try: + if len(R) > 0: + R = R[-1] + except TypeError: + pass + except AttributeError: + R = None + + iR = self.geometry.iR(na_iR, R=R) # Create eta-object - eta = progressbar(self.na, self.__class__.__name__ + '.construct', 'atom', eta) + eta = progressbar(self.na, f"{self.__class__.__name__ }.construct", "atom", eta) # Do the loop - for ias, idxs in self.geometry.iter_block(iR=iR, method=method): + for ias, idxs in self.geometry.iter_block(iR=iR, method=method, R=R): # Get all the indexed atoms... # This speeds up the searching for coordinates... - idxs_xyz = self.geometry[idxs, :] + idxs_xyz = self.geometry[idxs] # Loop the atoms inside for ia in ias: