lots of bug-fixes that had problems with pbc migration

The boundary conditions were rarely properly
shipped together with the Lattice methods.
Basically it boiled down to the copy method
which did not take all variables into account.

This commit also addes a Listify|listify
class|object that can be used to ensures
arguments are converted to a list.
It can be quite handy as it allows "piping".

Enabled a setter for pbc to fast changing
stuff. One can pass *none* to not change it.

Signed-off-by: Nick Papior <[email protected]>

CHANGELOG.md

@@ -9,6 +9,9 @@ we hit release version 1.0.0.
 ## [0.15.0] - YYYY-MM-DD
 
 ### Added
+- added `Listify` to ensure arguments behaves as *iterables*
+- setter for `Lattice.pbc` to specify it through an array
+- `Lattice.volumef` to calculate a subset volume based on axes
 - added `write_grid` to Siesta binary grid files
 - added the `goldene` 2D lattice, a `hexagonal` Gold 2D structure
 - added the `hexagonal` 2D lattice, close-packed FCC(111) surface
@@ -47,6 +50,9 @@ we hit release version 1.0.0.
 - A new `AtomicMatrixPlot` to plot sparse matrices, #668
 
 ### Fixed
+- lots of `Lattice` methods did not consistently copy over BC
+- `BrillouinZone.volume` fixed to actually return BZ volume
+  use `Lattice.volume` for getting the lattice volume.
 - xsf files now only respect `lattice.pbc` for determining PBC, #764
 - fixed `CHGCAR` spin-polarized density reads, #754
 - dispatch methods now searches the mro for best matches, #721
@@ -85,7 +91,8 @@ we hit release version 1.0.0.
 - removed `Selector` and `TimeSelector`, they were never used internally
 
 ### Changed
-- changed `Eigenmode.displacement` return shape, please read the documentation
+- deprecated `periodic` to `axes` argument in `BrillouinZone.volume`
+- changed `Eigenmode.displacement` shape, please read the documentation
 - bumped minimal Python version to 3.9, #640
 - documentation build system on RTD is updated, #745
 - `gauge` arguments now accept 'cell' and 'orbital' in replacements for 'R' and 'r', respectively

docs/api/utilities.misc.rst

@@ -26,6 +26,8 @@ Miscellaneous routines
 
    str_spec
    direction - abc/012 -> 012
+   Listify
+   listify
    angle - radian to degree
    iter_shape
    math_eval

src/sisl/_core/_ufuncs_lattice.py

@@ -26,31 +26,31 @@
 
 
 @register_sisl_dispatch(Lattice, module="sisl")
-def copy(lattice: Lattice, cell=None, origin: Optional[Coord] = None) -> Lattice:
+def copy(lattice: Lattice, cell=None, **kwargs) -> Lattice:
     """A deepcopy of the object
 
     Parameters
     ----------
     cell : array_like
        the new cell parameters
-    origin : array_like
-       the new origin
     """
     d = dict()
-    d["nsc"] = lattice.nsc.copy()
-    d["boundary_condition"] = lattice.boundary_condition.copy()
-    if origin is None:
-        d["origin"] = lattice.origin.copy()
-    else:
-        d["origin"] = origin
+    for key in ("nsc", "boundary_condition", "origin"):
+        if key in kwargs:
+            d[key] = kwargs.pop(key)
+        else:
+            d[key] = getattr(lattice, key).copy()
     if cell is None:
         d["cell"] = lattice.cell.copy()
     else:
         d["cell"] = np.array(cell)
+    assert len(kwargs) == 0, f"Unknown arguments passed to Lattice.copy {kwargs.keys()}"
 
     copy = lattice.__class__(**d)
     # Ensure that the correct super-cell information gets carried through
-    if not np.allclose(copy.sc_off, lattice.sc_off):
+    if np.allclose(copy.nsc, lattice.nsc) and not np.allclose(
+        copy.sc_off, lattice.sc_off
+    ):
         copy.sc_off = lattice.sc_off
     return copy
 
