diff --git a/CHANGELOG.md b/CHANGELOG.md
index f4597b0970..7be9e4abf2 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -52,6 +52,8 @@ we hit release version 1.0.0.
 - removed `Selector` and `TimeSelector`, they were never used internally
 
 ### Changed
+- `stdoutSileSiesta.read_*` now defaults to read the *next* entry, and not the last
+- `stdoutSileSiesta.read_*` changed MD output functionality, see #586 for details
 - `AtomNeighbours` changed name to `AtomNeighbor` to follow #393
 - removed `Lattice.translate|move`, they did not make sense, and so their
   usage should be deferred to `Lattice.add` instead.
@@ -66,7 +68,7 @@ we hit release version 1.0.0.
 
 ### Added
 - Creation of honeycomb flakes (`sisl.geom.honeycomb_flake`,
-  `sisl.geom.graphene_flake`). #636
+  `sisl.geom.graphene_flake`), #636
 - added `Geometry.as_supercell` to create the supercell structure,
   thanks to @pfebrer for the suggestion
 - added `Lattice.to` and `Lattice.new` to function the same
diff --git a/src/sisl/_core/lattice.py b/src/sisl/_core/lattice.py
index 99ec085072..4f37316ef3 100644
--- a/src/sisl/_core/lattice.py
+++ b/src/sisl/_core/lattice.py
@@ -23,7 +23,7 @@
 from sisl._dispatcher import AbstractDispatch, ClassDispatcher, TypeDispatcher
 from sisl._internal import set_module
 from sisl._math_small import cross3, dot3
-from sisl.messages import SislError, deprecate, deprecation, info, warn
+from sisl.messages import SislError, deprecate, deprecation, warn
 from sisl.shape.prism4 import Cuboid
 from sisl.typing import Axes, Axies, Axis
 from sisl.utils.mathematics import fnorm
@@ -274,7 +274,7 @@ def conv(v):
                     "must have that BC."
                 )
             if changed.any() and (~bc).all() and nsc > 1:
