fixed eigenmode reader, state order changed

Also ensured that atol is used (not ktol).

read_eigen* methods which can't find stuff based
specific k-values now raises LookupError,
not a KeyError, not an IndexError.

Cleaned tests for the vibra method.
Overall cleaning of the class as well.
Added a read_brillouinzone in the vibra output.

Removed sc argument allowance in eigenmode.
This is deprecated for other methods as well.

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

CHANGELOG.md

@@ -9,6 +9,7 @@ we hit release version 1.0.0.
 ## [0.15.0] - YYYY-MM-DD
 
 ### Added
+- added `vectorsSileSiesta` to read vibra eigenmode output
 - added `dihedral` to `Geometry`, #773
 - ability to retain sub-classes through `<class>.new` calls
 - added `Listify` to ensure arguments behaves as *iterables*

src/sisl/io/siesta/binaries.py

@@ -1919,7 +1919,16 @@ def yield_eigenstate(self):
             # The loop in which the generator was used has been broken.
             self._close_wfsx()
 
-    def read_eigenstate(self, k=(0, 0, 0), spin=0, ktol=1e-4) -> EigenstateElectron:
+    @deprecate_argument(
+        "ktol",
+        "atol",
+        "use atol instead of ktol",
+        "0.15",
+        "0.16",
+    )
+    def read_eigenstate(
+        self, k=(0, 0, 0), spin: int = 0, atol: float = 1e-4
+    ) -> EigenstateElectron:
         """Reads a specific eigenstate from the file.
 
         This method iterates over the states until it finds a match. Do not call
@@ -1929,10 +1938,10 @@ def read_eigenstate(self, k=(0, 0, 0), spin=0, ktol=1e-4) -> EigenstateElectron:
         ----------
         k: array-like of shape (3,), optional
             The k point of the state you want to find.
-        spin: integer, optional
+        spin:
             The spin index of the state you want to find. Only meaningful for polarized
             calculations.
-        ktol: float, optional
+        atol:
             The threshold value for considering two k-points the same (i.e. to match
             the query k point with the states k point).
 
@@ -1944,16 +1953,21 @@ def read_eigenstate(self, k=(0, 0, 0), spin=0, ktol=1e-4) -> EigenstateElectron:
         -------
         EigenstateElectron or None:
             If found, the state that was queried.
-            If not found, returns `None`. NOTE this may change to an exception in the future
+
+        Raises
+        ------
+        LookupError :
+            in case the requested k-point can not be found in the file.
         """
         # Iterate over all eigenstates in the file
         for state in self.yield_eigenstate():
-            if state.info.get("spin", 0) == spin and np.allclose(
-                state.info["k"], k, atol=ktol
-            ):
+            info = state.info
+            if info.get("spin", 0) == spin and np.allclose(info["k"], k, atol=atol):
                 # This is the state that the user requested
                 return state
-        return None
+        raise LookupError(
+            f"{self.__class__.__name__}.read_eigenstate could not find k-point: {k!s} eigenstate"
+        )
 
     def read_info(self):
         """Reads the information for all the k points contained in the file

src/sisl/io/siesta/tests/test_vibra.py

@@ -5,6 +5,7 @@
 
 import os.path as osp
 
+import numpy as np
 import pytest
 
 import sisl
@@ -26,15 +27,31 @@ def test_eigenmode_read(sisl_files):
     nlines = len(blines)
     iklast = [sum(int(l.split()[0]) for l in blines[: i + 1]) for i in range(nlines)]
     klast = [list(map(float, l.split()[1:4])) for l in blines]
-    print(iklast, klast)
 
     # yield modes
     nmodes = 0
     nmodes_total = 0
+    ks = []
     for state in vectors.yield_eigenmode():
         nmodes += 1
         nmodes_total += len(state)
+        ks.append(state.info["k"])
 
-    print(nmodes, nmodes_total)
     assert nmodes == nk
     assert nmodes_total == geometry.na * 3 * nk
+
+    bz = vectors.read_brillouinzone()
+    assert len(bz) == nk
+    assert np.allclose(bz.k, ks)
+
+
+def test_eigenmode_values(sisl_files):
+    fdf = sisl.get_sile(sisl_files(_dir, "h_chain_vibra.fdf"))
+    vectors = sisl.io.siesta.vectorsSileSiesta(
+        sisl_files(_dir, "h_chain_vibra.vectors"),
+        geometry=fdf.read_geometry(),
+    )
+
+    mode = vectors.read_eigenmode()
+    assert np.allclose(mode.state[0, 0:3], [-0.7071, -0.2400e-4, -0.2400e-4])
+    assert np.allclose(mode.state[0, 3:7], [0.7071, 0.2400e-4, 0.2400e-4])

src/sisl/io/siesta/vibra.py

@@ -4,13 +4,15 @@
 from __future__ import annotations
 
 import re
+from typing import Iterator
 
 import numpy as np
 
 import sisl._array as _a
 from sisl import Geometry, Lattice
 from sisl._internal import set_module
 from sisl.messages import warn
+from sisl.physics.brillouinzone import BrillouinZone
 from sisl.physics.phonon import EigenmodePhonon
 from sisl.unit.siesta import unit_convert
 
@@ -30,11 +32,10 @@ class vectorsSileSiesta(SileSiesta):
     Parameters
     ----------
     parent : obj, optional
-        a parent may contain a DynamicalMatrix, or Geometry
+        a parent may contain a DynamicalMatrix, Geometry or Lattice
     geometry : Geometry, optional
         a geometry contains a cell with corresponding lattice vectors
         used to convert k [1/Ang] -> [b], and the atomic masses
-
     """
 
     def _setup(self, *args, **kwargs):
