From dd9e0448b5e116d28c6cdbd81b76b6cb5eae710d Mon Sep 17 00:00:00 2001 From: Nick Papior Date: Mon, 6 May 2024 15:05:08 +0200 Subject: [PATCH] fixed pbc usage across sisl, fixes #764 Now pbc usage has been moved across the code base (I hope everything is managed). I have added a setter for pbc to enable fast setting pbc. It will silently ignore any dimensions not specified. This may be used as a shorthand for lattice.set_boundary_condition(...) Signed-off-by: Nick Papior --- CHANGELOG.md | 1 + src/sisl/_core/geometry.py | 22 +++++++++++----------- src/sisl/_core/lattice.py | 10 ++++++++++ src/sisl/_core/sparse_geometry.py | 2 +- src/sisl/io/tests/test_xsf.py | 22 ++++++++++++++++++++++ src/sisl/io/xsf.py | 7 ++++--- src/sisl/physics/brillouinzone.py | 17 ++++++++++++----- 7 files changed, 61 insertions(+), 20 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 630174629..e3a66b3b5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -47,6 +47,7 @@ we hit release version 1.0.0. - A new `AtomicMatrixPlot` to plot sparse matrices, #668 ### Fixed +- 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 - all `eps` arguments has changed to `atol` diff --git a/src/sisl/_core/geometry.py b/src/sisl/_core/geometry.py index 465c90fdf..c48a49719 100644 --- a/src/sisl/_core/geometry.py +++ b/src/sisl/_core/geometry.py @@ -1724,7 +1724,7 @@ def translate2uc( Notes ----- - By default only the periodic axes will be translated to the UC. If + By default only the periodic axes (``self.pbc``) will be translated to the UC. If translation is required for all axes, supply them directly. Parameters @@ -1738,7 +1738,7 @@ def translate2uc( directions. """ if axes is None: - axes = (self.lattice.nsc > 1).nonzero()[0] + axes = self.pbc.nonzero()[0] elif isinstance(axes, bool): if axes: axes = (0, 1, 2) @@ -3265,7 +3265,7 @@ def distance( "The internal `maxR()` is negative and thus not set. " "Set an explicit value for `R`." ) - elif np.any(self.nsc > 1): + elif np.any(self.pbc): maxR = fnorm(self.cell).max() # These loops could be leveraged if we look at angles... for i, j, k in product( @@ -3366,7 +3366,7 @@ def within_inf( Note this function is rather different from `close` and `within`. Specifically this routine is returning *all* indices for the infinite - periodic system (where ``self.nsc > 1`` or `periodic` is true). + periodic system (where ``self.pbc`` or `periodic` is true). Atomic coordinates lying on the boundary of the supercell will be duplicated on the neighboring supercell images. Thus performing `geom.within_inf(geom.lattice)` @@ -3383,7 +3383,7 @@ def within_inf( the supercell in which this geometry should be expanded into. periodic : explicitly define the periodic directions, by default the periodic - directions are only where ``self.nsc > 1 & self.pbc``. + directions are only where ``self.pbc``. tol : length tolerance for the fractional coordinates to be on a duplicate site (in Ang). This allows atoms within `tol` of the cell boundaries to be taken as *inside* the @@ -3402,7 +3402,7 @@ def within_inf( """ lattice = Lattice.new(lattice) if periodic is None: - periodic = np.logical_and(self.pbc, self.nsc > 1) + periodic = self.pbc else: periodic = list(periodic) @@ -4143,7 +4143,7 @@ def dispatch(self, **kwargs): symbols=geom.atoms.Z, positions=geom.xyz.tolist(), cell=geom.cell.tolist(), - pbc=geom.nsc > 1, + pbc=geom.pbc, **kwargs, ) @@ -4174,10 +4174,10 @@ def dispatch(self, **kwargs): xyz = geom.xyz species = [PT.Z_label(Z) for Z in geom.atoms.Z] - if all(self.nsc == 1): - # we define a molecule - return Molecule(species, xyz, **kwargs) - return Structure(lattice, species, xyz, coords_are_cartesian=True, **kwargs) + if np.any(self.pbc): + return Structure(lattice, species, xyz, coords_are_cartesian=True, **kwargs) + # we define a molecule + return Molecule(species, xyz, **kwargs) to_dispatch.register("pymatgen", GeometryTopymatgenDispatch) diff --git a/src/sisl/_core/lattice.py b/src/sisl/_core/lattice.py index 82cb5de93..a772e6fbc 100644 --- a/src/sisl/_core/lattice.py +++ b/src/sisl/_core/lattice.py @@ -179,6 +179,16 @@ def pbc(self) -> np.ndarray: # along the same lattice vector. So checking one should suffice return self._bc[:, 0] == BoundaryCondition.PERIODIC + @pbc.setter + def pbc(self, pbc) -> None: + """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 + assert len(pbc) == 3 + for axis, bc in enumerate(pbc): + if bc: + self._bc[axis] = BoundaryCondition.PERIODIC + @property def origin(self) -> ndarray: """Origin for the cell""" diff --git a/src/sisl/_core/sparse_geometry.py b/src/sisl/_core/sparse_geometry.py index faff2ff16..6d2a76cd9 100644 --- a/src/sisl/_core/sparse_geometry.py +++ b/src/sisl/_core/sparse_geometry.py @@ -182,7 +182,7 @@ def translate2uc(self, atoms: AtomsIndex = None, axes: Optional[CellAxes] = None """ # Sanitize the axes argument if axes is None: - axes = (self.lattice.nsc > 1).nonzero()[0] + axes = self.lattice.pbc.nonzero()[0] elif isinstance(axes, bool): if axes: diff --git a/src/sisl/io/tests/test_xsf.py b/src/sisl/io/tests/test_xsf.py index 43d4081cb..5a237d580 100644 --- a/src/sisl/io/tests/test_xsf.py +++ b/src/sisl/io/tests/test_xsf.py @@ -24,6 +24,28 @@ def test_default(sisl_tmp): assert grid.geometry is None +@pytest.mark.parametrize( + "pbc", + [ + (True, True, True), + (True, True, False), + (True, False, False), + (False, False, False), + ], +) +def test_pbc(sisl_tmp, pbc): + f = sisl_tmp("GRID_default.xsf", _dir) + geom = Geometry( + np.random.rand(10, 3), + np.random.randint(1, 70, 10), + lattice=[10, 10, 10, 45, 60, 90], + ) + geom.lattice.pbc = pbc + geom.write(f) + geom2 = geom.read(f) + assert all(geom.pbc == geom2.pbc) + + def test_default_size(sisl_tmp): f = sisl_tmp("GRID_default_size.xsf", _dir) grid = Grid(0.2, lattice=2.0) diff --git a/src/sisl/io/xsf.py b/src/sisl/io/xsf.py index a007696a6..eea6541e3 100644 --- a/src/sisl/io/xsf.py +++ b/src/sisl/io/xsf.py @@ -122,11 +122,12 @@ def write_lattice(self, lattice: Lattice, fmt: str = ".8f"): "# File created by: sisl {}\n#\n".format(strftime("%Y-%m-%d", gmtime())) ) - if all(lattice.nsc == 1): + pbc = lattice.pbc + if pbc.sum() == 0: self._write_once("MOLECULE\n#\n") - elif all(lattice.nsc[:2] > 1): + elif all(pbc == (True, True, False)): self._write_once("SLAB\n#\n") - elif lattice.nsc[0] > 1: + elif all(pbc == (True, False, False)): self._write_once("POLYMER\n#\n") else: self._write_once("CRYSTAL\n#\n") diff --git a/src/sisl/physics/brillouinzone.py b/src/sisl/physics/brillouinzone.py index d4fc20d66..4047a59df 100644 --- a/src/sisl/physics/brillouinzone.py +++ b/src/sisl/physics/brillouinzone.py @@ -389,7 +389,14 @@ def merge(bzs, weight_scale: Union[Sequence[float], float] = 1.0, parent=None): return BrillouinZone(parent, np.concatenate(k), np.concatenate(w)) - def volume(self, ret_dim: bool = False, periodic=None): + @deprecate_argument( + "periodic", + "axes", + "argument 'periodic' has been deprecated in favor of 'axes', please update your code.", + "0.15", + "0.16", + ) + def volume(self, ret_dim: bool = False, axes: Optional[AnyAxes] = None): """Calculate the volume of the full Brillouin zone of the parent This will return the volume depending on the dimensions of the system. @@ -400,11 +407,11 @@ def volume(self, ret_dim: bool = False, periodic=None): Parameters ---------- - ret_dim: + ret_dim : also return the dimensionality of the system - periodic : array_like of int, optional + axes : estimate the volume using only the directions indexed by this array. - The default value is `(self.parent.nsc > 1).nonzero()[0]`. + The default value is ``self.parent.pbc.nonzero()[0]``. Returns ------- @@ -416,7 +423,7 @@ def volume(self, ret_dim: bool = False, periodic=None): """ # default periodic array if periodic is None: - periodic = (self.parent.nsc > 1).nonzero()[0] + periodic = self.parent.pbc.nonzero()[0] dim = len(periodic) vol = 0.0