From bb4cb6776c53ba9fd013af1b3d2c6eeb7538aba0 Mon Sep 17 00:00:00 2001 From: Nick Papior Date: Fri, 20 Oct 2023 19:16:29 +0200 Subject: [PATCH] Lattice is now owning the boundary conditions This should be a first step towards moving the nsc attribute from Lattice to somewhere more appropriate. This also makes the boundary conditions in the Grid more manageable since it is the lattice which keeps the information. Simplified reading xyz files since the boundary conditions are now explicit. Signed-off-by: Nick Papior --- CHANGELOG.md | 3 + docs/api/basic.rst | 1 + src/sisl/grid.py | 136 +++++++-------- src/sisl/io/tests/test_xyz.py | 3 +- src/sisl/io/xyz.py | 114 ++++++------- src/sisl/lattice.py | 161 ++++++++++++++++-- src/sisl/physics/densitymatrix.py | 4 +- src/sisl/physics/electron.py | 5 +- src/sisl/tests/test_grid.py | 24 +-- src/sisl/tests/test_lattice.py | 23 ++- .../transiesta/poisson/fftpoisson_fix.py | 8 +- 11 files changed, 303 insertions(+), 179 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1c4c5ff03a..f28171fda3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,9 @@ we hit release version 1.0.0. debugging. - marked all `toSphere|toEllipsoid|...` as deprecated +### Changed +- `Lattice` now holds the boundary conditions (not `Grid`), see #626 + ## [0.14.2] - 2023-10-04 diff --git a/docs/api/basic.rst b/docs/api/basic.rst index 4823d461f9..25a3b78f48 100644 --- a/docs/api/basic.rst +++ b/docs/api/basic.rst @@ -30,6 +30,7 @@ Simple objects Atoms Geometry Lattice + BoundaryCondition Grid diff --git a/src/sisl/grid.py b/src/sisl/grid.py index 223de3bfb3..662620b45d 100644 --- a/src/sisl/grid.py +++ b/src/sisl/grid.py @@ -1,6 +1,7 @@ # This Source Code Form is subject to the terms of the Mozilla Public # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at https://mozilla.org/MPL/2.0/. +import logging from math import pi from numbers import Integral, Real @@ -14,8 +15,8 @@ from ._help import dtype_complex_to_real, wrap_filterwarnings from ._internal import set_module from .geometry import Geometry -from .lattice import LatticeChild -from .messages import SislError, deprecate_argument +from .lattice import BoundaryCondition, LatticeChild +from .messages import SislError, deprecate_argument, deprecation from .shape import Shape from .utils import ( cmd, @@ -30,6 +31,10 @@ __all__ = ['Grid', 'sgrid'] +_log = logging.getLogger("sisl") +_log.info(f"adding logger: {__name__}") +_log = logging.getLogger(__name__) + @set_module("sisl") class Grid(LatticeChild): @@ -76,26 +81,29 @@ class Grid(LatticeChild): """ #: Constant for defining a periodic boundary condition - PERIODIC = 1 + PERIODIC = BoundaryCondition.PERIODIC #: Constant for defining a Neumann boundary condition - NEUMANN = 2 + NEUMANN = BoundaryCondition.NEUMANN #: Constant for defining a Dirichlet boundary condition - DIRICHLET = 3 + DIRICHLET = BoundaryCondition.DIRICHLET #: Constant for defining an open boundary condition - OPEN = 4 + OPEN = BoundaryCondition.OPEN @deprecate_argument("sc", "lattice", "argument sc has been deprecated in favor of lattice, please update your code.", "0.15.0") + @deprecate_argument("bc", None, + "argument bc has been deprecated (removed) in favor of the boundary conditions in Lattice, please update your code.", + "0.15.0") def __init__(self, shape, bc=None, lattice=None, dtype=None, geometry=None): - if bc is None: - bc = [[self.PERIODIC] * 2] * 3 self.set_lattice(None) # Create the atomic structure in the grid, if possible self.set_geometry(geometry) if lattice is not None: + if bc is None: + bc = [[self.PERIODIC] * 2] * 3 self.set_lattice(lattice) if isinstance(shape, Real): @@ -106,7 +114,20 @@ def __init__(self, shape, bc=None, lattice=None, dtype=None, geometry=None): self.set_grid(shape, dtype=dtype) # Create the grid boundary conditions - self.set_bc(bc) + if bc is not None: + self.lattice.set_boundary_condition(bc) + + @deprecation("Grid.set_bc is deprecated since boundary conditions are moved to Lattice (see github issue #626", "0.15.0") + def set_bc(self, bc): + self.lattice.set_boundary_condition(bc) + + @deprecation("Grid.set_boundary is deprecated since boundary conditions are moved to Lattice (see github issue #626", "0.15.0") + def set_boundary(self, bc): + self.lattice.set_boundary_condition(bc) + + @deprecation("Grid.set_boundary_condition is deprecated since boundary conditions are moved to Lattice (see github issue #626", "0.15.0") + def set_boundary_condition(self, bc): + self.lattice.set_boundary_condition(bc) def __getitem__(self, key): """ Grid value at `key` """ @@ -131,20 +152,31 @@ def _is_commensurate(self): return False return np.all(abs(reps - np.round(reps)) < 1e-5) - def set_geometry(self, geometry): + def set_geometry(self, geometry, also_lattice: bool=True): """ Sets the `Geometry` for the grid. Setting the `Geometry` for the grid is a possibility to attach atoms to the grid. It is not a necessary entity, so passing `None` is a viable option. + + Parameters + ---------- + geometry : Geometry or None + specify the new geometry in the `Grid`. If ``None`` will + remove the geometry (but not the lattice) + also_lattice : bool, optional + whether to also set the lattice for the grid according to the + lattice of the `geometry`, if ``False`` it will keep the lattice + as it was. """ if geometry is None: # Fake geometry self.geometry = None else: self.geometry = geometry - self.set_lattice(geometry.lattice) + if also_lattice: + self.set_lattice(geometry.lattice) def fill(self, val): """ Fill the grid with this value @@ -350,51 +382,7 @@ def set_grid(self, shape, dtype=None): raise ValueError(f"{self.__class__.__name__}.set_grid requires shape to be of length 3") self.grid = np.zeros(shape, dtype=dtype) - def set_bc(self, boundary=None, a=None, b=None, c=None): - """ Set the boundary conditions on the grid - - Parameters - ---------- - boundary : (3, 2) or (3, ) or int, optional - boundary condition for all boundaries (or the same for all) - a : int or list of int, optional - boundary condition for the first unit-cell vector direction - b : int or list of int, optional - boundary condition for the second unit-cell vector direction - c : int or list of int, optional - boundary condition for the third unit-cell vector direction - - Raises - ------ - ValueError - if specifying periodic one one boundary, so must the opposite side. - """ - if not boundary is None: - if isinstance(boundary, Integral): - self.bc = _a.fulli([3, 2], boundary) - else: - self.bc = _a.asarrayi(boundary) - if not a is None: - self.bc[0, :] = a - if not b is None: - self.bc[1, :] = b - if not c is None: - self.bc[2, :] = c - - # shorthand for bc - bc = self.bc[:, :] - for i in range(3): - if (bc[i, 0] == self.PERIODIC and bc[i, 1] != self.PERIODIC) or \ - (bc[i, 0] != self.PERIODIC and bc[i, 1] == self.PERIODIC): - raise ValueError(f"{self.__class__.__name__}.set_bc has a one non-periodic and " - "one periodic direction. If one direction is periodic, both instances " - "must have that BC.") - - # Aliases - set_boundary = set_bc - set_boundary_condition = set_bc - - def __sc_geometry_dict(self): + def _sc_geometry_dict(self): """ Internal routine for copying the Lattice and Geometry """ d = dict() d['lattice'] = self.lattice.copy() @@ -404,12 +392,12 @@ def __sc_geometry_dict(self): def copy(self, dtype=None): r""" Copy the object, possibly changing the data-type """ - d = self.__sc_geometry_dict() + d = self._sc_geometry_dict() if dtype is None: d['dtype'] = self.dtype else: d['dtype'] = dtype - grid = self.__class__([1] * 3, bc=np.copy(self.bc), **d) + grid = self.__class__([1] * 3, **d) # This also ensures the shape is copied! grid.grid = self.grid.astype(dtype=d['dtype']) return grid @@ -429,10 +417,10 @@ def swapaxes(self, a, b): idx[b] = a idx[a] = b s = np.copy(self.shape) - d = self.__sc_geometry_dict() + d = self._sc_geometry_dict() d['lattice'] = d['lattice'].swapaxes(a, b) d['dtype'] = self.dtype - grid = self.__class__(s[idx], bc=self.bc[idx], **d) + grid = self.__class__(s[idx], **d) # We need to force the C-order or we loose the contiguity grid.grid = np.copy(np.swapaxes(self.grid, a, b), order='C') return grid @@ -457,7 +445,7 @@ def _copy_sub(self, n, axis, scale_geometry=False): shape[axis] = n if n < 1: raise ValueError('You cannot retain no indices.') - grid = self.__class__(shape, bc=np.copy(self.bc), dtype=self.dtype, **self.__sc_geometry_dict()) + grid = self.__class__(shape, dtype=self.dtype, **self._sc_geometry_dict()) # Update cell shape (the cell is smaller now) grid.set_lattice(cell) if scale_geometry and not self.geometry is None: @@ -913,14 +901,15 @@ def append(self, other, axis): """ Appends other `Grid` to this grid along axis """ shape = list(self.shape) shape[axis] += other.shape[axis] - d = self.__sc_geometry_dict() + d = self._sc_geometry_dict() if 'geometry' in d: if not other.geometry is None: d['geometry'] = d['geometry'].append(other.geometry, axis) else: d['geometry'] = other.geometry + d['lattice'] = self.lattice.append(other.lattice, axis) d['dtype'] = self.dtype - return self.__class__(shape, bc=np.copy(self.bc), **d) + return self.__class__(shape, **d) @staticmethod def read(sile, *args, **kwargs): @@ -970,15 +959,7 @@ def write(self, sile, *args, **kwargs) -> None: def __str__(self): """ String of object """ - s = self.__class__.__name__ + '{{kind: {kind}, shape: [{shape[0]} {shape[1]} {shape[2]}],\n'.format(kind=self.dkind, shape=self.shape) - bc = {self.PERIODIC: 'periodic', - self.NEUMANN: 'neumann', - self.DIRICHLET: 'dirichlet', - self.OPEN: 'open' - } - 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]]) + s = '{name}{{kind: {kind}, shape: [{shape[0]} {shape[1]} {shape[2]}],\n'.format(kind=self.dkind, shape=self.shape, name=self.__class__.__name__) 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]}]" @@ -1280,13 +1261,14 @@ def sl2idx(sl): new_sl = sl[:] # LOWER BOUNDARY - bc = self.bc[i, 0] + bci = self.boundary_condition[i] new_sl[i] = slice(0, 1) idx1 = sl2idx(new_sl) # lower - if self.bc[i, 0] == self.PERIODIC or \ - self.bc[i, 1] == self.PERIODIC: - if self.bc[i, 0] != self.bc[i, 1]: + bc = bci[0] + if bci[0] == self.PERIODIC or \ + bci[1] == self.PERIODIC: + if bci[0] != bci[1]: raise ValueError(f"{self.__class__.__name__}.pyamg_boundary_condition found a periodic and non-periodic direction in the same direction!") new_sl[i] = slice(self.shape[i]-1, self.shape[i]) idx2 = sl2idx(new_sl) # upper @@ -1304,7 +1286,7 @@ def sl2idx(sl): Dirichlet(idx1) # UPPER BOUNDARY - bc = self.bc[i, 1] + bc = bci[1] new_sl[i] = slice(self.shape[i]-1, self.shape[i]) idx1 = sl2idx(new_sl) # upper diff --git a/src/sisl/io/tests/test_xyz.py b/src/sisl/io/tests/test_xyz.py index c3e576214c..fa98fe91e1 100644 --- a/src/sisl/io/tests/test_xyz.py +++ b/src/sisl/io/tests/test_xyz.py @@ -64,7 +64,8 @@ def test_xyz_ase(sisl_tmp): assert np.allclose(g.xyz[:, 0], [0, 1, 2]) assert np.allclose(g.xyz[:, 1], 0.) assert np.allclose(g.xyz[:, 2], 0.) - assert np.allclose(g.nsc, [1, 1, 3]) + assert np.allclose(g.nsc, [1, 1, 1]) + assert np.allclose(g.pbc, [False, False, True]) def test_xyz_arbitrary(sisl_tmp): diff --git a/src/sisl/io/xyz.py b/src/sisl/io/xyz.py index c08f4197c4..0c276399bf 100644 --- a/src/sisl/io/xyz.py +++ b/src/sisl/io/xyz.py @@ -7,7 +7,7 @@ import numpy as np import sisl._array as _a -from sisl import Geometry, Lattice +from sisl import BoundaryCondition, Geometry, Lattice from sisl._internal import set_module from sisl.messages import deprecate_argument, warn @@ -23,6 +23,46 @@ class xyzSile(Sile): """ XYZ file object """ + def _parse_lattice(self, header, xyz, lattice): + """Internal helper routine for extracting the lattice """ + if lattice is not None: + return lattice + + # Parse + nsc = None + if "nsc" in header: + nsc = list(map(int, header.pop("nsc").split())) + + BC = BoundaryCondition + bc = BC.UNKNOWN + if "pbc" in header: + bc = [] + for pbc in header.pop("pbc").split(): + if pbc == "T": + bc.append(BoundaryCondition.PERIODIC) + else: + bc.append(BoundaryCondition.UNKNOWN) + if "boundary_condition" in header: + bc = [] + for b in header.pop("boundary_condition").split(): + bc.append(getattr(BC, b.upper())) + bc = _a.arrayi(bc).reshape(3, 2) + + if "Lattice" in header: + cell = _a.fromiterd(header.pop("Lattice").split()).reshape(3, 3) + elif "cell" in header: + cell = _a.fromiterd(header.pop("cell").split()).reshape(3, 3) + else: + cell = xyz.max(0) - xyz.min(0) + 10 + + origin = None + if "Origin" in header: + origin = _a.fromiterd(header.pop("Origin").strip('"').split()).reshape(3) + + return Lattice(cell, nsc=nsc, + origin=origin, + boundary_condition=bc) + @sile_fh_open() def write_geometry(self, geometry, fmt='.8f', comment=None): """ Writes the geometry to the contained file @@ -39,6 +79,7 @@ def write_geometry(self, geometry, fmt='.8f', comment=None): """ # Check that we can write to the file sile_raise_write(self) + lattice = geometry.lattice # Write the number of atoms in the geometry self._write(' {}\n'.format(len(geometry))) @@ -50,60 +91,24 @@ def write_geometry(self, geometry, fmt='.8f', comment=None): fields.append(('Lattice="' + f'{{:{fmt}}} ' * 9 + '"').format(*geometry.cell.ravel())) nsc = geometry.nsc[:] fields.append('nsc="{} {} {}"'.format(*nsc)) - pbc = ['T' if n else 'F' for n in nsc] + pbc = ['T' if n else 'F' for n in lattice.pbc] fields.append('pbc="{} {} {}"'.format(*pbc)) + BC = BoundaryCondition.getitem + bc = [f"{BC(n[0]).name} {BC(n[1]).name}" for n in lattice.boundary_condition] + fields.append('boundary_condition="{} {} {}"'.format(*bc)) + if comment is not None: + fields.append(f'Comment="{comment}"') - if comment is None: - self._write(' '.join(fields) + "\n") - else: - self._write(' '.join(fields) + f'Comment="{comment}"\n') + self._write(' '.join(fields) + "\n") fmt_str = '{{0:2s}} {{1:{0}}} {{2:{0}}} {{3:{0}}}\n'.format(fmt) for ia, a, _ in geometry.iter_species(): - s = {'fa': 'Ds'}.get(a.symbol, a.symbol) + s = a.symbol + s = {'fa': 'Ds'}.get(s, s) self._write(fmt_str.format(s, *geometry.xyz[ia, :])) - def _r_geometry_sisl(self, na, header, sp, xyz, lattice): - """ Read the geometry as though it was created with sisl """ - # Default version of the header is 1 - #v = int(header.get("sisl-version", 1)) - nsc = list(map(int, header.pop("nsc").split())) - cell = _a.fromiterd(header.pop("cell").split()).reshape(3, 3) - if lattice is None: - lattice = Lattice(cell, nsc=nsc) - return Geometry(xyz, atoms=sp, lattice=lattice) - - def _r_geometry_ase(self, na, header, sp, xyz, lattice): - """ Read the geometry as though it was created with ASE """ - # Convert F T to nsc - # F = 1 - # T = 3 - nsc = header.pop("nsc", "").strip('"') - if nsc: - nsc = list(map(int, nsc.split())) - else: - nsc = list(map(lambda x: "FT".index(x) * 2 + 1, header.pop("pbc").strip('"').split())) - cell = _a.fromiterd(header.pop("Lattice").strip('"').split()).reshape(3, 3) - if "Origin" in header: - origin = _a.fromiterd(header.pop("Origin").strip('"').split()).reshape(3) - else: - origin = None - if lattice is None: - lattice = Lattice(cell, nsc=nsc, origin=origin) - - return Geometry(xyz, atoms=sp, lattice=lattice) - - def _r_geometry(self, na, sp, xyz, lattice): - """ Read the geometry for a generic xyz file (not sisl, nor ASE) """ - # The cell dimensions isn't defined, we are going to create a molecule box - cell = xyz.max(0) - xyz.min(0) + 10. - if lattice is None: - lattice = Lattice(cell, nsc=[1] * 3) - return Geometry(xyz, atoms=sp, lattice=lattice) - def _r_geometry_skip(self, *args, **kwargs): """ Read the geometry for a generic xyz file (not sisl, nor ASE) """ - # The cell dimensions isn't defined, we are going to create a molecule box line = self.readline() if line == '': return None @@ -136,7 +141,9 @@ def read_geometry(self, atoms=None, lattice=None): # Read header, and try and convert to dictionary header = self.readline() - kv = header_to_dict(header) + header = {k: v.strip('"') for k, v in + header_to_dict(header).items() + } # Read atoms and coordinates sp = [None] * na @@ -150,17 +157,8 @@ def read_geometry(self, atoms=None, lattice=None): if atoms is not None: sp = atoms - def _has_keys(d, *keys): - for key in keys: - if not key in d: - return False - return True - - if _has_keys(kv, "cell", "nsc"): - return self._r_geometry_sisl(na, kv, sp, xyz, lattice) - elif _has_keys(kv, "Lattice", "pbc"): - return self._r_geometry_ase(na, kv, sp, xyz, lattice) - return self._r_geometry(na, sp, xyz, lattice) + lattice = self._parse_lattice(header, xyz, lattice) + return Geometry(xyz, atoms=sp, lattice=lattice) def ArgumentParser(self, p=None, *args, **kwargs): """ Returns the arguments that is available for this Sile """ diff --git a/src/sisl/lattice.py b/src/sisl/lattice.py index adb2ffc88a..23c0aa69cf 100644 --- a/src/sisl/lattice.py +++ b/src/sisl/lattice.py @@ -10,9 +10,10 @@ import logging import math import warnings +from enum import IntEnum, auto from numbers import Integral from pathlib import Path -from typing import TYPE_CHECKING, Tuple, Union +from typing import TYPE_CHECKING, Optional, Sequence, Tuple, Union import numpy as np from numpy import dot, ndarray @@ -24,18 +25,36 @@ from ._internal import set_module from ._lattice import cell_invert, cell_reciprocal from ._math_small import cross3, dot3 -from .messages import deprecate, deprecate_argument, deprecation +from .messages import SislError, deprecate, deprecate_argument, deprecation from .quaternion import Quaternion from .shape.prism4 import Cuboid from .utils.mathematics import fnorm -__all__ = ["Lattice", "SuperCell", "LatticeChild"] +__all__ = ["Lattice", "SuperCell", "LatticeChild", "BoundaryCondition"] _log = logging.getLogger("sisl") _log.info(f"adding logger: {__name__}") _log = logging.getLogger(__name__) +class BoundaryCondition(IntEnum): + UNKNOWN = auto() + PERIODIC = auto() + DIRICHLET = auto() + NEUMANN = auto() + OPEN = auto() + + @classmethod + def getitem(cls, key): + """Search for a specific integer entry by value, and not by name """ + for bc in cls: + if bc == key: + return bc + raise KeyError(f"{cls.__name__}.getitem could not find key={key}") + +BoundaryConditionType = Union[int, Sequence[int]] + + @set_module("sisl") class Lattice(_Dispatchs, dispatchs=[ @@ -64,12 +83,19 @@ class Lattice(_Dispatchs, number of supercells along each lattice vector origin : (3,) of float, optional the origin of the supercell. + boundary_condition : list of int (3, 2) or (3, ), optional + the boundary conditions for each of the cell's planes. Defaults to periodic boundary condition. + See `BoundaryCondition` for valid enumerations. """ # We limit the scope of this Lattice object. - __slots__ = ('cell', '_origin', 'nsc', 'n_s', '_sc_off', '_isc_off') + __slots__ = ('cell', '_origin', 'nsc', 'n_s', '_sc_off', '_isc_off', '_bc') + + #: Internal reference to `BoundaryCondition` for simpler short-hands + BC = BoundaryCondition - def __init__(self, cell, nsc=None, origin=None): + def __init__(self, cell, nsc=None, origin=None, + boundary_condition: Sequence[BoundaryConditionType] =BoundaryCondition.PERIODIC): if nsc is None: nsc = [1, 1, 1] @@ -88,6 +114,7 @@ def __init__(self, cell, nsc=None, origin=None): self.nsc = _a.onesi(3) # Set the super-cell self.set_nsc(nsc=nsc) + self.set_boundary_condition(boundary_condition) @property def length(self) -> ndarray: @@ -103,6 +130,22 @@ def area(self, ax0, ax1): """ Calculate the area spanned by the two axis `ax0` and `ax1` """ return (cross3(self.cell[ax0, :], self.cell[ax1, :]) ** 2).sum() ** 0.5 + @property + def boundary_condition(self) -> np.ndarray: + """Boundary conditions for each lattice vector (lower/upper) sides ``(3, 2)``""" + return self._bc + + @boundary_condition.setter + def boundary_condition(self, boundary_condition): + self.set_boundary_condition(boundary_condition) + + @property + def pbc(self) -> np.ndarray: + """Boolean array to specify whether the boundary conditions are periodic`""" + # set_boundary_condition does not allow to have PERIODIC and non-PERIODIC + # along the same lattice vector. So checking one should suffice + return self._bc[:, 0] == BoundaryCondition.PERIODIC + @property def origin(self) -> ndarray: """ Origin for the cell """ @@ -113,7 +156,7 @@ def origin(self, origin): """ Set origin for the cell """ self._origin[:] = origin - @deprecation("toCuboid is deprecated, please use lattice.to['cuboid'](...) instead.") + @deprecation("toCuboid is deprecated, please use lattice.to['cuboid'](...) instead.", "0.15.0") def toCuboid(self, *args, **kwargs): """ A cuboid with vectors as this unit-cell and center with respect to its origin @@ -124,7 +167,51 @@ def toCuboid(self, *args, **kwargs): """ return self.to[Cuboid](*args, **kwargs) - def parameters(self, rad=False) -> Tuple[float, float, float, float, float, float]: + def set_boundary_condition(self, + boundary: Optional[Sequence[BoundaryConditionType]] =None, + a: Sequence[BoundaryConditionType] =None, + b: Sequence[BoundaryConditionType] =None, + c: Sequence[BoundaryConditionType] =None): + """ Set the boundary conditions on the grid + + Parameters + ---------- + boundary : (3, 2) or (3, ) or int, optional + boundary condition for all boundaries (or the same for all) + a : int or list of int, optional + boundary condition for the first unit-cell vector direction + b : int or list of int, optional + boundary condition for the second unit-cell vector direction + c : int or list of int, optional + boundary condition for the third unit-cell vector direction + + Raises + ------ + ValueError + if specifying periodic one one boundary, so must the opposite side. + """ + if not boundary is None: + self._bc = _a.emptyi([3, 2]) + if isinstance(boundary, Integral): + self._bc[:, :] = boundary + else: + for i, bc in enumerate(boundary): + self._bc[i] = bc + if not a is None: + self._bc[0, :] = a + if not b is None: + self._bc[1, :] = b + if not c is None: + self._bc[2, :] = c + + # shorthand for bc + for bc in self._bc == BoundaryCondition.PERIODIC: + if bc.any() and not bc.all(): + raise ValueError(f"{self.__class__.__name__}.set_boundary_condition has a one non-periodic and " + "one periodic direction. If one direction is periodic, both instances " + "must have that BC.") + + def parameters(self, rad: bool=False) -> Tuple[float, float, float, float, float, float]: r""" Cell parameters of this cell in 3 lengths and 3 angles Notes @@ -304,12 +391,19 @@ def copy(self, cell=None, origin=None): origin : array_like the new origin """ + d = dict() + d["nsc"] = self.nsc.copy() + d["boundary_condition"] = self.boundary_condition.copy() if origin is None: - origin = self.origin.copy() + d["origin"] = self.origin.copy() + else: + d["origin"] = origin if cell is None: - copy = self.__class__(np.copy(self.cell), nsc=np.copy(self.nsc), origin=origin) + d["cell"] = self.cell.copy() else: - copy = self.__class__(np.copy(cell), nsc=np.copy(self.nsc), origin=origin) + d["cell"] = np.array(cell) + + copy = self.__class__(**d) # Ensure that the correct super-cell information gets carried through if not np.allclose(copy.sc_off, self.sc_off): copy.sc_off = self.sc_off @@ -454,6 +548,7 @@ def swapaxes(self, axes_a: Union[int, str], cell = self.cell nsc = self.nsc origin = self.origin + bc = self.boundary_condition if len(axes_a) != len(axes_b): raise ValueError(f"{self.__class__.__name__}.swapaxes expects axes_a and axes_b to have the same lengeth {len(axes_a)}, {len(axes_b)}.") @@ -472,6 +567,7 @@ def swapaxes(self, axes_a: Union[int, str], # we are dealing with lattice vectors cell = cell[idx] nsc = nsc[idx] + bc = bc[idx] else: aidx -= 3 @@ -481,8 +577,12 @@ def swapaxes(self, axes_a: Union[int, str], # we are dealing with cartesian coordinates cell = cell[:, idx] origin = origin[idx] + bc = bc[idx] - return self.__class__(cell.copy(), nsc=nsc.copy(), origin=origin.copy()) + return self.__class__(cell.copy(), + nsc=nsc.copy(), + origin=origin.copy(), + boundary_condition=bc) def plane(self, ax1, ax2, origin=True): """ Query point and plane-normal for the plane spanning `ax1` and `ax2` @@ -889,7 +989,7 @@ def prepend(self, other, axis): For a `Lattice` object this is equivalent to `append`. """ - return self.append(other, axis) + return other.append(self, axis) def translate(self, v): """ Appends additional space to the object """ @@ -1098,13 +1198,31 @@ def equal(self, other, tol=1e-4): def __str__(self): """ Returns a string representation of the object """ # Create format for lattice vectors + def bcstr(bc): + left = BoundaryCondition.getitem(bc[0]).name.capitalize() + if bc[0] == bc[1]: + # single string + return left + right = BoundaryCondition.getitem(bc[1]).name.capitalize() + return f"[{left}, {right}]" s = ',\n '.join(['ABC'[i] + '=[{:.4f}, {:.4f}, {:.4f}]'.format(*self.cell[i]) for i in (0, 1, 2)]) origin = "{:.4f}, {:.4f}, {:.4f}".format(*self.origin) - return f"{self.__class__.__name__}{{nsc: {self.nsc},\n origin={origin},\n {s}\n}}" + bc = ",\n ".join(map(bcstr, self.boundary_condition)) + return f"{self.__class__.__name__}{{nsc: {self.nsc},\n origin={origin},\n {s},\n bc=[{bc}]\n}}" def __repr__(self): a, b, c, alpha, beta, gamma = map(lambda r: round(r, 4), self.parameters()) - return f"<{self.__module__}.{self.__class__.__name__} a={a}, b={b}, c={c}, α={alpha}, β={beta}, γ={gamma}, nsc={self.nsc}>" + BC = BoundaryCondition + bc = self.boundary_condition + def bcstr(bc): + left = BC.getitem(bc[0]).name[0] + if bc[0] == bc[1]: + # single string + return left + right = BC.getitem(bc[1]).name[0] + return f"[{left}, {right}]" + bc = ", ".join(map(bcstr, self.boundary_condition)) + return f"<{self.__module__}.{self.__class__.__name__} a={a}, b={b}, c={c}, α={alpha}, β={beta}, γ={gamma}, bc=[{bc}], nsc={self.nsc}>" def __eq__(self, other): """ Equality check """ @@ -1399,3 +1517,18 @@ def sc_index(self, *args, **kwargs): """ Call local `Lattice` object `sc_index` function """ return self.lattice.sc_index(*args, **kwargs) + + @property + def boundary_condition(self) -> np.ndarray: + f"{Lattice.boundary_condition.__doc__}" + return self.lattice.boundary_condition + + @boundary_condition.setter + def boundary_condition(self, boundary_condition: Sequence[BoundaryConditionType]): + raise SislError(f"Cannot use property to set boundary conditions of LatticeChild") + + @property + def pbc(self) -> np.ndarray: + f"{Lattice.pbc.__doc__}" + return self.lattice.pbc + diff --git a/src/sisl/physics/densitymatrix.py b/src/sisl/physics/densitymatrix.py index f848bbb228..99013f255d 100644 --- a/src/sisl/physics/densitymatrix.py +++ b/src/sisl/physics/densitymatrix.py @@ -11,6 +11,7 @@ from scipy.sparse import tril, triu import sisl._array as _a +from sisl import BoundaryCondition as BC from sisl import Geometry, Lattice from sisl._indices import indices_fabs_le, indices_le from sisl._internal import set_module @@ -537,7 +538,8 @@ def density(self, grid, spinor=None, tol=1e-7, eta=None): # 1. Ensure the grid has a geometry associated with it 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 == BC.PERIODIC or geometry.nsc[i] > 1 + for i, bc in enumerate(grid.lattice.boundary_condition[:, 0])] if grid.geometry is None: # Create the actual geometry that encompass the grid ia, xyz, _ = geometry.within_inf(lattice, periodic=pbc) diff --git a/src/sisl/physics/electron.py b/src/sisl/physics/electron.py index 1e9e4929e3..0e6e7fa92d 100644 --- a/src/sisl/physics/electron.py +++ b/src/sisl/physics/electron.py @@ -68,6 +68,7 @@ from scipy.sparse import csr_matrix, hstack, identity, issparse import sisl._array as _a +from sisl import BoundaryCondition as BC from sisl import Geometry, Grid, Lattice, constant, units from sisl._help import dtype_complex_to_real, dtype_real_to_complex from sisl._indices import indices_le @@ -1253,8 +1254,8 @@ def idx2spherical(ix, iy, iz, offset, dc, R): 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 == BC.PERIODIC or geometry.nsc[i] > 1 + for i, bc in enumerate(grid.lattice.boundary_condition[:, 0])] if grid.geometry is None: # Create the actual geometry that encompass the grid ia, xyz, _ = geometry.within_inf(lattice, periodic=pbc) diff --git a/src/sisl/tests/test_grid.py b/src/sisl/tests/test_grid.py index a70cb970e0..1d569ee54b 100644 --- a/src/sisl/tests/test_grid.py +++ b/src/sisl/tests/test_grid.py @@ -386,32 +386,12 @@ def test_set_grid2(self, setup): with pytest.raises(ValueError): g.set_grid([2, 2, 2, 4]) - def test_bc1(self, setup): - assert np.all(setup.g.bc == setup.g.PERIODIC) - setup.g.set_bc(a=setup.g.NEUMANN) - setup.g.set_bc(b=setup.g.NEUMANN, c=setup.g.NEUMANN) - assert np.all(setup.g.bc == setup.g.NEUMANN) - setup.g.set_bc(setup.g.PERIODIC) - assert np.all(setup.g.bc == setup.g.PERIODIC) - - def test_bc2(self, setup): - g = setup.g.copy() - P = g.PERIODIC - D = g.DIRICHLET - N = g.NEUMANN - assert np.all(g.bc == P) - g.set_bc(N) - assert np.all(g.bc == setup.g.NEUMANN) - bc = [[P, P], [N, D], [D, N]] - g.set_bc(bc) - assert np.all(g.bc == bc) - def test_argumentparser(self, setup): setup.g.ArgumentParser() def test_pyamg1(self, setup): g = setup.g.copy() - g.set_bc(g.PERIODIC) # periodic boundary conditions + g.lattice.set_boundary_condition(g.PERIODIC) # periodic boundary conditions n = np.prod(g.shape) A = csr_matrix((n, n)) b = np.zeros(A.shape[0]) @@ -439,7 +419,7 @@ def test_pyamg2(self, setup): bc = [[g.PERIODIC] * 2, [g.NEUMANN, g.DIRICHLET], [g.DIRICHLET, g.NEUMANN]] - g.set_bc(bc) + g.lattice.set_boundary_condition(bc) n = np.prod(g.shape) A = csr_matrix((n, n)) b = np.zeros(A.shape[0]) diff --git a/src/sisl/tests/test_lattice.py b/src/sisl/tests/test_lattice.py index f99d69dec4..3584e04448 100644 --- a/src/sisl/tests/test_lattice.py +++ b/src/sisl/tests/test_lattice.py @@ -9,7 +9,7 @@ import sisl import sisl.linalg as lin -from sisl import Lattice, LatticeChild, SuperCell +from sisl import BoundaryCondition, Lattice, LatticeChild, SuperCell from sisl.geom import graphene pytestmark = [pytest.mark.supercell, pytest.mark.sc, pytest.mark.lattice] @@ -487,3 +487,24 @@ def test_lattice_indices(): def test_supercell_warn(): with pytest.warns(sisl.SislDeprecation): lattice = SuperCell([1] * 3, nsc=[3, 5, 7]) + + +##### +# boundary-condition tests +##### + +@pytest.mark.parametrize("bc", list(BoundaryCondition)) +def test_lattice_bc_init(bc): + Lattice(1, boundary_condition=bc) + Lattice(1, boundary_condition=[bc, bc, bc]) + Lattice(1, boundary_condition=[[bc, bc], + [bc, bc], + [bc, bc]]) + +def test_lattice_bc_set(): + lat = Lattice(1, boundary_condition=Lattice.BC.PERIODIC) + assert lat.pbc.all() + lat.boundary_condition = Lattice.BC.UNKNOWN + assert not lat.pbc.any() + assert (lat.boundary_condition == Lattice.BC.UNKNOWN).all() + diff --git a/src/sisl_toolbox/transiesta/poisson/fftpoisson_fix.py b/src/sisl_toolbox/transiesta/poisson/fftpoisson_fix.py index d280f73b04..bdd87f3788 100644 --- a/src/sisl_toolbox/transiesta/poisson/fftpoisson_fix.py +++ b/src/sisl_toolbox/transiesta/poisson/fftpoisson_fix.py @@ -213,11 +213,13 @@ def dtype(self): # Change boundaries to always use dirichlet # This ensures that once we set the boundaries we don't # get any side-effects - periodic = grid.bc[:, 0] == grid.PERIODIC - bc = np.repeat(np.array([grid.DIRICHLET], np.int32), 6).reshape(3, 2) + BC = si.BoundaryCondition + periodic = [bc == BC.PERIODIC or geometry.nsc[i] > 1 + for i, bc in enumerate(grid.lattice.boundary_condition[:, 0])] + bc = np.repeat(np.array([BC.DIRICHLET], np.int32), 6).reshape(3, 2) for i in (0, 1, 2): if periodic[i]: - bc[i, :] = grid.PERIODIC + bc[i, :] = BC.PERIODIC grid.set_bc(bc) A, b = grid.topyamg()