@@ -62,7 +63,7 @@ def _setup(self, *args, **kwargs):
         if geometry is not None and lattice is None:
             lattice = geometry.lattice
 
-        lattice = kwargs.get("lattice", kwargs.get("sc", lattice))
+        lattice = kwargs.get("lattice", lattice)
         if lattice is None and geometry is not None:
             raise ValueError(
                 f"{self.__class__.__name__}(geometry=Geometry, lattice=None) is not an allowed argument combination."
@@ -82,7 +83,7 @@ def conv(k):
                 if not np.allclose(k, 0.0):
                     warn(
                         f"{self.__class__.__name__} cannot convert stored k-points from 1/Ang to reduced coordinates. "
-                        "Please ensure parent=DynamicalMatrix, geometry=Geometry, or lattice=Lattice to ensure reduced k."
+                        "Please ensure parent=DynamicalMatrix, geometry=Geometry, or lattice=Lattice to calculate reduced k."
                     )
                 return k / _Bohr2Ang
 
@@ -93,14 +94,15 @@ def conv(k):
 
         self._convert_k = conv
 
-    def _open(self, rewind=False):
+    def _open(self, rewind: bool = False):
         """Open the file
 
         Here we also initialize some variables to keep track of the state of the read.
         """
         super()._open()
         if rewind:
             self.seek(0)
+
         # Here we initialize the variables that will keep track of the state of the read.
         # The process for identification is done on this basis:
         #  _ik is the current (Python) index for the k-point to be read
@@ -110,27 +112,14 @@ def _open(self, rewind=False):
         self._state = -1
         self._ik = 0
 
-    def _setup_parsing(self, close=True):
-        """Gets all the things needed to parse the wfsx file.
-
-        Parameters
-        -----------
-        close: bool, optional
-            Whether the file unit should be closed afterwards.
-        """
-        self._open()
-        # Read the sizes relevant to the file.
-        # We also read whether there's only gamma point information or there are multiple points
-        if self._state == -1:
-            self._sizes = self._r_next_sizes()
-            self._state = 0
-
-        if close:
-            self._close()
+        # read in sizes
+        self._sizes = self._r_next_sizes()
+        self._state = 0
 
-    def _r_next_sizes(self):
+    def _r_next_sizes(self) -> None:
         """Determine the dimension of the data stored in the vectors file"""
         pos = self.fh.tell()
+
         # Skip past header
         self.readline()  # Empty line
         self.readline()  # K point
@@ -161,7 +150,7 @@ def _r_next_sizes(self):
         # Rewind
         self.seek(pos)
 
-    def _r_next_eigenmode(self):
+    def _r_next_eigenmode(self) -> EigenmodePhonon:
         """Reads the next phonon eigenmode in the vectors file.
 
         Returns
@@ -172,14 +161,14 @@ def _r_next_eigenmode(self):
         # Skip empty line at the head of the file
         self.readline()
 
-        k = _a.asarrayd(list(map(float, self.readline().split()[2:])))
+        k = _a.fromiterd(map(float, self.readline().split()[2:]))
         if len(k) == 0:
-            raise GeneratorExit
+            return None
         # Read first eigenvector index
 
         # Determine number of atoms (i.e. rows per mode)
-        mode = np.empty((self._nmodes, 3, self._natoms), dtype=np.complex128)
-        c = np.empty(self._nmodes, dtype=np.float64)
+        mode = _a.emptyz([self._nmodes, self._natoms, 3])
+        c = _a.emptyd(self._nmodes)
         for imode in range(self._nmodes):
             self.readline()
 
@@ -191,39 +180,36 @@ def _r_next_eigenmode(self):
             self.readline()
             for iatom in range(self._natoms):
                 line = self.readline()
-                mode[imode, :, iatom].real = list(map(float, line.split()))
+                mode[imode, iatom].real = list(map(float, line.split()))
 
             # Read imaginary part of eigenmode
             # Skip eigenmode header
             self.readline()
             for iatom in range(self._natoms):
                 line = self.readline()
-                mode[imode, :, iatom].imag = list(map(float, line.split()))
+                mode[imode, iatom].imag = list(map(float, line.split()))
 
         info = dict(k=self._convert_k(k), parent=self._parent, gauge="orbital")
         return EigenmodePhonon(mode.reshape(self._nmodes, -1), c * _cm1_eV, **info)
 
-    def yield_eigenmode(self):
+    def yield_eigenmode(self) -> Iterator[EigenmodePhonon]:
         """Iterates over the modes in the vectors file
 
         Yields
         ------
         EigenmodePhonon
         """
         # Open file and get parsing information
-        self._setup_parsing(close=False)
-
-        try:
-            # Iterate over all eigenmodes in the WFSX file, yielding control to the caller at
-            # each iteration.
-            while True:
-                yield self._r_next_eigenmode()
-            # We ran out of eigenmodes
-        except GeneratorExit:
-            # The loop in which the generator was used has been broken.
-            pass
-
-    def read_eigenmode(self, k=(0, 0, 0), ktol: float = 1e-4) -> EigenmodePhonon:
+        self._open()
+
+        # Iterate over all eigenmodes in the vectors file, yielding control to the caller at
+        # each iteration.
+        mode = self._r_next_eigenmode()
+        while mode is not None:
+            yield mode
+            mode = self._r_next_eigenmode()
+
+    def read_eigenmode(self, k=(0, 0, 0), atol: float = 1e-4) -> EigenmodePhonon:
         """Reads a specific eigenmode from the file.
 
         This method iterates over the modes until it finds a match. Do not call
@@ -233,7 +219,7 @@ def read_eigenmode(self, k=(0, 0, 0), ktol: float = 1e-4) -> EigenmodePhonon:
         ----------
         k: array-like of shape (3,), optional
             The k point of the mode you want to find.
-        ktol:
+        atol:
             The threshold value for considering two k-points the same (i.e. to match
             the query k point with the modes k point).
 
@@ -245,14 +231,25 @@ def read_eigenmode(self, k=(0, 0, 0), ktol: float = 1e-4) -> EigenmodePhonon:
         -------
         EigenmodePhonon or None:
             If found, the mode that was queried.
-            If not found, returns `None`. NOTE this may change to an exception in the future
+
+        Raises
+        ------
+        LookupError :
+            in case the requested k-point can not be found in the file.
         """
         # Iterate over all eigenmodes in the file
         for mode in self.yield_eigenmode():
-            if np.allclose(mode.info["k"], k, atol=ktol):
+            if np.allclose(mode.info["k"], k, atol=atol):
                 # This is the mode that the user requested
                 return mode
-        return None
+        raise LookupError(
+            f"{self.__class__.__name__}.read_eigenmode could not find k-point: {k!s} eigenmode"
+        )
+
+    def read_brillouinzone(self) -> BrillouinZone:
+        """Read the brillouin zone object (only containing the k-points)"""
+        k = [mode.info["k"] for mode in self.yield_eigenmode()]
+        return BrillouinZone(self._parent, k=k)
 
 
 add_sile("vectors", vectorsSileSiesta, gzip=True)