Skip to content

Commit

Permalink
fixed #633 due to weird combination of boundary lookups
Browse files Browse the repository at this point in the history
The problem of close arose when the searched sphere
was just crossing a border of the unit-cell (or supercell).
In those cases could the translated to unit-cell indices
miss a neighbour.

The fix was rather simple; once the atoms closests
has been found, only loop over those in the primary unit-cell.

There are many other small changes here, but they
are more cosmetic and should provide more stability
when using orbital distances very close to atomic coordinates.

Signed-off-by: Nick Papior <[email protected]>
  • Loading branch information
zerothi committed Oct 30, 2023
1 parent ef5eb1e commit 223898c
Show file tree
Hide file tree
Showing 3 changed files with 67 additions and 55 deletions.
3 changes: 3 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,9 @@ we hit release version 1.0.0.
debugging.
- marked all `toSphere|toEllipsoid|...` as deprecated

### Fixed
- fixed cases where `Geometry.close` would not catch all neighbours, #633

### Changed
- `Lattice` now holds the boundary conditions (not `Grid`), see #626

Expand Down
117 changes: 63 additions & 54 deletions src/sisl/geometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -196,9 +196,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
Expand Down Expand Up @@ -840,7 +840,7 @@ def iR(self, na=1000, iR=20, R=None):

# default block iterator
if R is None:
R = self.maxR()
R = self.maxR() * 1.01
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?")

Expand All @@ -858,13 +858,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)
Expand All @@ -873,11 +872,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.01
# 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:
Expand All @@ -890,6 +887,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

Expand All @@ -898,15 +896,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
Expand Down Expand Up @@ -942,18 +945,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.01
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')

Expand Down Expand Up @@ -1002,7 +1005,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])):
Expand All @@ -1018,10 +1020,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]]
Expand Down Expand Up @@ -1064,7 +1070,7 @@ def iter_block(self, iR=20, R=None, atoms: Optional[AtomsArgument]=None, method:
iR : int, optional
the number of `R` ranges taken into account when doing the iterator
R : float, optional
enables overwriting the local R quantity. Defaults to ``self.maxR()``
enables overwriting the local R quantity. Defaults to ``self.maxR() + 0.01``
atoms : array_like, optional
enables only effectively looping a subset of the full geometry
method : {'rand', 'sphere', 'cube'}
Expand All @@ -1086,21 +1092,24 @@ def iter_block(self, iR=20, R=None, atoms: Optional[AtomsArgument]=None, method:
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.01

# 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:
""" A copy of the object. """
Expand Down Expand Up @@ -3138,7 +3147,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)
Expand Down Expand Up @@ -3385,7 +3394,7 @@ def close_sc(self, xyz_ia, isc=(0, 0, 0), R=None,
If a single float it will return ``x <= R``.
atoms : array_like of int, optional
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 : bool, optional
Expand All @@ -3405,7 +3414,7 @@ def close_sc(self, xyz_ia, isc=(0, 0, 0), R=None,
distance of the indexed atoms to the center coordinate (only for true `ret_rij`)
"""
if R is None:
R = np.array([self.maxR()], np.float64)
R = np.array([self.maxR() + 0.01], np.float64)
elif not isndarray(R):
R = _a.asarrayd(R).ravel()

Expand All @@ -3420,25 +3429,25 @@ def close_sc(self, xyz_ia, isc=(0, 0, 0), R=None,
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
Expand All @@ -3449,20 +3458,20 @@ def close_sc(self, xyz_ia, isc=(0, 0, 0), R=None,
# 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...
Expand All @@ -3477,7 +3486,7 @@ def close_sc(self, xyz_ia, isc=(0, 0, 0), R=None,
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.
Expand Down Expand Up @@ -3762,13 +3771,13 @@ def close(self, xyz_ia, R=None,
integer lattice offsets for the couplings (related to `rij` without atomic coordinates)
"""
if R is None:
R = self.maxR()
R = self.maxR() + 0.01
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)

Expand Down Expand Up @@ -3797,7 +3806,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)
Expand Down Expand Up @@ -4239,7 +4248,7 @@ def sparserij(self, dtype=np.float64, na_iR=1000, method='rand'):
rij = SparseAtom(self, nnzpr=20, dtype=dtype)

# Get R
R = (0.1, self.maxR())
R = (0.1, self.maxR() + 0.01)
iR = self.iR(na_iR)

# Do the loop
Expand Down Expand Up @@ -4335,9 +4344,9 @@ def distance(self, atoms: Optional[AtomsArgument]=None,
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:
Expand Down
2 changes: 1 addition & 1 deletion src/sisl/sparse_geometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -721,7 +721,7 @@ def construct(self, func, na_iR=1000, method='rand', eta=None):
iR = self.geometry.iR(na_iR)

# 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):
Expand Down

0 comments on commit 223898c

Please sign in to comment.