@@ -195,7 +195,7 @@ def swapaxes(
             origin = origin[idx]
             bc = bc[idx]
 
-    return lattice.__class__(
+    return lattice.copy(
         cell.copy(), nsc=nsc.copy(), origin=origin.copy(), boundary_condition=bc
     )
 
@@ -269,7 +269,7 @@ def add(lattice: Lattice, other) -> Lattice:
     cell = lattice.cell + other.cell
     origin = lattice.origin + other.origin
     nsc = np.where(lattice.nsc > other.nsc, lattice.nsc, other.nsc)
-    return lattice.__class__(cell, nsc=nsc, origin=origin)
+    return lattice.copy(cell, nsc=nsc, origin=origin)
 
 
 @register_sisl_dispatch(Lattice, module="sisl")
@@ -290,15 +290,14 @@ def tile(lattice: Lattice, reps: int, axis: CellAxis) -> Lattice:
     axis = direction(axis)
     cell = np.copy(lattice.cell)
     nsc = np.copy(lattice.nsc)
-    origin = np.copy(lattice.origin)
     cell[axis] *= reps
     # Only reduce the size if it is larger than 5
     if nsc[axis] > 3 and reps > 1:
         # This is number of connections for the primary cell
         h_nsc = nsc[axis] // 2
         # The new number of supercells will then be
         nsc[axis] = max(1, int(math.ceil(h_nsc / reps))) * 2 + 1
-    return lattice.__class__(cell, nsc=nsc, origin=origin)
+    return lattice.copy(cell, nsc=nsc)
 
 
 @register_sisl_dispatch(Lattice, module="sisl")

src/sisl/_core/geometry.py

@@ -49,13 +49,21 @@
 from sisl._namedindex import NamedIndex
 from sisl.messages import SislError, deprecate_argument, info, warn
 from sisl.shape import Cube, Shape, Sphere
-from sisl.typing import ArrayLike, AtomsIndex, NDArray, OrbitalsIndex, SileLike
+from sisl.typing import (
+    ArrayLike,
+    AtomsIndex,
+    CellAxes,
+    NDArray,
+    OrbitalsIndex,
+    SileLike,
+)
 from sisl.utils import (
     angle,
     cmd,
     default_ArgumentParser,
     default_namespace,
     direction,
+    listify,
     lstranges,
     str_spec,
     strmap,
@@ -1746,6 +1754,7 @@ def translate2uc(
                 raise ValueError(
                     "translate2uc with a bool argument can only be True to signal all axes"
                 )
+        axes = map(direction, listify(axes)) | listify
 
         fxyz = self.fxyz
         # move to unit-cell
@@ -3358,18 +3367,19 @@ def func(lst):
     def within_inf(
         self,
         lattice: Lattice,
-        periodic: Optional[Sequence[bool]] = None,
+        periodic: Optional[Union[Sequence[bool], CellAxes]] = None,
         tol: float = 1e-5,
         origin: Sequence[float] = (0.0, 0.0, 0.0),
     ) -> Tuple[ndarray, ndarray, ndarray]:
         """Find all atoms within a provided supercell
 
         Note this function is rather different from `close` and `within`.
         Specifically this routine is returning *all* indices for the infinite
-        periodic system (where ``self.pbc`` or `periodic` is true).
+        periodic system. The default periodic directions are ``self.pbc``,
+        unless `periodic` is provided.
 
         Atomic coordinates lying on the boundary of the supercell will be duplicated
-        on the neighboring supercell images. Thus performing `geom.within_inf(geom.lattice)`
+        on the neighboring supercell images. Thus performing ``geom.within_inf(geom.lattice)``
         may result in more atoms than in the structure.
 
         Notes
@@ -3402,15 +3412,25 @@ def within_inf(
         """
         lattice = Lattice.new(lattice)
         if periodic is None:
-            periodic = self.pbc
+            periodic = self.pbc.nonzero()[0]
+        elif isinstance(periodic, bool):
+            periodic = (0, 1, 2)
         else:
-            periodic = list(periodic)
+            try:
+                periodic = map(direction, listify(periodic)) | listify
+            except:
+                periodic = np.asarray(periodic).nonzero()[0]
+
+        # extract the non-periodic directions
+        non_periodic = filter(lambda i: i not in periodic, range(3)) | listify
 
         if origin is None:
             origin = _a.zerosd(3)
 
         # Our first task is to construct a geometry large
         # enough to fully encompass the supercell
+        # The supercell here defines how big `self` needs to be
+        # to be fully located inside `lattice`.
 
         # 1. Number of times each lattice vector must be expanded to fit
         #    inside the "possibly" larger `lattice`.
@@ -3419,7 +3439,7 @@ def within_inf(
         tile_max = ceil(idx.max(0)).astype(dtype=int32)
 
         # Intrinsic offset (when atomic coordinates are outside primary unit-cell)
-        idx = dot(self.xyz, self.icell.T)
+        idx = self.fxyz
         tmp = floor(idx.min(0))
         tile_min = np.where(tile_min < tmp, tile_min, tmp).astype(dtype=int32)
         tmp = ceil(idx.max(0))
@@ -3433,8 +3453,8 @@ def within_inf(
         tile_min = np.where(tile_min < idx, tile_min, idx).astype(dtype=int32)
 
         # 2. Reduce tiling along non-periodic directions
-        tile_min = np.where(periodic, tile_min, 0)
-        tile_max = np.where(periodic, tile_max, 1)
+        tile_min[non_periodic] = 0
+        tile_max[non_periodic] = 1
 
         # 3. Find the *new* origin according to the *negative* tilings.
         #    This is important for skewed cells as the placement of the new
@@ -3443,14 +3463,17 @@ def within_inf(
 
         # The xyz geometry that fully encompass the (possibly) larger supercell
         tile = tile_max - tile_min
+
         full_geom = (self * tile).translate(big_origin - origin)
 
         # Now we have to figure out all atomic coordinates within
         cuboid = lattice.to.Cuboid()
 
         # Make sure that full_geom doesn't return coordinates outside the unit cell
         # for non periodic directions
-        full_geom.set_nsc([full_geom.nsc[i] if periodic[i] else 1 for i in range(3)])
+        nsc = full_geom.nsc.copy()
+        nsc[non_periodic] = 1
+        full_geom.set_nsc(nsc)
 
         # Now retrieve all atomic coordinates from the full geometry
         xyz = full_geom.axyz(_a.arangei(full_geom.na_s))

src/sisl/_core/lattice.py

@@ -26,7 +26,7 @@
 from sisl.shape.prism4 import Cuboid
 from sisl.typing import CellAxes, CellAxis
 from sisl.utils.mathematics import fnorm
-from sisl.utils.misc import direction
+from sisl.utils.misc import direction, listify
 
 from ._lattice import cell_invert, cell_reciprocal
 
@@ -151,16 +151,64 @@ def length(self) -> ndarray:
         """Length of each lattice vector"""
         return fnorm(self.cell)
 
+    def lengthf(self, axes: CellAxes = (0, 1, 2)) -> ndarray:
+        """Length of specific lattice vectors (as a function)
+        Parameters
+        ----------
+        axes:
+            only calculate the volume based on a subset of axes
+
+        Examples
+        --------
+        Only get lengths of two lattice vectors:
+
+        >>> lat = Lattice(1)
+        >>> lat.lengthf([0, 1])
+        """
+        axes = map(direction, listify(axes)) | listify
+        return fnorm(self.cell[axes])
+
     @property
     def volume(self) -> float:
         """Volume of cell"""
-        return abs(dot3(self.cell[0], cross3(self.cell[1], self.cell[2])))
+        return self.volumef((0, 1, 2))
+
+    def volumef(self, axes: CellAxes = (0, 1, 2)) -> float:
+        """Volume of cell (as a function)
+
+        Default to the 3D volume.
+        For `axes` with only 2 elements, it corresponds to an area.
+        For `axes` with only 1 element, it corresponds to a length.
+
+        Parameters
+        ----------
+        axes:
+            only calculate the volume based on a subset of axes
+
+        Examples
+        --------
+        Only get the volume of the periodic directions:
+
+        >>> lat = Lattice(1)
+        >>> lat.pbc = (True, False, True)
+        >>> lat.volumef(lat.pbc.nonzero()[0])
+        """
+        axes = map(direction, listify(axes)) | listify
+
+        cell = self.cell
+        if len(axes) == 3:
+            return abs(dot3(cell[axes[0]], cross3(cell[axes[1]], cell[axes[2]])))
+        if len(axes) == 2:
+            return fnorm(cross3(cell[axes[0]], cell[axes[1]]))
+        if len(axes) == 1:
+            return fnorm(cell[axes])
+        return 0.0
 
     def area(self, axis1: CellAxis, axis2: CellAxis) -> float:
         """Calculate the area spanned by the two axis `ax0` and `ax1`"""
         axis1 = direction(axis1)
         axis2 = direction(axis2)
-        return (cross3(self.cell[axis1], self.cell[axis2]) ** 2).sum() ** 0.5
+        return fnorm(cross3(self.cell[axis1], self.cell[axis2]))
 
     @property
     def boundary_condition(self) -> np.ndarray:
@@ -185,9 +233,18 @@ def pbc(self, pbc) -> None:
         # set_boundary_condition does not allow to have PERIODIC and non-PERIODIC
         # along the same lattice vector. So checking one should suffice
         assert len(pbc) == 3
+
+        PERIODIC = BoundaryCondition.PERIODIC
         for axis, bc in enumerate(pbc):
+
+            # Simply skip those that are not T|F
+            if not isinstance(bc, bool):
+                continue
+
             if bc:
-                self._bc[axis] = BoundaryCondition.PERIODIC
+                self._bc[axis] = PERIODIC
+            elif self._bc[axis, 0] == PERIODIC:
+                self._bc[axis] = BoundaryCondition.UNKNOWN
 
     @property
     def origin(self) -> ndarray:
@@ -532,13 +589,9 @@ def fit(self, xyz, axes: CellAxes = (0, 1, 2), tol: float = 0.05) -> Lattice:
         # Reduce to total repetitions
         ireps = np.amax(ix, axis=0) - np.amin(ix, axis=0) + 1
 
-        # Only repeat the axis requested
-        if isinstance(axes, Integral):
-            axes = [axes]
-
         # Reduce the non-set axis
         if not axes is None:
-            axes = list(map(direction, axes))
+            axes = map(direction, listify(axes))
             for ax in (0, 1, 2):
                 if ax not in axes:
                     ireps[ax] = 1
@@ -606,6 +659,7 @@ def plane(
         """
         axis1 = direction(axis1)
         axis2 = direction(axis2)
+
         cell = self.cell
         n = cross3(cell[axis1], cell[axis2])
         # Normalize
@@ -683,9 +737,7 @@ def cell2length(self, length, axes: CellAxes = (0, 1, 2)) -> ndarray:
         numpy.ndarray
              cell-vectors with prescribed length, same order as `axes`
         """
-        if isinstance(axes, Integral):
-            axes = [axes]
-        axes = list(map(direction, axes))
+        axes = map(direction, listify(axes)) | listify
 
         length = _a.asarray(length).ravel()
         if len(length) != len(axes):
@@ -696,6 +748,7 @@ def cell2length(self, length, axes: CellAxes = (0, 1, 2)) -> ndarray:
                     f"{self.__class__.__name__}.cell2length length parameter should be a single "
                     "float, or an array of values according to axes argument."
                 )
+
         return self.cell[axes] * (length / self.length[axes]).reshape(-1, 1)
 
     def offset(self, isc=None) -> Tuple[float, float, float]:
@@ -942,9 +995,8 @@ def parallel(self, other, axes: CellAxes = (0, 1, 2)) -> bool:
         axes :
            only check the specified axes (default to all)
         """
-        if isinstance(axes, Integral):
-            axes = [axes]
-        axis = list(map(direction, axes))
+        axis = map(direction, listify(axes))
+
         # Convert to unit-vector cell
         for i in axis:
             a = self.cell[i] / fnorm(self.cell[i])