From 1ef16ed437174d84c510f6872f01f0613ffcbe9a Mon Sep 17 00:00:00 2001 From: Nick Papior Date: Thu, 13 Jul 2023 10:09:28 +0200 Subject: [PATCH] fixed lattice handling in Grid Previously one could tile non-commensurate grids. A commensurate grid has a lattice and a geometry that are integer repetitions of the original geometry. This is now checked before one can tile the grid. The Grid.__str__ has been updated to highlight this configuration where the commensurability is shown (instead of lattice). However, when they are not commensurate, the full lattice + full geometry is shown. Removed some set_lattice when calling this on LatticeChild. It had the side-effect of setting the lattice for some hamiltonians etc when calling wavefunction. Signed-off-by: Nick Papior --- CHANGELOG.md | 1 + src/sisl/grid.py | 59 +++++++++++++++++++++++++++--------- src/sisl/lattice.py | 29 +++++------------- src/sisl/physics/electron.py | 12 +++++--- src/sisl/tests/test_grid.py | 35 ++++++++++++++++++++- src/sisl/viz/__init__.py | 1 + 6 files changed, 96 insertions(+), 41 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ad93cf2e0e..e00d88f7de 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -82,6 +82,7 @@ we hit release version 1.0.0. ### Fixed +- tiling Grid now only possible for commensurate grids (grid.lattice % grid.geometry.lattice) - rare cases for non-Gamma calculations with actual Gamma matrices resulted in crashes #572 - `MonkhorstPack.replace` now checks for symmetry k-points if the BZ is using diff --git a/src/sisl/grid.py b/src/sisl/grid.py index 3d2379e8da..bbeda6fea7 100644 --- a/src/sisl/grid.py +++ b/src/sisl/grid.py @@ -15,7 +15,7 @@ from ._internal import set_module from .geometry import Geometry from .lattice import LatticeChild -from .messages import deprecate_argument +from .messages import SislError, deprecate_argument from .shape import Shape from .utils import ( cmd, @@ -65,6 +65,14 @@ class Grid(LatticeChild): >>> grid = Grid(0.1, lattice=10, dtype=np.complex128) >>> grid == grid1 False + + It is possible to provide a geometry *and* a different lattice to make a smaller (or bigger) lattice + based on a geometry. This might be useful when creating wavefunctions or expanding densities to grids. + Here we create a square grid based on a hexagonal graphene lattice. Expanding wavefunctions from + this ``geometry`` will automatically convert to the ``lattice`` size. + >>> lattice = Lattice(10) # square lattice 10x10x10 Ang + >>> geometry = geom.graphene() + >>> grid = Grid(0.1, lattice=lattice, geometry=geometry) """ #: Constant for defining a periodic boundary condition @@ -83,10 +91,12 @@ def __init__(self, shape, bc=None, lattice=None, dtype=None, geometry=None): if bc is None: bc = [[self.PERIODIC] * 2] * 3 - self.set_lattice(lattice) + self.set_lattice(None) # Create the atomic structure in the grid, if possible self.set_geometry(geometry) + if lattice is not None: + self.set_lattice(lattice) if isinstance(shape, Real): d = (self.cell ** 2).sum(1) ** 0.5 @@ -98,12 +108,6 @@ def __init__(self, shape, bc=None, lattice=None, dtype=None, geometry=None): # Create the grid boundary conditions self.set_bc(bc) - # If the user sets the super-cell, that has precedence. - if lattice is not None: - if not self.geometry is None: - self.geometry.set_lattice(lattice) - self.set_lattice(lattice) - def __getitem__(self, key): """ Grid value at `key` """ return self.grid[key] @@ -112,6 +116,21 @@ def __setitem__(self, key, val): """ Updates the grid contained """ self.grid[key] = val + def _is_commensurate(self): + """Determine whether the contained geometry and lattice are commensurate """ + if self.geometry is None: + return True + # ideally this should be checked that they are integer equivalent + l_len = self.lattice.length + g_len = self.geometry.lattice.length + reps = np.ones(3) + for i, (l, g) in enumerate(zip(l_len, g_len)): + if l >= g: + reps[i] = l / g + else: + return False + return np.all(abs(reps - np.round(reps)) < 1e-5) + def set_geometry(self, geometry): """ Sets the `Geometry` for the grid. @@ -629,19 +648,26 @@ def tile(self, reps, axis): axis : int direction of tiling, 0, 1, 2 according to the cell-direction + Raises + ------ + SislError : when the lattice is not commensurate with the geometry + See Also -------- Geometry.tile : equivalent method for Geometry class """ + if not self._is_commensurate(): + raise SislError(f"{self.__class__.__name__} cannot tile the grid since the contained" + " Geometry and Lattice are not commensurate.") grid = self.copy() grid.grid = None reps_all = [1, 1, 1] reps_all[axis] = reps grid.grid = np.tile(self.grid, reps_all) - if self.geometry is None: - grid.set_lattice(self.lattice.tile(reps, axis)) - else: + lattice = self.lattice.tile(reps, axis) + if self.geometry is not None: grid.set_geometry(self.geometry.tile(reps, axis)) + grid.set_lattice(lattice) return grid def index2xyz(self, index): @@ -953,11 +979,14 @@ def __str__(self): s += ' bc: [{}, {},\n {}, {},\n {}, {}],\n '.format(bc[self.bc[0, 0]], bc[self.bc[0, 1]], bc[self.bc[1, 0]], bc[self.bc[1, 1]], bc[self.bc[2, 0]], bc[self.bc[2, 1]]) - if not self.geometry is None: - s += '{}\n}}'.format(str(self.geometry).replace('\n', '\n ')) + if self._is_commensurate() and self.geometry is not None: + l = np.round(self.lattice.length / self.geometry.lattice.length).astype(np.int32) + s += f"commensurate: [{l[0]} {l[1]} {l[2]}]" else: - s += '{}\n}}'.format(str(self.lattice).replace('\n', '\n ')) - return s + s += '{}'.format(str(self.lattice).replace('\n', '\n ')) + if not self.geometry is None: + s += ',\n {}'.format(str(self.geometry).replace('\n', '\n ')) + return f"{s}\n}}" def _check_compatibility(self, other, msg): """ Internal check for asserting two grids are commensurable """ diff --git a/src/sisl/lattice.py b/src/sisl/lattice.py index c8f31a7783..df7757c965 100644 --- a/src/sisl/lattice.py +++ b/src/sisl/lattice.py @@ -111,12 +111,13 @@ def find_min_max(cmin, cmax, new): for i in range(3): cmin[i] = min(cmin[i], new[i]) cmax[i] = max(cmax[i], new[i]) - cmin = self.cell.min(0) - cmax = self.cell.max(0) - find_min_max(cmin, cmax, self.cell[[0, 1], :].sum(0)) - find_min_max(cmin, cmax, self.cell[[0, 2], :].sum(0)) - find_min_max(cmin, cmax, self.cell[[1, 2], :].sum(0)) - find_min_max(cmin, cmax, self.cell.sum(0)) + cell = self.cell + cmin = cell.min(0) + cmax = cell.max(0) + find_min_max(cmin, cmax, cell[[0, 1], :].sum(0)) + find_min_max(cmin, cmax, cell[[0, 2], :].sum(0)) + find_min_max(cmin, cmax, cell[[1, 2], :].sum(0)) + find_min_max(cmin, cmax, cell.sum(0)) return Cuboid(cmax - cmin, self.center() + self.origin) def parameters(self, rad=False) -> Tuple[float, float, float, float, float, float]: @@ -1060,7 +1061,7 @@ def read(sile, *args, **kwargs): return fh.read_lattice(*args, **kwargs) def equal(self, other, tol=1e-4): - """ Check whether two supercell are equivalent + """ Check whether two lattices are equivalent Parameters ---------- @@ -1209,20 +1210,6 @@ def set_lattice(self, lattice): # The supercell is given as a cell self.lattice = Lattice(lattice) - # Loop over attributes in this class - # if it inherits LatticeChild, we call - # set_sc on that too. - # Sadly, getattr fails for @property methods - # which forces us to use try ... except - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - for a in dir(self): - try: - if isinstance(getattr(self, a), LatticeChild): - getattr(self, a).set_lattice(self.lattice) - except Exception: - pass - set_sc = deprecation("set_sc is deprecated; please use set_lattice instead", "0.14")(set_lattice) set_supercell = deprecation("set_sc is deprecated; please use set_lattice instead", "0.15")(set_lattice) diff --git a/src/sisl/physics/electron.py b/src/sisl/physics/electron.py index db10a54b5d..5ebb6be7b5 100644 --- a/src/sisl/physics/electron.py +++ b/src/sisl/physics/electron.py @@ -1114,8 +1114,7 @@ def wavefunction(v, grid, geometry=None, k=None, spinor=0, spin=None, eta=None): if geometry is None: raise SislError("wavefunction: did not find a usable Geometry through keywords or the Grid!") # Ensure coordinates are in the primary unit-cell, regardless of origin etc. - geometry = geometry.copy() - geometry.xyz = (geometry.fxyz % 1) @ geometry.lattice.cell + geometry = geometry.translate2uc(axes=(0, 1, 2)) # In case the user has passed several vectors we sum them to plot the summed state if v.ndim == 2: @@ -1159,6 +1158,7 @@ def wavefunction(v, grid, geometry=None, k=None, spinor=0, spin=None, eta=None): # Extract sub variables used throughout the loop shape = _a.asarrayi(grid.shape) dcell = grid.dcell + dlen = (dcell ** 2).sum(1) ** 0.5 ic_shape = grid.lattice.icell * shape.reshape(3, 1) # Convert the geometry (hosting the wavefunction coefficients) coordinates into @@ -1188,6 +1188,8 @@ def idx2spherical(ix, iy, iz, offset, dc, R): return n, idx, rx, ry, rz # Figure out the max-min indices with a spacing of 1 radian + # calculate based on the minimum length of the grid-spacing + dlen_min = dlen.min() rad1 = pi / 180 theta, phi = ogrid[-pi:pi:rad1, 0:pi:rad1] cphi, sphi = cos(phi), sin(phi) @@ -1244,8 +1246,10 @@ def idx2spherical(ix, iy, iz, offset, dc, R): # We can *perhaps* easily attach a geometry with the given # atoms in the unit-cell lattice = grid.lattice.copy() + # Find the periodic directions - pbc = [bc == grid.PERIODIC or geometry.nsc[i] > 1 for i, bc in enumerate(grid.bc[:, 0])] + pbc = [bc == grid.PERIODIC or geometry.nsc[i] > 1 + for i, bc in enumerate(grid.bc[:, 0])] if grid.geometry is None: # Create the actual geometry that encompass the grid ia, xyz, _ = geometry.within_inf(lattice, periodic=pbc) @@ -1260,7 +1264,7 @@ def idx2spherical(ix, iy, iz, offset, dc, R): # For extremely skewed lattices this will be way too much, hence we make # them square. - o = lattice.toCuboid(True) + o = lattice.toCuboid(orthogonal=True) lattice = Lattice(o._v + np.diag(2 * add_R), origin=o.origin - add_R) # Retrieve all atoms within the grid supercell diff --git a/src/sisl/tests/test_grid.py b/src/sisl/tests/test_grid.py index ecadbf51de..a70cb970e0 100644 --- a/src/sisl/tests/test_grid.py +++ b/src/sisl/tests/test_grid.py @@ -7,7 +7,19 @@ import pytest from scipy.sparse import csr_matrix -from sisl import Atom, Cuboid, Ellipsoid, Geometry, Grid, Lattice, SphericalOrbital +from sisl import ( + Atom, + Cuboid, + Ellipsoid, + Geometry, + Grid, + Lattice, + SislError, + SphericalOrbital, + geom, +) + +pytestmark = [pytest.mark.grid] @pytest.fixture @@ -485,3 +497,24 @@ def test_grid_tile_geom(): assert grid.shape[2] == grid4.shape[2] // 2 assert grid.volume * 4 == pytest.approx(grid4.volume) assert grid.geometry.na * 4 == grid4.geometry.na + + +def test_grid_tile_commensurate(): + gr = geom.graphene() + gr_lat = gr.lattice.tile(2, 0).tile(2, 1) + grid = Grid([4, 5, 6], geometry=gr, lattice=gr_lat) + str(grid) + grid2 = grid.tile(2, 2) + assert grid.shape[:2] == grid2.shape[:2] + assert grid.shape[2] == grid2.shape[2] // 2 + assert grid.volume * 2 == pytest.approx(grid2.volume) + assert grid.geometry.na * 2 == grid2.geometry.na + + +def test_grid_tile_in_commensurate(): + gr = geom.graphene() + lat = Lattice(4.) + grid = Grid([4, 5, 6], geometry=gr, lattice=lat) + str(grid) + with pytest.raises(SislError): + grid.tile(2, 2) diff --git a/src/sisl/viz/__init__.py b/src/sisl/viz/__init__.py index 3ae41ed598..39644e8ec2 100644 --- a/src/sisl/viz/__init__.py +++ b/src/sisl/viz/__init__.py @@ -31,6 +31,7 @@ # isort: split from .plot import Animation, MultiplePlot, Plot, SubPlots + # isort: split from ._plotables import register_plotable