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 <[email protected]>

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

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 """

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)
 

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

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)

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