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 <nickpapior@gmail.com>

CHANGELOG.md

@@ -48,6 +48,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`

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)
@@ -3279,7 +3279,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(
@@ -3387,7 +3387,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)`
@@ -3404,7 +3404,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``.
         atol :
             length tolerance for the coordinates to be on a duplicate site (in Ang).
             This allows atoms within `atol` of the cell boundaries to be taken as *inside* the
@@ -3423,7 +3423,7 @@ def within_inf(
         """
         lattice = self.lattice.__class__.new(lattice)
         if periodic is None:
-            periodic = np.logical_and(self.pbc, self.nsc > 1)
+            periodic = self.pbc
         else:
             periodic = list(periodic)
 
@@ -4172,7 +4172,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,
         )
 
@@ -4203,10 +4203,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)

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

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:

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)

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

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