-                info(
+                warn(
                     f"{self.__class__.__name__}.set_boundary_condition is having image connections (nsc={nsc}>1) "
                     "while having a non-periodic boundary condition."
                 )
diff --git a/src/sisl/_core/tests/test_lattice.py b/src/sisl/_core/tests/test_lattice.py
index fb0b9318c4..f91f757cfd 100644
--- a/src/sisl/_core/tests/test_lattice.py
+++ b/src/sisl/_core/tests/test_lattice.py
@@ -563,7 +563,7 @@ def test_lattice_bc_fail():
 
 def test_lattice_info():
     lat = Lattice(1, nsc=[3, 3, 3])
-    with pytest.warns(sisl.SislInfo) as record:
+    with pytest.warns(sisl.SislWarning) as record:
         lat.set_boundary_condition(b=Lattice.BC.DIRICHLET)
         lat.set_boundary_condition(c=Lattice.BC.PERIODIC)
         assert len(record) == 1
diff --git a/src/sisl/io/_multiple.py b/src/sisl/io/_multiple.py
index 89b69d0021..11aed9c4b8 100644
--- a/src/sisl/io/_multiple.py
+++ b/src/sisl/io/_multiple.py
@@ -2,6 +2,7 @@
 # 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/.
 from functools import reduce, update_wrapper
+from itertools import zip_longest
 from numbers import Integral
 from textwrap import dedent
 from typing import Any, Callable, Optional, Type
@@ -9,6 +10,35 @@
 Func = Callable[..., Optional[Any]]
 
 
+def postprocess_tuple(*funcs):
+    """Post-processes the returned value according to the funcs for multiple data
+
+    The internal algorithm will use zip_longest to apply the *last* `funcs[-1]` to
+    the remaining return functions. Useful if there are many tuples that can
+    be handled equivalently:
+
+    Examples
+    --------
+    >>> postprocess_tuple(np.array) == postprocess_tuple(np.array, np.array)
+    """
+
+    def post(ret):
+        nonlocal funcs
+        if isinstance(ret[0], tuple):
+            return tuple(
+                func(r)
+                for r, func in zip_longest(zip(*ret), funcs, fillvalue=funcs[-1])
+            )
+        return funcs[0](ret)
+
+    return post
+
+
+def is_sliceable(method: Func):
+    """Check whether a function is implemented with the `SileBinder` decorator"""
+    return isinstance(method, (SileBinder, SileBound))
+
+
 class SileSlicer:
     """Handling io-methods in sliced behaviour for multiple returns
 
@@ -26,6 +56,7 @@ def __init__(
         func: Func,
         key: Type[Any],
         *,
+        check_empty: Optional[Func] = None,
         skip_func: Optional[Func] = None,
         postprocess: Optional[Callable[..., Any]] = None,
     ):
@@ -39,6 +70,16 @@ def __init__(
             self.skip_func = func
         else:
             self.skip_func = skip_func
+
+        if check_empty is None:
+
+            def check_empty(r):
+                if isinstance(r, tuple):
+                    return reduce(lambda x, y: x and y is None, r, True)
+                return r is None
+
+        self.check_empty = check_empty
+
         if postprocess is None:
 
             def postprocess(ret):
@@ -60,11 +101,6 @@ def __call__(self, *args, **kwargs):
 
         inf = 100000000000000
 
-        def check_none(r):
-            if isinstance(r, tuple):
-                return reduce(lambda x, y: x and y is None, r, True)
-            return r is None
-
         # Determine whether we can reduce the call overheads
         start = 0
         stop = inf
@@ -100,7 +136,7 @@ def check_none(r):
 
             # now do actual parsing
             retval = func(obj, *args, **kwargs)
-            while not check_none(retval):
+            while not self.check_empty(retval):
                 append(retval)
                 if len(retvals) >= stop:
                     # quick exit
@@ -112,11 +148,10 @@ def check_none(r):
                 return None
 
         # ensure the next call won't use this key
-        # This will prohibit the use
+        # This will enable the use
         # tmp = sile.read_geometry[:10]
-        # tmp() # will return the first 10
-        # tmp() # will return the default (single) item
-        self.key = None
+        # tmp() # will returns the first 10
+        # tmp() # will returns the next 10
         if isinstance(key, Integral):
             return retvals[key]
 
diff --git a/src/sisl/io/cube.py b/src/sisl/io/cube.py
index cf20473591..ef8750f1c0 100644
--- a/src/sisl/io/cube.py
+++ b/src/sisl/io/cube.py
@@ -32,11 +32,11 @@ class cubeSile(Sile):
     )
     def write_lattice(
         self,
-        lattice,
-        fmt="15.10e",
+        lattice: Lattice,
+        fmt: str = "15.10e",
         size=None,
         origin=None,
-        unit="Bohr",
+        unit: str = "Bohr",
         *args,
         **kwargs,
     ):
@@ -44,15 +44,15 @@ def write_lattice(
 
         Parameters
         ----------
-        lattice : Lattice
+        lattice :
             lattice to be written
-        fmt : str, optional
+        fmt :
             floating point format for stored values
         size : (3, ), optional
             shape of the stored grid (``[1, 1, 1]``)
         origin : (3, ), optional
             origin of the cell (``[0, 0, 0]``)
-        unit: str, optional
+        unit:
             what length unit should the cube file data be written in
         """
         sile_raise_write(self)
@@ -83,11 +83,11 @@ def write_lattice(
     @sile_fh_open()
     def write_geometry(
         self,
-        geometry,
-        fmt="15.10e",
+        geometry: Geometry,
+        fmt: str = "15.10e",
         size=None,
         origin=None,
-        unit="Bohr",
+        unit: str = "Bohr",
         *args,
         **kwargs,
     ):
@@ -95,15 +95,15 @@ def write_geometry(
 
         Parameters
         ----------
-        geometry : Geometry
+        geometry :
             geometry to be written
-        fmt : str, optional
+        fmt :
             floating point format for stored values
         size : (3, ), optional
             shape of the stored grid (``[1, 1, 1]``)
         origin : (3, ), optional
             origin of the cell (``[0, 0, 0]``)
-        unit: str, optional
+        unit:
             what length unit should the cube file data be written in
         """
         sile_raise_write(self)
@@ -140,24 +140,32 @@ def write_geometry(
             )
 
     @sile_fh_open()
-    def write_grid(self, grid, fmt=".5e", imag=False, unit="Bohr", *args, **kwargs):
+    def write_grid(
+        self,
+        grid: Grid,
+        fmt: str = ".5e",
+        imag: bool = False,
+        unit: str = "Bohr",
+        *args,
+        **kwargs,
+    ):
         """Write `Grid` to the contained file
 
         Parameters
         ----------
-        grid : Grid
+        grid :
            the grid to be written in the CUBE file
-        fmt : str, optional
+        fmt :
            format used for precision output
-        imag : bool, optional
+        imag :
            write only imaginary part of the grid, default to only writing the
            real part.
-        buffersize : int, optional
-           size of the buffer while writing the data, (6144)
-        unit: str, optional
+        unit:
             what length unit should the cube file data be written in.
             The grid data is assumed to be unit-less, this unit only refers
             to the lattice vectors and atomic coordinates.
+        buffersize : int, optional
+           size of the buffer while writing the data, (6144)
         """
         # Check that we can write to the file
         sile_raise_write(self)
@@ -220,7 +228,7 @@ def _r_header_dict(self):
         return header
 
     @sile_fh_open()
-    def read_lattice(self, na=False):
+    def read_lattice(self, na: bool = False):
         """Returns `Lattice` object from the CUBE file
 
         Parameters
diff --git a/src/sisl/io/siesta/bands.py b/src/sisl/io/siesta/bands.py
index 820e2bf463..bf0adbd0d6 100644
--- a/src/sisl/io/siesta/bands.py
+++ b/src/sisl/io/siesta/bands.py
@@ -23,7 +23,7 @@ def read_fermi_level(self):
         return float(self.readline())
 
     @sile_fh_open()
-    def read_data(self, as_dataarray=False):
+    def read_data(self, as_dataarray: bool = False):
         """Returns data associated with the bands file
 
         The energy levels are shifted with respect to the Fermi-level.
diff --git a/src/sisl/io/siesta/fa.py b/src/sisl/io/siesta/fa.py
index eb3fb879e1..20ee40e17f 100644
--- a/src/sisl/io/siesta/fa.py
+++ b/src/sisl/io/siesta/fa.py
@@ -28,7 +28,7 @@ def read_force(self):
         return f
 
     @sile_fh_open()
-    def write_force(self, f, fmt=".9e"):
+    def write_force(self, f, fmt: str = ".9e"):
         """Write forces to file
 
         Parameters
diff --git a/src/sisl/io/siesta/fc.py b/src/sisl/io/siesta/fc.py
index 0ab07cdfa6..c1f6e66ada 100644
--- a/src/sisl/io/siesta/fc.py
+++ b/src/sisl/io/siesta/fc.py
@@ -1,6 +1,8 @@
 # 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/.
+from typing import Optional
+
 import numpy as np
 
 from sisl._internal import set_module
@@ -18,7 +20,9 @@ class fcSileSiesta(SileSiesta):
     """Force constant file"""
 
     @sile_fh_open()
-    def read_force(self, displacement=None, na=None):
+    def read_force(
+        self, displacement: Optional[float] = None, na: Optional[int] = None
+    ):
         """Reads all displacement forces by multiplying with the displacement value
 
         Since the force constant file does not contain the non-displaced configuration
@@ -33,12 +37,12 @@ def read_force(self, displacement=None, na=None):
 
         Parameters
         ----------
-        displacement : float, optional
+        displacement :
            the used displacement in the calculation, since Siesta 4.1-b4 this value
            is written in the FC file and hence not required.
            If prior Siesta versions are used and this is not supplied the 0.04 Bohr displacement
            will be assumed.
-        na : int, optional
+        na :
            number of atoms in geometry (for returning correct number of atoms), since Siesta 4.1-b4
            this value is written in the FC file and hence not required.
            If prior Siesta versions are used then the file is expected to only contain 1-atom displacement.
@@ -67,12 +71,12 @@ def read_force(self, displacement=None, na=None):
         return self.read_force_constant(na) * displacement.reshape(1, 3, 2, 1, 1)
 
     @sile_fh_open()
-    def read_force_constant(self, na=None):
+    def read_force_constant(self, na: Optional[int] = None):
         """Reads the force-constant stored in the FC file
 
         Parameters
         ----------
-        na : int, optional
+        na :
            number of atoms in the unit-cell, if not specified it will guess on only
            one atom displacement.
 
diff --git a/src/sisl/io/siesta/fdf.py b/src/sisl/io/siesta/fdf.py
index f799ef7293..289657b0dc 100644
--- a/src/sisl/io/siesta/fdf.py
+++ b/src/sisl/io/siesta/fdf.py
@@ -7,6 +7,7 @@
 import warnings
 from datetime import datetime
 from os.path import isfile
+from typing import Any, Optional
 
 import numpy as np
 import scipy as sp
@@ -18,6 +19,7 @@
     Atoms,
     DynamicalMatrix,
     Geometry,
+    Grid,
     Lattice,
     Orbital,
     SphericalOrbital,
@@ -25,8 +27,8 @@
 )
 from sisl._indices import indices_only
 from sisl._internal import set_module
-from sisl.messages import SislError, info, warn
-from sisl.physics import BandStructure
+from sisl.messages import SislError, deprecate_argument, info, warn
+from sisl.physics import BandStructure, BrillouinZone
 from sisl.unit.siesta import unit_convert, unit_default, unit_group, units
 from sisl.utils.cmd import default_ArgumentParser, default_namespace
 from sisl.utils.mathematics import fnorm
@@ -232,7 +234,7 @@ def add(f):
         return includes
 
     @sile_fh_open()
-    def _read_label(self, label):
+    def _read_label(self, label: str):
         """Try and read the first occurence of a key
 
         This will take care of blocks, labels and piped in labels
@@ -329,7 +331,7 @@ def process_line(line):
         return l
 
     @classmethod
-    def _type(cls, value):
+    def _type(cls, value: Any):
         """Determine the type by the value
 
         Parameters
@@ -377,7 +379,7 @@ def _type(cls, value):
         return "n"
 
     @sile_fh_open()
-    def type(self, label):
+    def type(self, label: str):
         """Return the type of the fdf-keyword
 
         Parameters
@@ -389,7 +391,13 @@ def type(self, label):
         return self._type(self._read_label(label))
 
     @sile_fh_open()
-    def get(self, label, default=None, unit=None, with_unit=False):
+    def get(
+        self,
+        label: str,
+        default: Optional[Any] = None,
+        unit: Optional[str] = None,
+        with_unit: bool = False,
+    ):
         """Retrieve fdf-keyword from the file
 
         Parameters
@@ -475,7 +483,7 @@ def get(self, label, default=None, unit=None, with_unit=False):
 
         return value
 
-    def set(self, key, value, keep=True):
+    def set(self, key: str, value: Any, keep: bool = True):
         """Add the key and value to the FDF file
 
         Parameters
@@ -543,7 +551,7 @@ def write(fh, value):
         self._open()
 
     @staticmethod
-    def print(key, value):
+    def print(key: str, value: Any):
         """Return a string which is pretty-printing the key+value"""
         if isinstance(value, list):
             s = f"%block {key}"
@@ -568,14 +576,17 @@ def print(key, value):
         return s
 
     @sile_fh_open()
-    def write_lattice(self, sc, fmt=".8f", *args, **kwargs):
+    @deprecate_argument(
+        "sc", "lattice", "use lattice= instead of sc=", from_version="0.15"
+    )
+    def write_lattice(self, lattice: Lattice, fmt: str = ".8f", *args, **kwargs):
         """Writes the supercell
 
         Parameters
         ----------
-        sc : Lattice
+        lattice:
            supercell object to write
-        fmt : str, optional
+        fmt :
            precision used to store the lattice vectors
         unit : {"Ang", "Bohr"}
            the unit used when writing the data.
@@ -594,20 +605,20 @@ def write_lattice(self, sc, fmt=".8f", *args, **kwargs):
         # Write out the cell
         self._write(f"LatticeConstant 1.0 {unit}\n")
         self._write("%block LatticeVectors\n")
-        self._write(fmt_str.format(*sc.cell[0, :] * conv))
-        self._write(fmt_str.format(*sc.cell[1, :] * conv))
-        self._write(fmt_str.format(*sc.cell[2, :] * conv))
+        self._write(fmt_str.format(*lattice.cell[0, :] * conv))
+        self._write(fmt_str.format(*lattice.cell[1, :] * conv))
+        self._write(fmt_str.format(*lattice.cell[2, :] * conv))
         self._write("%endblock LatticeVectors\n")
 
     @sile_fh_open()
-    def write_geometry(self, geometry, fmt=".8f", *args, **kwargs):
+    def write_geometry(self, geometry: Geometry, fmt: str = ".8f", *args, **kwargs):
         """Writes the geometry
 
         Parameters
         ----------
-        geometry : Geometry
+        geometry :
            geometry object to write
-        fmt : str, optional
+        fmt :
            precision used to store the atomic coordinates
         unit : {"Ang", "Bohr", "fractional", "frac"}
            the unit used when writing the data.
@@ -683,7 +694,7 @@ def write_block(atoms, append, write_block):
             self._write("%endblock\n")
 
     @sile_fh_open()
-    def write_brillouinzone(self, bz):
+    def write_brillouinzone(self, bz: BrillouinZone):
         r"""Writes Brillouinzone information to the fdf input file
 
         The `bz` object will be written as options in the input file.
@@ -692,7 +703,7 @@ def write_brillouinzone(self, bz):
 
         Parameters
         ----------
-        bz : BrillouinZone
+        bz :
             which setting to write to the file
 
         Notes
@@ -754,7 +765,7 @@ def _r_lattice_nsc_orb_indx(self, *args, **kwargs):
             return orbindxSileSiesta(f).read_lattice_nsc()
         return None
 
-    def read_lattice(self, output=False, *args, **kwargs):
+    def read_lattice(self, output: bool = False, *args, **kwargs) -> Lattice:
         """Returns Lattice object by reading fdf or Siesta output related files.
 
         One can limit the tried files to only one file by passing
@@ -762,7 +773,7 @@ def read_lattice(self, output=False, *args, **kwargs):
 
         Parameters
         ----------
-        output: bool, optional
+        output:
             whether to read supercell from output files (default to read from
             the fdf file).
         order: list of str, optional
@@ -1401,7 +1412,7 @@ def _dynamical_matrix_from_fc(self, geom, FC, FC_atoms, *args, **kwargs):
 
         return D
 
-    def read_geometry(self, output=False, *args, **kwargs):
+    def read_geometry(self, output: bool = False, *args, **kwargs) -> Geometry:
         """Returns Geometry object by reading fdf or Siesta output related files.
 
         One can limit the tried files to only one file by passing
@@ -1409,7 +1420,7 @@ def read_geometry(self, output=False, *args, **kwargs):
 
         Parameters
         ----------
-        output: bool, optional
+        output:
             whether to read geometry from output files (default to read from
             the fdf file).
         order: list of str, optional
@@ -1647,7 +1658,7 @@ def _r_geometry_fdf(self, *args, **kwargs):
 
         return geom
 
-    def read_grid(self, name, *args, **kwargs):
+    def read_grid(self, name: str, *args, **kwargs) -> Grid:
         """Read grid related information from any of the output files
 
         The order of the readed data is shown below.
diff --git a/src/sisl/io/siesta/kp.py b/src/sisl/io/siesta/kp.py
index bcce731be3..edb0e23211 100644
--- a/src/sisl/io/siesta/kp.py
+++ b/src/sisl/io/siesta/kp.py
@@ -1,10 +1,14 @@
 # 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/.
+from typing import Optional, Union
+
 import numpy as np
 
+from sisl import Lattice, LatticeChild
 from sisl._internal import set_module
 from sisl.messages import deprecate_argument
+from sisl.physics.brillouinzone import BrillouinZone
 from sisl.unit.siesta import unit_convert
 
 from ..sile import add_sile, sile_fh_open, sile_raise_write
@@ -14,6 +18,8 @@
 
 Bohr2Ang = unit_convert("Bohr", "Ang")
 
+TLattice = Optional[Union[Lattice, LatticeChild]]
+
 
 @set_module("sisl.io.siesta")
 class kpSileSiesta(SileSiesta):
@@ -23,12 +29,12 @@ class kpSileSiesta(SileSiesta):
     @deprecate_argument(
         "sc", "lattice", "use lattice= instead of sc=", from_version="0.15"
     )
-    def read_data(self, lattice=None):
+    def read_data(self, lattice: TLattice = None):
         """Returns K-points from the file (note that these are in reciprocal units)
 
         Parameters
         ----------
-        lattice : LatticeChild, optional
+        lattice :
            if supplied the returned k-points will be in reduced coordinates
 
         Returns
@@ -53,7 +59,7 @@ def read_data(self, lattice=None):
         return np.dot(k, lattice.cell.T / (2 * np.pi)), w
 
     @sile_fh_open()
-    def write_data(self, k, weight, fmt=".9e"):
+    def write_data(self, k, weight, fmt: str = ".9e"):
         """Writes K-points to file
 
         Parameters
@@ -78,7 +84,7 @@ def write_data(self, k, weight, fmt=".9e"):
     @deprecate_argument(
         "sc", "lattice", "use lattice= instead of sc=", from_version="0.15"
     )
-    def read_brillouinzone(self, lattice):
+    def read_brillouinzone(self, lattice: TLattice):
         """Returns K-points from the file (note that these are in reciprocal units)
 
         Parameters
@@ -99,7 +105,7 @@ def read_brillouinzone(self, lattice):
         return bz
 
     @sile_fh_open()
-    def write_brillouinzone(self, bz, fmt=".9e"):
+    def write_brillouinzone(self, bz: BrillouinZone, fmt: str = ".9e"):
         """Writes BrillouinZone-points to file
 
         Parameters
@@ -146,7 +152,7 @@ def read_data(self):
     @deprecate_argument(
         "sc", "lattice", "use lattice= instead of sc=", from_version="0.15"
     )
-    def read_brillouinzone(self, lattice):
+    def read_brillouinzone(self, lattice: TLattice):
         """Returns K-points from the file
 
         Parameters
@@ -167,7 +173,7 @@ def read_brillouinzone(self, lattice):
         return bz
 
     @sile_fh_open()
-    def write_brillouinzone(self, bz, fmt=".9e"):
+    def write_brillouinzone(self, bz: BrillouinZone, fmt: str = ".9e"):
         """Writes BrillouinZone-points to file
 
         Parameters
diff --git a/src/sisl/io/siesta/pdos.py b/src/sisl/io/siesta/pdos.py
index 2256b513d0..9be33428fa 100644
--- a/src/sisl/io/siesta/pdos.py
+++ b/src/sisl/io/siesta/pdos.py
@@ -53,7 +53,7 @@ def read_fermi_level(self):
         return Ef
 
     @sile_fh_open(True)
-    def read_data(self, as_dataarray=False):
+    def read_data(self, as_dataarray: bool = False):
         r"""Returns data associated with the PDOS file
 
         For spin-polarized calculations the returned values are up/down, orbitals, energy.
@@ -390,7 +390,7 @@ def Eindex(e):
             "-E",
             action=ERange,
             help="""Denote the sub-section of energies that are extracted: "-1:0,1:2" [eV]
-                       
+
                        This flag takes effect on all energy-resolved quantities and is reset whenever --plot or --out is called""",
         )
 
diff --git a/src/sisl/io/siesta/stdout.py b/src/sisl/io/siesta/stdout.py
index f15327e645..e5fa65f72f 100644
--- a/src/sisl/io/siesta/stdout.py
+++ b/src/sisl/io/siesta/stdout.py
@@ -3,12 +3,14 @@
 # file, You can obtain one at https://mozilla.org/MPL/2.0/.
 import os
 from functools import lru_cache
+from typing import Optional
 
 import numpy as np
 
 import sisl._array as _a
 from sisl import Atom, Geometry, Lattice
 from sisl._common import Opt
+from sisl._help import voigt_matrix
 from sisl._internal import set_module
 from sisl.messages import deprecation, warn
 from sisl.physics import Spin
@@ -16,6 +18,7 @@
 from sisl.utils import PropertyDict
 from sisl.utils.cmd import *
 
+from .._multiple import SileBinder, postprocess_tuple
 from ..sile import SileError, add_sile, sile_fh_open
 from .sile import SileSiesta
 
@@ -48,6 +51,56 @@ def _parse_spin(attr, match):
     return Spin()
 
 
+def _read_scf_empty(scf):
+    if isinstance(scf, tuple):
+        return len(scf[0]) == 0
+    return len(scf) == 0
+
+
+def _read_scf_md_process(scfs):
+
+    if len(scfs) == 0:
+        return None
+
+    if not isinstance(scfs, list):
+        # single MD request either as:
+        #  - np.ndarray
+        #  - np.ndarray, tuple
+        #  - pd.DataFrame
+        return scfs
+
+    has_props = isinstance(scfs[0], tuple)
+    if has_props:
+        my_len = lambda scf: len(scf[0])
+    else:
+        my_len = len
+
+    scf_len1 = np.all(_a.fromiterd(map(my_len, scfs)) == 1)
+    if isinstance(scfs[0], (np.ndarray, tuple)):
+
+        if has_props:
+            props = scfs[0][1]
+            scfs = [scf[0] for scf in scfs]
+
+        if scf_len1:
+            scfs = np.array(scfs)
+        if has_props:
+            return scfs, props
+        return scfs
+
+    # We are dealing with a dataframe
+    import pandas as pd
+
+    df = pd.concat(
+        scfs,
+        keys=_a.arangei(1, len(scfs) + 1),
+        names=["imd"],
+    )
+    if scf_len1:
+        df.reset_index("iscf", inplace=True)
+    return df
+
+
 @set_module("sisl.io.siesta")
 class stdoutSileSiesta(SileSiesta):
     """Output file from Siesta
@@ -67,6 +120,12 @@ class stdoutSileSiesta(SileSiesta):
             r"^redata: Spin configuration",
             _parse_spin,
         ),
+        _A(
+            "_final_analysis",
+            r"^siesta: Final energy",
+            lambda attr, match: lambda: True,
+            default=lambda: False,
+        ),
     ]
 
     @deprecation(
@@ -236,17 +295,17 @@ def _r_geometry_atomic(self, line, atoms=None):
 
         return Geometry(xyz, atoms, lattice=cell)
 
+    @SileBinder()
     @sile_fh_open()
-    def read_geometry(self, last=True, all=False):
+    def read_geometry(self, skip_input: bool = True):
         """Reads the geometry from the Siesta output file
 
         Parameters
         ----------
-        last: bool, optional
-           only read the last geometry
-        all: bool, optional
-           return a list of all geometries (like an MD)
-           If `True` `last` is ignored
+        skip_input :
+            the input geometry may be contained as a print-out.
+            This is not part of an MD calculation, and hence is per
+            default not returned.
 
         Returns
         -------
@@ -255,62 +314,30 @@ def read_geometry(self, last=True, all=False):
              a list of geometries corresponding to the MD-runs.
         """
         atoms = self.read_basis()
-        if all:
-            # force last to be false
-            last = False
 
-        def func_none(*args, **kwargs):
+        def func(*args, **kwargs):
             """Wrapper to return None"""
             return None
 
-        def next_geom():
-            coord, func = 0, func_none
-            line = " "
-            while coord == 0 and line != "":
-                line = self.readline()
-                if "outcoor" in line and "coordinates" in line:
-                    coord, func = 1, self._r_geometry_outcoor
-                elif "siesta: Atomic coordinates" in line:
-                    coord, func = 2, self._r_geometry_atomic
-            return coord, func(line, atoms)
-
-        # Read until a coordinate block is found
-        geom0 = None
-        mds = []
-
-        if all or last:
-            # we need to read through all things!
-            while True:
-                coord, geom = next_geom()
-                if coord == 0:
-                    break
-                if coord == 2:
-                    geom0 = geom
-                else:
-                    mds.append(geom)
-
-            # Since the user requests only the MD geometries
-            # we only return those
-            if last:
-                if len(mds) > 0:
-                    return mds[-1]
-                return geom0
-            return mds
+        line = " "
+        while line != "":
+            line = self.readline()
+            if "outcoor" in line and "coordinates" in line:
+                func = self._r_geometry_outcoor
+                break
+            elif "siesta: Atomic coordinates" in line and not skip_input:
+                func = self._r_geometry_atomic
+                break
 
-        # just read the next geometry we hit
-        return next_geom()[1]
+        return func(line, atoms)
 
-    @sile_fh_open(True)
-    def read_force(self, last=True, all=False, total=False, max=False, key="siesta"):
+    @SileBinder(postprocess=postprocess_tuple(_a.arrayd))
+    @sile_fh_open()
+    def read_force(self, total: bool = False, max: bool = False, key: str = "siesta"):
         """Reads the forces from the Siesta output file
 
         Parameters
         ----------
-        last: bool, optional
-           only read the last force
-        all: bool, optional
-           return a list of all forces (like an MD)
-           If `True` `last` is ignored
         total: bool, optional
             return the total forces instead of the atomic forces.
         max: bool, optional
@@ -337,97 +364,68 @@ def read_force(self, last=True, all=False, total=False, max=False, key="siesta")
                 - total: (nMDsteps, 3)
                 - max: (nMDsteps, )
 
-            If `all` is `False`, the first dimension does not exist. In the case of max, the returned value
-            will therefore be just a float, not an array.
-
             If `total` and `max` are both `True`, they are returned separately as a tuple: ``(total, max)``
         """
-        if all:
-            last = False
-
         # Read until forces are found
-        def next_force():
-            found, line = self.step_to(f"{key}: Atomic forces", allow_reread=False)
-            if not found:
-                return None
+        found, line = self.step_to(f"{key}: Atomic forces", allow_reread=False)
+        if not found:
+            return None
 
-            # Now read data
-            F = []
-            line = self.readline()
-            if "siesta:" in line:
-                # This is the final summary, we don't need to read it as it does not contain new information
-                # and also it make break things since max forces are not written there
-                return None
-
-            # First, we encounter the atomic forces
-            while "---" not in line:
-                line = line.split()
-                if not (total or max):
-                    F.append([float(x) for x in line[-3:]])
-                line = self.readline()
-                if line == "":
-                    break
+        # Now read data
+        line = self.readline()
+        if "siesta:" in line:
+            # This is the final summary, we don't need to read it as it does not contain new information
+            # and also it make break things since max forces are not written there
+            return None
 
+        F = []
+        # First, we encounter the atomic forces
+        while "---" not in line:
+            line = line.split()
+            if not (total or max):
+                F.append([float(x) for x in line[-3:]])
             line = self.readline()
-            # Then, the total forces
-            if total:
-                F = [float(x) for x in line.split()[-3:]]
+            if line == "":
+                break
 
-            line = self.readline()
-            # And after that we can read the max force
-            if max and len(line.split()) != 0:
-                line = self.readline()
-                maxF = float(line.split()[1])
-
-                # In case total is also requested, we are going to store it all in the same variable
-                # It will be separated later
-                if total:
-                    F = (*F, maxF)
-                else:
-                    F = maxF
-
-            return _a.arrayd(F)
-
-        def return_forces(Fs):
-            # Handle cases where we can't now if they are found
-            if Fs is None:
-                return None
-            Fs = _a.arrayd(Fs)
-            if max and total:
-                return (Fs[..., :-1], Fs[..., -1])
-            elif max and not all:
-                return Fs.ravel()[0]
-            return Fs
-
-        if all or last:
-            # list of all forces
-            Fs = []
-            while True:
-                F = next_force()
-                if F is None:
-                    break
-                Fs.append(F)
+        if not F:
+            F = None
 
-            if last:
-                return return_forces(Fs[-1])
+        line = self.readline().split()
+        # Parse total forces if requested
+        if total and line:
+            F = _a.arrayd([float(x) for x in line[-3:]])
+
+        # And after that we can read the max force
+        line = self.readline()
+        if max and line:
+            line = self.readline().split()
+            maxF = _a.arrayd(float(line[1]))
 
-            return return_forces(Fs)
+            # In case total is also requested, we are going to
+            # store it all in the same variable.
+            # It will be separated later
+            if total:
+                # create tuple
+                F = (F, maxF)
+            else:
+                F = maxF
 
-        return return_forces(next_force())
+        return F
 
-    @sile_fh_open(True)
-    def read_stress(self, key="static", last=True, all=False):
+    @SileBinder(postprocess=_a.arrayd)
+    @sile_fh_open()
+    def read_stress(self, key: str = "static", skip_final: bool = True):
         """Reads the stresses from the Siesta output file
 
         Parameters
         ----------
-        key : {'static', 'total'}
+        key : {static, total, Voigt}
            which stress to read from the output.
-        last: bool, optional
-           only read the last stress
-        all: bool, optional
-           return a list of all stresses (like an MD)
-           If `True` `last` is ignored
+        skip_final:
+            the static stress tensor is duplicated in the output when running
+            MD simulations. This flag is used in case one wish to get the final
+            one.
 
         Returns
         -------
@@ -435,46 +433,42 @@ def read_stress(self, key="static", last=True, all=False):
             returns ``None`` if the stresses are not found in the
             output, otherwise stresses will be returned
         """
-        if all:
-            last = False
 
-        # Read until stress are found
-        def next_stress():
-            found, line = self.step_to(f"siesta: Stress tensor", allow_reread=False)
+        is_voigt = key.lower() == "voigt"
+        if is_voigt:
+            key = "Voigt"
+            search = "Stress tensor Voigt"
+        else:
+            search = "siesta: Stress tensor"
+
+        found = False
+        line = " "
+        while not found and line != "":
+            found, line = self.step_to(search, allow_reread=False)
             found = found and key in line
-            while not found and line != "":
-                found, line = self.step_to(f"siesta: Stress tensor", allow_reread=False)
-                found = found and key in line
-            if not found:
-                return None
 
-            # Now read data
+        if not found:
+            return None
+
+        # Now read data
+        if is_voigt:
+            Svoigt = _a.arrayd([float(x) for x in line.split()[-6:]])
+            S = voigt_matrix(Svoigt, False)
+        else:
             S = []
             for _ in range(3):
                 line = self.readline().split()
                 S.append([float(x) for x in line[-3:]])
+                if skip_final:
+                    if line[0].startswith("siesta:"):
+                        return None
+            S = _a.arrayd(S)
 
-            return _a.arrayd(S)
-
-        if all or last:
-            # list of all stresses
-            Ss = []
-            while True:
-                S = next_stress()
-                if S is None:
-                    break
-                Ss.append(S)
-
-            if last:
-                return Ss[-1]
-            if self.completed() and key == "static":
-                return Ss[:-1]
-            return Ss
+        return S
 
-        return next_stress()
-
-    @sile_fh_open(True)
-    def read_moment(self, orbitals=False, quantity="S", last=True, all=False):
+    @SileBinder(postprocess=_a.arrayd)
+    @sile_fh_open()
+    def read_moment(self, orbitals=False, quantity="S"):
         """Reads the moments from the Siesta output file
 
         These will only be present in case of spin-orbit coupling.
@@ -486,13 +480,7 @@ def read_moment(self, orbitals=False, quantity="S", last=True, all=False):
            moments.
         quantity: {'S', 'L'}, optional
            return the spin-moments or the L moments
-        last: bool, optional
-           only read the last force
-        all: bool, optional
-           return a list of all forces (like an MD)
-           If `True` `last` is ignored
         """
-
         # Read until outcoor is found
         if not self.step_to("moments: Atomic", allow_reread=False)[0]:
             return None
@@ -509,6 +497,7 @@ def read_moment(self, orbitals=False, quantity="S", last=True, all=False):
         while True:
             next(itt)  # ""
             next(itt)  # Atom    Orb ...
+
             # Loop atoms in this species list
             while True:
                 line = next(itt)
@@ -523,12 +512,14 @@ def read_moment(self, orbitals=False, quantity="S", last=True, all=False):
                         ia = int(line[0])
                     elif ia != int(line[0]):
                         raise ValueError("Error in moments formatting.")
+
                     # Track maximum number of atoms
                     na = max(ia, na)
                     if quantity == "S":
                         atom.append([float(x) for x in line[4:7]])
                     elif quantity == "L":
                         atom.append([float(x) for x in line[7:10]])
+
                 line = next(itt).split()  # Total ...
                 if not orbitals:
                     ia = int(line[0])
@@ -547,9 +538,7 @@ def read_moment(self, orbitals=False, quantity="S", last=True, all=False):
         for ia, atom in tbl:
             moments[ia - 1] = atom
 
-        if not all:
-            return _a.arrayd(moments)
-        return moments
+        return _a.arrayd(moments)
 
     @sile_fh_open(True)
     def read_energy(self):
@@ -709,9 +698,15 @@ def read_data(self, *args, **kwargs):
         for name in run:
             kwargs.pop(name)
 
+        slice = kwargs.pop("slice", None)
         val = []
         for name in run:
-            val.append(getattr(self, f"read_{name.lower()}")(*args, **kwargs))
+            if slice is None:
+                val.append(getattr(self, f"read_{name.lower()}")(*args, **kwargs))
+            else:
+                val.append(
+                    getattr(self, f"read_{name.lower()}")[slice](*args, **kwargs)
+                )
 
         if len(val) == 0:
             return None
@@ -719,9 +714,16 @@ def read_data(self, *args, **kwargs):
             val = val[0]
         return val
 
-    @sile_fh_open(True)
+    @SileBinder(
+        default_slice=-1, check_empty=_read_scf_empty, postprocess=_read_scf_md_process
+    )
+    @sile_fh_open()
     def read_scf(
-        self, key="scf", iscf=-1, imd=None, as_dataframe=False, ret_header=False
+        self,
+        key: str = "scf",
+        iscf: Optional[int] = -1,
+        as_dataframe: bool = False,
+        ret_header: bool = False,
     ):
         r"""Parse SCF information and return a table of SCF information depending on what is requested
 
@@ -729,18 +731,15 @@ def read_scf(
         ----------
         key : {'scf', 'ts-scf'}
             parse SCF information from Siesta SCF or TranSiesta SCF
-        iscf : int, optional
+        iscf :
             which SCF cycle should be stored. If ``-1`` only the final SCF step is stored,
             for None *all* SCF cycles are returned. When `iscf` values queried are not found they
             will be truncated to the nearest SCF step.
-        imd: int or None, optional
-            whether only a particular MD step is queried, if None, all MD steps are
-            parsed and returned. A negative number wraps for the last MD steps.
-        as_dataframe: boolean, optional
+        as_dataframe:
             whether the information should be returned as a `pandas.DataFrame`. The advantage of this
             format is that everything is indexed and therefore you know what each value means.You can also
             perform operations very easily on a dataframe.
-        ret_header: bool, optional
+        ret_header:
             whether to also return the headers that define each value in the returned array,
             will have no effect if `as_dataframe` is true.
         """
@@ -753,11 +752,6 @@ def read_scf(
                 raise ValueError(
                     f"{self.__class__.__name__}.read_scf requires iscf argument to *not* be 0!"
                 )
-        if not imd is None:
-            if imd == 0:
-                raise ValueError(
-                    f"{self.__class__.__name__}.read_scf requires imd argument to *not* be 0!"
-                )
 
         def reset_d(d, line):
             if line.startswith("SCF cycle converged") or line.startswith(
@@ -913,7 +907,6 @@ def construct_data(d, data):
                 data.extend(d[key])
             d["data"] = data
 
-        md = []
         scf = []
         for line in self:
             parse_next(line, d)
@@ -925,6 +918,7 @@ def construct_data(d, data):
 
                 if iscf is None or iscf < 0:
                     scf.append(data)
+
                 elif data[0] <= iscf:
                     # this ensures we will retain the latest iscf in
                     # case the requested iscf is too big
@@ -956,80 +950,25 @@ def construct_data(d, data):
                     # truncate to 0
                     scf = scf[max(len(scf) + iscf, 0)]
 
-                # Populate md
-                md.append(np.array(scf))
-                # Reset SCF data
-                scf = []
-
-                # In case we wanted a given MD step and it's this one, just stop reading
-                # We are going to return the last MD (see below)
-                if imd == len(md):
-                    break
+                # found a full MD
+                break
 
         # Define the function that is going to convert the information of a MDstep to a Dataset
         if as_dataframe:
             import pandas as pd
 
-            def MDstep_dataframe(scf):
-                scf = np.atleast_2d(scf)
-                return pd.DataFrame(
-                    scf[..., 1:],
-                    index=pd.Index(scf[..., 0].ravel().astype(np.int32), name="iscf"),
-                    columns=props[1:],
-                )
-
-        # Now we know how many MD steps there are
-
-        # We will return stuff based on what the user requested
-        # For pandas DataFrame this will be dependent
-        #  1. all MD steps requested => imd == index, iscf == column (regardless of iscf==none|int)
-        #  2. 1 MD step requested => iscf == index
-
-        if imd is None:
-            if as_dataframe:
-                if len(md) == 0:
-                    # return an empty dataframe (with imd as index)
-                    return pd.DataFrame(index=pd.Index([], name="imd"), columns=props)
-                # Regardless of what the user requests we will always have imd == index
-                # and iscf a column, a user may easily change this.
-                df = pd.concat(
-                    map(MDstep_dataframe, md),
-                    keys=_a.arangei(1, len(md) + 1),
-                    names=["imd"],
-                )
-                if iscf is not None:
-                    df.reset_index("iscf", inplace=True)
-                return df
-
-            if iscf is not None:
-                # since each MD step may be a different number of SCF steps
-                # we can only convert for a specific entry
-                md = np.array(md)
-            if ret_header:
-                return md, props
-            return md
-
-        # correct imd to ensure we check against the final size
-        imd = min(len(md) - 1, max(len(md) + imd, 0))
-        if len(md) == 0:
-            # no data collected
-            if as_dataframe:
+            if len(scf) == 0:
                 return pd.DataFrame(index=pd.Index([], name="iscf"), columns=props[1:])
-            md = np.array(md[imd])
-            if ret_header:
-                return md, props
-            return md
 
-        if imd > len(md):
-            raise ValueError(
-                f"{self.__class__.__name__}.read_scf could not find requested MD step ({imd})."
+            scf = np.atleast_2d(scf)
+            return pd.DataFrame(
+                scf[..., 1:],
+                index=pd.Index(scf[..., 0].ravel().astype(np.int32), name="iscf"),
+                columns=props[1:],
             )
 
-        # If a certain imd was requested, get it
-        # Remember that if imd is positive, we stopped reading at the moment we reached it
-        scf = np.array(md[imd])
-        if as_dataframe:
-            return MDstep_dataframe(scf)
+        # Convert to numpy array
+        scf = np.array(scf)
         if ret_header:
             return scf, props
         return scf
diff --git a/src/sisl/io/siesta/struct.py b/src/sisl/io/siesta/struct.py
index fdca6f0888..a1b82fee7a 100644
--- a/src/sisl/io/siesta/struct.py
+++ b/src/sisl/io/siesta/struct.py
@@ -21,15 +21,15 @@ class structSileSiesta(SileSiesta):
     """Geometry file"""
 
     @sile_fh_open()
-    def write_geometry(self, geometry, fmt=".9f"):
+    def write_geometry(self, geometry: Geometry, fmt: str = ".9f"):
         """Writes the geometry to the contained file
 
         Parameters
         ----------
-        geometry : Geometry
-           geometry to write in the XV file
-        fmt : str, optional
-           the precision used for writing the XV file
+        geometry :
+           geometry to write in the STRUCT file
+        fmt :
+           the precision used for writing the coordinates in the file
         """
         # Check that we can write to the file
         sile_raise_write(self)
@@ -51,7 +51,7 @@ def write_geometry(self, geometry, fmt=".9f"):
                 self._write(fmt_str.format(ips + 1, a.Z, *fxyz[ia]))
 
     @sile_fh_open()
-    def read_lattice(self):
+    def read_lattice(self) -> Lattice:
         """Returns `Lattice` object from the STRUCT file"""
 
         cell = np.empty([3, 3], np.float64)
@@ -61,7 +61,7 @@ def read_lattice(self):
         return Lattice(cell)
 
     @sile_fh_open()
-    def read_geometry(self, species_Z=False):
+    def read_geometry(self, species_Z: bool = False) -> Geometry:
         """Returns a `Geometry` object from the STRUCT file
 
         Parameters
diff --git a/src/sisl/io/siesta/tests/test_stdout.py b/src/sisl/io/siesta/tests/test_stdout.py
index 00b21103dc..61912196f7 100644
--- a/src/sisl/io/siesta/tests/test_stdout.py
+++ b/src/sisl/io/siesta/tests/test_stdout.py
@@ -19,10 +19,9 @@ def test_md_nose_out(sisl_files):
     f = sisl_files(_dir, "md_nose.out")
     out = stdoutSileSiesta(f)
 
-    # nspin, nk, nb
-    geom0 = out.read_geometry(last=False)
-    geom = out.read_geometry()
-    geom1 = out.read_data(geometry=True)
+    geom0 = out.read_geometry()
+    geom = out.read_geometry[-1]()
+    geom1 = out.read_data(geometry=True, slice=-1)
 
     # assert it works correct
     assert isinstance(geom0, sisl.Geometry)
@@ -35,20 +34,20 @@ def test_md_nose_out(sisl_files):
     # try and read all outputs
     # there are 5 outputs in this output file.
     nOutputs = 5
-    assert len(out.read_geometry(all=True)) == nOutputs
-    assert len(out.read_force(all=True)) == nOutputs
-    assert len(out.read_stress(all=True)) == nOutputs
-    f0 = out.read_force(last=False)
-    f = out.read_force()
-    f1 = out.read_data(force=True)
+    assert len(out.read_geometry[:]()) == nOutputs
+    assert len(out.read_force[:]()) == nOutputs
+    assert len(out.read_stress[:]()) == nOutputs
+    f0 = out.read_force()
+    f = out.read_force[-1]()
+    f1 = out.read_data(force=True, slice=-1)
     assert not np.allclose(f0, f)
     assert np.allclose(f1, f)
 
     # Check that we can read the different types of forces
     nAtoms = 10
-    atomicF = out.read_force(all=True)
-    totalF = out.read_force(all=True, total=True)
-    maxF = out.read_force(all=True, max=True)
+    atomicF = out.read_force[:]()
+    totalF = out.read_force[:](total=True)
+    maxF = out.read_force[:](max=True)
     assert atomicF.shape == (nOutputs, nAtoms, 3)
     assert totalF.shape == (nOutputs, 3)
     assert maxF.shape == (nOutputs,)
@@ -56,31 +55,39 @@ def test_md_nose_out(sisl_files):
     assert totalF.shape == (3,)
     assert maxF.shape == ()
 
-    s0 = out.read_stress(last=False)
-    s = out.read_stress()
+    s0 = out.read_stress()
+    s = out.read_stress[-1]()
     assert not np.allclose(s0, s)
 
-    sstatic = out.read_stress("static", all=True)
-    stotal = out.read_stress("total", all=True)
-    sdata = out.read_data("total", all=True, stress=True)
+    sstatic = out.read_stress[:]("static")
+    stotal = out.read_stress[:]("total")
 
-    for S, T, D in zip(sstatic, stotal, sdata):
+    for S, T in zip(sstatic, stotal):
         assert not np.allclose(S, T)
-        assert np.allclose(D, T)
+
+
+def test_md_nose_out_scf(sisl_files):
+    f = sisl_files(_dir, "md_nose.out")
+    out = stdoutSileSiesta(f)
 
     # Ensure SCF reads are consistent
-    scf_last = out.read_scf()
-    scf = out.read_scf(imd=-1)
+    scf_last = out.read_scf[:]()
+    scf_last2, props = out.read_scf[:](ret_header=True)
+    scf = out.read_scf[-1]()
+    scf_props = out.read_scf[-1](ret_header=True)
+    assert scf_props[1] == props
+    assert np.allclose(scf, scf_props[0])
     assert np.allclose(scf_last[-1], scf)
     for i in range(len(scf_last)):
-        scf = out.read_scf(imd=i + 1)
+        scf = out.read_scf[i]()
         assert np.allclose(scf_last[i], scf)
+        assert np.allclose(scf_last[i], scf_last2[i])
 
-    scf_all = out.read_scf(iscf=None, imd=-1)
-    scf = out.read_scf(imd=-1)
+    scf_all = out.read_scf[-1](iscf=None)
+    scf = out.read_scf[-1]()
     assert np.allclose(scf_all[-1], scf)
     for i in range(len(scf_all)):
-        scf = out.read_scf(iscf=i + 1, imd=-1)
+        scf = out.read_scf[-1](iscf=i + 1)
         assert np.allclose(scf_all[i], scf)
 
 
@@ -111,15 +118,15 @@ def test_md_nose_out_dataframe(sisl_files):
     f = sisl_files(_dir, "md_nose.out")
     out = stdoutSileSiesta(f)
 
-    data = out.read_scf()
-    df = out.read_scf(as_dataframe=True)
+    data = out.read_scf[:]()
+    df = out.read_scf[:](as_dataframe=True)
     # this will read all MD-steps and only latest iscf
     assert len(data) == len(df)
     assert df.index.names == ["imd"]
 
-    df = out.read_scf(iscf=None, as_dataframe=True)
+    df = out.read_scf[:](iscf=None, as_dataframe=True)
     assert df.index.names == ["imd", "iscf"]
-    df = out.read_scf(iscf=None, imd=-1, as_dataframe=True)
+    df = out.read_scf(iscf=None, as_dataframe=True)
     assert df.index.names == ["iscf"]
 
 
diff --git a/src/sisl/io/siesta/xv.py b/src/sisl/io/siesta/xv.py
index bec19e73f6..725e468519 100644
--- a/src/sisl/io/siesta/xv.py
+++ b/src/sisl/io/siesta/xv.py
@@ -21,14 +21,14 @@ class xvSileSiesta(SileSiesta):
     """Geometry file"""
 
     @sile_fh_open()
-    def write_geometry(self, geometry, fmt=".9f", velocity=None):
+    def write_geometry(self, geometry: Geometry, fmt: str = ".9f", velocity=None):
         """Writes the geometry to the contained file
 
         Parameters
         ----------
-        geometry : Geometry
+        geometry :
            geometry to write in the XV file
-        fmt : str, optional
+        fmt :
            the precision used for writing the XV file
         velocity : numpy.ndarray, optional
            velocities to write in the XV file (will be zero if not specified).
@@ -68,7 +68,7 @@ def write_geometry(self, geometry, fmt=".9f", velocity=None):
                 self._write(fmt_str.format(ips + 1, a.Z, *tmp))
 
     @sile_fh_open()
-    def read_lattice(self):
+    def read_lattice(self) -> Lattice:
         """Returns `Lattice` object from the XV file"""
 
         cell = np.empty([3, 3], np.float64)
@@ -79,15 +79,17 @@ def read_lattice(self):
         return Lattice(cell)
 
     @sile_fh_open()
-    def read_geometry(self, velocity=False, species_Z=False):
+    def read_geometry(
+        self, velocity: bool = False, species_Z: bool = False
+    ) -> Geometry:
         """Returns a `Geometry` object from the XV file
 
         Parameters
         ----------
-        species_Z : bool, optional
+        species_Z :
            if ``True`` the atomic numbers are the species indices (useful when
            reading the ChemicalSpeciesLabel block simultaneously).
-        velocity : bool, optional
+        velocity :
            also return the velocities in the file
 
         Returns
diff --git a/src/sisl/io/tests/test_object.py b/src/sisl/io/tests/test_object.py
index 99630da557..16075bff11 100644
--- a/src/sisl/io/tests/test_object.py
+++ b/src/sisl/io/tests/test_object.py
@@ -22,6 +22,7 @@
 from sisl.io.siesta.binaries import _gfSileSiesta
 from sisl.io.tbtrans._cdf import *
 from sisl.io.vasp import chgSileVASP
+from sisl.io.xsf import xsfSile
 
 pytestmark = pytest.mark.io
 _dir = osp.join("sisl", "io")
@@ -275,7 +276,7 @@ def test_read_write_lattice(self, sisl_tmp, sisl_system, sile):
         L.set_nsc([1, 1, 1])
         f = sisl_tmp("test_read_write_geom.win", _dir)
         # These files does not store the atomic species
-        if issubclass(sile, (_ncSileTBtrans, deltancSileTBtrans)):
+        if issubclass(sile, (xsfSile, _ncSileTBtrans, deltancSileTBtrans)):
             return
         if sys.platform.startswith("win") and issubclass(sile, chgSileVASP):
             pytest.xfail(
@@ -308,7 +309,7 @@ def test_read_write_lattice(self, sisl_tmp, sisl_system, sile):
         L.set_nsc([1, 1, 1])
         f = sisl_tmp("test_read_write_geom.win", _dir)
         # These files does not store the atomic species
-        if issubclass(sile, (_ncSileTBtrans, deltancSileTBtrans)):
+        if issubclass(sile, (xsfSile, _ncSileTBtrans, deltancSileTBtrans)):
             return
         if sys.platform.startswith("win") and issubclass(sile, chgSileVASP):
             pytest.xfail(
diff --git a/src/sisl/io/xsf.py b/src/sisl/io/xsf.py
index 0e8bfe5c30..5d46e53fd9 100644
--- a/src/sisl/io/xsf.py
+++ b/src/sisl/io/xsf.py
@@ -3,6 +3,7 @@
 # file, You can obtain one at https://mozilla.org/MPL/2.0/.
 import os.path as osp
 from numbers import Integral
+from typing import Optional
 
 import numpy as np
 
@@ -12,7 +13,7 @@
 from sisl.messages import deprecate_argument
 from sisl.utils import str_spec
 
-from ._multiple import SileBinder
+from ._multiple import SileBinder, postprocess_tuple
 
 # Import sile objects
 from .sile import *
@@ -20,7 +21,7 @@
 __all__ = ["xsfSile"]
 
 
-def _get_kw_index(key):
+def _get_kw_index(key: str):
     # Get the integer in a line like 'ATOMS 2', converted to 0-indexing, and with -1 if no int is there
     kl = key.split()
     if len(kl) == 1:
@@ -28,7 +29,7 @@ def _get_kw_index(key):
     return int(kl[1]) - 1
 
 
-def reset_values(*names_values, animsteps=False):
+def reset_values(*names_values, animsteps: bool = False):
     if animsteps:
 
         def reset(self):
@@ -47,19 +48,6 @@ def reset(self):
     return reset
 
 
-def postprocess(*funcs):
-    """Post-processes the returned value according to the funcs for multiple data"""
-
-    def post(ret):
-        nonlocal funcs
-        if isinstance(ret[0], tuple):
-            return tuple(func(r) for r, func in zip(zip(*ret), funcs))
-        return funcs[0](ret)
-
-    return post
-
-
-# Implementation notice!
 # The XSF files are compatible with Vesta, but ONLY
 # if there are no empty lines!
 @set_module("sisl.io")
@@ -92,17 +80,17 @@ def _setup(self, *args, **kwargs):
         self._read_type = None
         self._read_cell = None
 
-    def _write_key(self, key):
+    def _write_key(self, key: str):
         self._write(f"{key}\n")
 
-    def _write_key_index(self, key):
+    def _write_key_index(self, key: str):
         # index is 1-based in file
         if self._geometry_max > 1:
             self._write(f"{key} {self._geometry_write + 1}\n")
         else:
             self._write(f"{key}\n")
 
-    def _write_once(self, string):
+    def _write_once(self, string: str):
         if self._geometry_write <= 0:
             self._write(string)
 
@@ -114,14 +102,14 @@ def _write_animsteps(self):
     @deprecate_argument(
         "sc", "lattice", "use lattice= instead of sc=", from_version="0.15"
     )
-    def write_lattice(self, lattice, fmt=".8f"):
+    def write_lattice(self, lattice: Lattice, fmt: str = ".8f"):
         """Writes the supercell to the contained file
 
         Parameters
         ----------
-        lattice : Lattice
+        lattice :
            the supercell to be written
-        fmt : str, optional
+        fmt :
            used format for the precision of the data
         """
         # Check that we can write to the file
@@ -161,14 +149,14 @@ def write_lattice(self, lattice, fmt=".8f"):
         #    self._write(fmt_str.format(*convcell[i, :]))
 
     @sile_fh_open(reset=reset_values(("_geometry_write", 0), animsteps=True))
-    def write_geometry(self, geometry, fmt=".8f", data=None):
+    def write_geometry(self, geometry: Geometry, fmt: str = ".8f", data=None):
         """Writes the geometry to the contained file
 
         Parameters
         ----------
-        geometry : Geometry
+        geometry :
            the geometry to be written
-        fmt : str, optional
+        fmt :
            used format for the precision of the data
         data : (geometry.na, 3), optional
            auxiliary data associated with the geometry to be saved
@@ -202,7 +190,9 @@ def write_geometry(self, geometry, fmt=".8f", data=None):
                 self._write(fmt_str.format(geometry.atoms[ia].Z, *geometry.xyz[ia, :]))
 
     @sile_fh_open()
-    def _r_geometry_next(self, lattice=None, atoms=None, ret_data=False):
+    def _r_geometry_next(
+        self, lattice: Optional[Lattice] = None, atoms=None, ret_data: bool = False
+    ):
         if lattice is None:
             # fetch the prior read cell value
             lattice = self._read_cell
@@ -291,23 +281,27 @@ def next():
 
         typ = self._read_type
         if typ == "CRYSTAL":
-            nsc = [3, 3, 3]
+            bc = ["periodic", "periodic", "periodic"]
         elif typ == "SLAB":
-            nsc = [3, 3, 1]
+            bc = ["periodic", "periodic", "unknown"]
         elif typ == "POLYMER":
-            nsc = [3, 1, 1]
+            bc = ["periodic", "unknown", "unknown"]
+        elif typ == "MOLECULE":
+            bc = ["unknown", "unknown", "unknown"]
         else:
-            nsc = [1, 1, 1]
+            bc = ["unknown", "unknown", "unknown"]
 
         cell = None
 
         if primvec is not None:
-            cell = Lattice(primvec, nsc=nsc)
+            cell = Lattice(primvec, boundary_condition=bc)
         elif lattice is not None:
             cell = lattice
 
         elif typ == "MOLECULE":
-            cell = Lattice(np.diag(xyz.max(0) - xyz.min(0) + 10.0))
+            cell = Lattice(
+                np.diag(xyz.max(0) - xyz.min(0) + 10.0), boundary_condition=bc
+            )
 
         if cell is None:
             raise ValueError(
@@ -332,25 +326,31 @@ def next():
             return geom, _a.arrayd(data)
         return geom
 
-    @SileBinder(postprocess=postprocess(list, list))
+    @SileBinder(postprocess=postprocess_tuple(list))
+    def read_lattice(self) -> Lattice:
+        """Lattice contained in file"""
+        ret = self._r_geometry_next()
+        if ret is None:
+            return ret
+        return ret.lattice
+
+    @SileBinder(postprocess=postprocess_tuple(list))
     @deprecate_argument(
         "sc", "lattice", "use lattice= instead of sc=", from_version="0.15"
     )
-    def read_geometry(self, lattice=None, atoms=None, ret_data=False):
+    def read_geometry(
+        self, lattice: Optional[Lattice] = None, atoms=None, ret_data: bool = False
+    ):
         """Geometry contained in file, and optionally the associated data
 
-        If the file contains more geometries, one can read multiple geometries
-        by using the arguments `start`, `stop` and `step`.
-        The default is to read the first geometry, only.
-
         Parameters
         ----------
-        lattice : Lattice, optional
+        lattice :
             the supercell in case the lattice vectors are not present in the current
             block.
         atoms : Atoms, optional
             atomic species used regardless of the contained atomic species
-        ret_data : bool, optional
+        ret_data :
            in case the the file has auxiliary data, return that as well.
         """
         return self._r_geometry_next(lattice=lattice, atoms=atoms, ret_data=ret_data)
@@ -474,7 +474,7 @@ def ArgumentParser_out(self, p, *args, **kwargs):
 
         Parameters
         ----------
-        p : ``argparse.ArgumentParser``
+        p : `argparse.ArgumentParser`
            the parser which gets amended the additional output options.
         """
         import argparse
diff --git a/src/sisl/io/xyz.py b/src/sisl/io/xyz.py
index 7ba3b6bc8b..6c753f7bf3 100644
--- a/src/sisl/io/xyz.py
+++ b/src/sisl/io/xyz.py
@@ -4,6 +4,8 @@
 """
 Sile object for reading/writing XYZ files
 """
+from typing import Optional
+
 import numpy as np
 
 import sisl._array as _a
@@ -23,7 +25,7 @@
 class xyzSile(Sile):
     """XYZ file object"""
 
-    def _parse_lattice(self, header, xyz, lattice):
+    def _parse_lattice(self, header: str, xyz, lattice: Optional[Lattice]):
         """Internal helper routine for extracting the lattice"""
         if lattice is not None:
             return lattice
@@ -42,6 +44,7 @@ def _parse_lattice(self, header, xyz, lattice):
                     bc.append(BoundaryCondition.PERIODIC)
                 else:
                     bc.append(BoundaryCondition.UNKNOWN)
+
         if "boundary_condition" in header:
             bc = []
             for b in header.pop("boundary_condition").split():
@@ -62,7 +65,9 @@ def _parse_lattice(self, header, xyz, lattice):
         return Lattice(cell, nsc=nsc, origin=origin, boundary_condition=bc)
 
     @sile_fh_open()
-    def write_geometry(self, geometry, fmt=".8f", comment=None):
+    def write_geometry(
+        self, geometry: Geometry, fmt: str = ".8f", comment: Optional[str] = None
+    ):
         """Writes the geometry to the contained file
 
         Parameters
@@ -124,7 +129,7 @@ def _r_geometry_skip(self, *args, **kwargs):
     @deprecate_argument(
         "sc", "lattice", "use lattice= instead of sc=", from_version="0.15"
     )
-    def read_geometry(self, atoms=None, lattice=None):
+    def read_geometry(self, atoms=None, lattice: Optional[Lattice] = None):
         """Returns Geometry object from the XYZ file
 
         Parameters