From 15bdcc93a328af87af9aab92be7a99fc26d8d41f Mon Sep 17 00:00:00 2001
From: Nils Wittemeier <nils.wittemeier@gmail.com>
Date: Wed, 10 Apr 2024 10:55:10 +0200
Subject: [PATCH 01/10] draft vectorsSileSIESTA

---
 src/sisl/io/siesta/__init__.py |   1 +
 src/sisl/io/siesta/vibra.py    | 257 +++++++++++++++++++++++++++++++++
 2 files changed, 258 insertions(+)
 create mode 100644 src/sisl/io/siesta/vibra.py

diff --git a/src/sisl/io/siesta/__init__.py b/src/sisl/io/siesta/__init__.py
index 3a1faf404..36940a2d4 100644
--- a/src/sisl/io/siesta/__init__.py
+++ b/src/sisl/io/siesta/__init__.py
@@ -64,4 +64,5 @@
 from .stdout import *
 from .struct import *
 from .transiesta_grid import *
+from .vibra import *
 from .xv import *
diff --git a/src/sisl/io/siesta/vibra.py b/src/sisl/io/siesta/vibra.py
new file mode 100644
index 000000000..da2d74304
--- /dev/null
+++ b/src/sisl/io/siesta/vibra.py
@@ -0,0 +1,257 @@
+# 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 __future__ import annotations
+
+import numpy as np
+from itertools import product
+
+from sisl import Geometry, Lattice, constant
+from sisl._internal import set_module
+from sisl.messages import warn
+from sisl.physics.phonon import EigenmodePhonon
+from sisl.unit.siesta import unit_convert
+
+from ..sile import add_sile, sile_fh_open, sile_raise_write
+from .sile import SileSiesta
+
+
+__all__ = ["vectorsSileSiesta"]
+
+_Bohr2Ang = unit_convert("Bohr", "Ang")
+_J2eV = unit_convert("J", "eV")
+_hc_eV = constant.h * constant.c * _J2eV
+
+@set_module("sisl.io.siesta")
+class vectorsSileSiesta(SileSiesta):
+    """Phonon eigenmode file"""
+
+    def _setup(self, *args, **kwargs):
+        """Simple setup that needs to be overwritten
+
+        All _r_next_* methods expect the fortran file unit to be handled
+        and that the position in the file is correct.
+        """
+        super()._setup(*args, **kwargs)
+
+        # default lattice
+        lattice = None
+
+        # In case the instantiation was called with wfsxSileSiesta("path", geometry=geometry)
+        parent = kwargs.get("parent")
+        if parent is None:
+            geometry = None
+        elif isinstance(parent, Geometry):
+            geometry = parent
+        elif isinstance(parent, Lattice):
+            lattice = parent
+        else:
+            geometry = parent.geometry
+
+        geometry = kwargs.get("geometry", geometry)
+        if geometry is not None and lattice is None:
+            lattice = geometry.lattice
+
+        lattice = kwargs.get("lattice", kwargs.get("sc", 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."
+            )
+
+        if parent is None:
+            parent = geometry
+        if parent is None:
+            parent = lattice
+
+        self._parent = parent
+        self._geometry = geometry
+        self._lattice = lattice
+        if self._parent is None and self._geometry is None and self._lattice is None:
+
+            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=Hamiltonian, geometry=Geometry, or lattice=Lattice to ensure reduced k."
+                    )
+                return k / _Bohr2Ang
+
+        else:
+
+            def conv(k):
+                return (k @ lattice.cell.T) / (2 * np.pi * _Bohr2Ang)
+
+        self._convert_k = conv
+
+    def _open(self, rewind=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
+        #  _state is:
+        #        -1 : the file-descriptor has just been opened (i.e. in front of header)
+        #         0 : it means that the file-descriptor is in front of k point information
+        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_sizes()
+            self._state = 0
+
+        if close:
+            self._close()
+
+    def _r_sizes(self):
+        """ 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
+        self.readline() # Eigenmode index
+        self.readline() # Frequency
+        self.readline() # Eigenmode header
+        
+        # Determine number of atoms (i.e. rows per mode)
+        natoms = 0
+        while True:
+            line = self.readline()
+            if line.startswith('Eigenmode'):
+                break
+            else:
+                natoms += 1
+        self._natoms = natoms
+
+        # Determine number of modes
+        nmodes = 1
+        while True:
+            line = self.readline()
+            if line.startswith('k            ='):
+                break 
+            elif line.startswith('Eigenvector'):
+                nmodes += 1
+        self._nmodes = nmodes
+
+        # Rewind
+        self.seek(pos)
+
+    def _r_next_eigenmode(self):
+        """Reads the next phonon eigenmode in the vectors file.
+
+        Parameters
+        ----------
+        ik: integer
+            the (python) k index of the next eigenstate.
+
+        Returns
+        --------
+        EigenmodePhonon:
+            The next eigenmode.
+        """
+        # Skip empty line at the head of the file
+        self.readline()
+            
+        k = list(map(float, self.readline().split()[2:]))
+        if len(k) == 0:
+            raise GeneratorExit
+        # Read first eigenvector index 
+        
+        # Determine number of atoms (i.e. rows per mode)
+        state = np.empty((self._nmodes, self._natoms, 3), dtype=np.complex128)
+        c = np.empty(self._nmodes, dtype=np.float64)
+        for imode in range(self._nmodes):
+            line = self.readline()
+            if line.strip() == '':
+                # End of block for current k is indicated by empty line
+                break
+            
+            # Read frequency
+            c[imode] = float(self.readline().split('=')[1].strip()) * _hc_eV
+            
+            # Read real part of eigenmode
+            # Skip eigenmode header
+            self.readline()
+            for iatom in range(self._natoms):
+                line = self.readline().strip().lower()
+                state[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().strip().lower()
+                state[imode,iatom,:].imag = list(map(float, line.split()))
+
+        info = dict(k=self._convert_k(k))
+        return EigenmodePhonon(state, c, **info)
+
+    def yield_eigenmode(self):
+        """Iterates over the states in the vectors file
+
+        Yields
+        ------
+        EigenmodePhonon
+        """
+        # Open file and get parsing information
+        self._setup_parsing(close=False)
+
+        try:
+            # Iterate over all eigenstates in the WFSX file, yielding control to the caller at
+            # each iteration.
+            while True:
+                yield self._r_next_eigenmode()
+            # We ran out of eigenstates
+        except GeneratorExit:
+            # The loop in which the generator was used has been broken.
+            pass
+
+    def read_eigenmode(self, k=(0,0,0), ktol=1e-4) -> EigenmodePhonon:
+        """Reads a specific eigenstate from the file.
+
+        This method iterates over the modes until it finds a match. Do not call
+        this method repeatedly. If you want to loop eigenstates, use `yield_eigenstate`.
+
+        Parameters
+        ----------
+        k: array-like of shape (3,), optional
+            The k point of the state you want to find.
+        ktol: float, optional
+            The threshold value for considering two k-points the same (i.e. to match
+            the query k point with the states k point).
+
+        See Also
+        --------
+        yield_eigenstate
+
+        Returns
+        -------
+        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
+        """
+        # Iterate over all eigenstates in the file
+        for state in self.yield_eigenstate():
+            if np.allclose(state.info["k"], k, atol=ktol):
+                # This is the state that the user requested
+                return state
+        return None
+
+
+add_sile("vectors", vectorsSileSiesta, gzip=True)

From 3a4269fe69733e8992977e9842f2ac0a1be6096c Mon Sep 17 00:00:00 2001
From: Nils Wittemeier <nils.wittemeier@gmail.com>
Date: Wed, 10 Apr 2024 14:51:33 +0200
Subject: [PATCH 02/10] minor changes according to first review

---
 src/sisl/io/siesta/vibra.py | 84 +++++++++++++++++++------------------
 1 file changed, 43 insertions(+), 41 deletions(-)

diff --git a/src/sisl/io/siesta/vibra.py b/src/sisl/io/siesta/vibra.py
index da2d74304..07a67a423 100644
--- a/src/sisl/io/siesta/vibra.py
+++ b/src/sisl/io/siesta/vibra.py
@@ -4,7 +4,6 @@
 from __future__ import annotations
 
 import numpy as np
-from itertools import product
 
 from sisl import Geometry, Lattice, constant
 from sisl._internal import set_module
@@ -12,7 +11,7 @@
 from sisl.physics.phonon import EigenmodePhonon
 from sisl.unit.siesta import unit_convert
 
-from ..sile import add_sile, sile_fh_open, sile_raise_write
+from ..sile import add_sile
 from .sile import SileSiesta
 
 
@@ -22,14 +21,27 @@
 _J2eV = unit_convert("J", "eV")
 _hc_eV = constant.h * constant.c * _J2eV
 
+
 @set_module("sisl.io.siesta")
 class vectorsSileSiesta(SileSiesta):
-    """Phonon eigenmode file"""
+    """Phonon eigenmode file
+
+    Parameters
+    ----------
+    parent : obj, optional
+        a parent may contain a geometry, and/or a supercell
+    geometry : Geometry, optional
+        a geometry contains a cell with corresponding lattice vectors
+        used to convert k [1/Ang] -> [b]
+    lattice : Lattice, optional
+        a supercell contains the lattice vectors to convert k
+
+    """
 
     def _setup(self, *args, **kwargs):
         """Simple setup that needs to be overwritten
 
-        All _r_next_* methods expect the fortran file unit to be handled
+        All _r_next_* methods expect the file unit to be handled
         and that the position in the file is correct.
         """
         super()._setup(*args, **kwargs)
@@ -37,7 +49,6 @@ def _setup(self, *args, **kwargs):
         # default lattice
         lattice = None
 
-        # In case the instantiation was called with wfsxSileSiesta("path", geometry=geometry)
         parent = kwargs.get("parent")
         if parent is None:
             geometry = None
@@ -112,28 +123,27 @@ def _setup_parsing(self, close=True):
         # 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_sizes()
+            self._sizes = self._r_next_sizes()
             self._state = 0
 
         if close:
             self._close()
 
-    def _r_sizes(self):
-        """ Determine the dimension of the data stored in the vectors file
-        """
+    def _r_next_sizes(self):
+        """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
-        self.readline() # Eigenmode index
-        self.readline() # Frequency
-        self.readline() # Eigenmode header
-        
+        self.readline()  # Empty line
+        self.readline()  # K point
+        self.readline()  # Eigenmode index
+        self.readline()  # Frequency
+        self.readline()  # Eigenmode header
+
         # Determine number of atoms (i.e. rows per mode)
         natoms = 0
         while True:
             line = self.readline()
-            if line.startswith('Eigenmode'):
+            if line.startswith("Eigenmode"):
                 break
             else:
                 natoms += 1
@@ -143,9 +153,9 @@ def _r_sizes(self):
         nmodes = 1
         while True:
             line = self.readline()
-            if line.startswith('k            ='):
-                break 
-            elif line.startswith('Eigenvector'):
+            if line.startswith("k            ="):
+                break
+            elif line.startswith("Eigenvector"):
                 nmodes += 1
         self._nmodes = nmodes
 
@@ -155,11 +165,6 @@ def _r_sizes(self):
     def _r_next_eigenmode(self):
         """Reads the next phonon eigenmode in the vectors file.
 
-        Parameters
-        ----------
-        ik: integer
-            the (python) k index of the next eigenstate.
-
         Returns
         --------
         EigenmodePhonon:
@@ -167,40 +172,37 @@ def _r_next_eigenmode(self):
         """
         # Skip empty line at the head of the file
         self.readline()
-            
+
         k = list(map(float, self.readline().split()[2:]))
         if len(k) == 0:
             raise GeneratorExit
-        # Read first eigenvector index 
-        
+        # Read first eigenvector index
+
         # Determine number of atoms (i.e. rows per mode)
         state = np.empty((self._nmodes, self._natoms, 3), dtype=np.complex128)
         c = np.empty(self._nmodes, dtype=np.float64)
         for imode in range(self._nmodes):
-            line = self.readline()
-            if line.strip() == '':
-                # End of block for current k is indicated by empty line
-                break
-            
+            self.readline()
+
             # Read frequency
-            c[imode] = float(self.readline().split('=')[1].strip()) * _hc_eV
-            
+            c[imode] = float(self.readline().split("=")[1])
+
             # Read real part of eigenmode
             # Skip eigenmode header
             self.readline()
             for iatom in range(self._natoms):
-                line = self.readline().strip().lower()
-                state[imode,iatom,:].real = list(map(float, line.split()))
+                line = self.readline()
+                state[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().strip().lower()
-                state[imode,iatom,:].imag = list(map(float, line.split()))
+                line = self.readline()
+                state[imode, iatom, :].imag = list(map(float, line.split()))
 
-        info = dict(k=self._convert_k(k))
-        return EigenmodePhonon(state, c, **info)
+        info = dict(k=self._convert_k(k), parent=self.parent)
+        return EigenmodePhonon(state, c * _hc_eV, **info)
 
     def yield_eigenmode(self):
         """Iterates over the states in the vectors file
@@ -222,7 +224,7 @@ def yield_eigenmode(self):
             # The loop in which the generator was used has been broken.
             pass
 
-    def read_eigenmode(self, k=(0,0,0), ktol=1e-4) -> EigenmodePhonon:
+    def read_eigenmode(self, k=(0, 0, 0), ktol=1e-4) -> EigenmodePhonon:
         """Reads a specific eigenstate from the file.
 
         This method iterates over the modes until it finds a match. Do not call

From 132a25c072c867c08d9241b074e6629f5c9ba1b8 Mon Sep 17 00:00:00 2001
From: Nils Wittemeier <nils.wittemeier@gmail.com>
Date: Thu, 11 Apr 2024 09:55:48 +0200
Subject: [PATCH 03/10] minor updates

---
 src/sisl/io/siesta/vibra.py | 44 ++++++++++++++++++-------------------
 1 file changed, 21 insertions(+), 23 deletions(-)

diff --git a/src/sisl/io/siesta/vibra.py b/src/sisl/io/siesta/vibra.py
index 07a67a423..96d866b8d 100644
--- a/src/sisl/io/siesta/vibra.py
+++ b/src/sisl/io/siesta/vibra.py
@@ -5,6 +5,7 @@
 
 import numpy as np
 
+import sisl._array as _a
 from sisl import Geometry, Lattice, constant
 from sisl._internal import set_module
 from sisl.messages import warn
@@ -14,12 +15,11 @@
 from ..sile import add_sile
 from .sile import SileSiesta
 
-
 __all__ = ["vectorsSileSiesta"]
 
 _Bohr2Ang = unit_convert("Bohr", "Ang")
 _J2eV = unit_convert("J", "eV")
-_hc_eV = constant.h * constant.c * _J2eV
+_cm1_eV = constant.h * constant.c * _J2eV
 
 
 @set_module("sisl.io.siesta")
@@ -29,12 +29,10 @@ class vectorsSileSiesta(SileSiesta):
     Parameters
     ----------
     parent : obj, optional
-        a parent may contain a geometry, and/or a supercell
+        a parent may contain a DynamicalMatrix, or Geometry
     geometry : Geometry, optional
         a geometry contains a cell with corresponding lattice vectors
-        used to convert k [1/Ang] -> [b]
-    lattice : Lattice, optional
-        a supercell contains the lattice vectors to convert k
+        used to convert k [1/Ang] -> [b], and the atomic masses
 
     """
 
@@ -83,7 +81,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=Hamiltonian, geometry=Geometry, or lattice=Lattice to ensure reduced k."
+                        "Please ensure parent=DynamicalMatrix, geometry=Geometry, or lattice=Lattice to ensure reduced k."
                     )
                 return k / _Bohr2Ang
 
@@ -173,13 +171,13 @@ def _r_next_eigenmode(self):
         # Skip empty line at the head of the file
         self.readline()
 
-        k = list(map(float, self.readline().split()[2:]))
+        k = _a.asarrayd(list(map(float, self.readline().split()[2:])))
         if len(k) == 0:
             raise GeneratorExit
         # Read first eigenvector index
 
         # Determine number of atoms (i.e. rows per mode)
-        state = np.empty((self._nmodes, self._natoms, 3), dtype=np.complex128)
+        state = np.empty((self._nmodes, 3, self._natoms), dtype=np.complex128)
         c = np.empty(self._nmodes, dtype=np.float64)
         for imode in range(self._nmodes):
             self.readline()
@@ -192,17 +190,17 @@ def _r_next_eigenmode(self):
             self.readline()
             for iatom in range(self._natoms):
                 line = self.readline()
-                state[imode, iatom, :].real = list(map(float, line.split()))
+                state[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()
-                state[imode, iatom, :].imag = list(map(float, line.split()))
+                state[imode, :, iatom].imag = list(map(float, line.split()))
 
-        info = dict(k=self._convert_k(k), parent=self.parent)
-        return EigenmodePhonon(state, c * _hc_eV, **info)
+        info = dict(k=self._convert_k(k), parent=self._parent, gauge="r")
+        return EigenmodePhonon(state.reshape(self._nmodes, -1), c * _cm1_eV, **info)
 
     def yield_eigenmode(self):
         """Iterates over the states in the vectors file
@@ -215,41 +213,41 @@ def yield_eigenmode(self):
         self._setup_parsing(close=False)
 
         try:
-            # Iterate over all eigenstates in the WFSX file, yielding control to the caller at
+            # 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 eigenstates
+            # 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=1e-4) -> EigenmodePhonon:
-        """Reads a specific eigenstate from the file.
+    def read_eigenmode(self, k=(0, 0, 0), ktol: 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
-        this method repeatedly. If you want to loop eigenstates, use `yield_eigenstate`.
+        this method repeatedly. If you want to loop eigenmodes, use `yield_eigenmode`.
 
         Parameters
         ----------
         k: array-like of shape (3,), optional
             The k point of the state you want to find.
-        ktol: float, optional
+        ktol:
             The threshold value for considering two k-points the same (i.e. to match
             the query k point with the states k point).
 
         See Also
         --------
-        yield_eigenstate
+        yield_eigenmode
 
         Returns
         -------
-        EigenstateElectron or None:
+        eigenmodeElectron or None:
             If found, the state that was queried.
             If not found, returns `None`. NOTE this may change to an exception in the future
         """
-        # Iterate over all eigenstates in the file
-        for state in self.yield_eigenstate():
+        # Iterate over all eigenmodes in the file
+        for state in self.yield_eigenmode():
             if np.allclose(state.info["k"], k, atol=ktol):
                 # This is the state that the user requested
                 return state

From 5a416daf3fc2726f648626b42adc12065cbbed18 Mon Sep 17 00:00:00 2001
From: Nils Wittemeier <nils.wittemeier@gmail.com>
Date: Wed, 22 May 2024 12:34:38 +0200
Subject: [PATCH 04/10] Make unit conversion easier to read

---
 src/sisl/io/siesta/vibra.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/sisl/io/siesta/vibra.py b/src/sisl/io/siesta/vibra.py
index 96d866b8d..d4ce827c7 100644
--- a/src/sisl/io/siesta/vibra.py
+++ b/src/sisl/io/siesta/vibra.py
@@ -18,8 +18,7 @@
 __all__ = ["vectorsSileSiesta"]
 
 _Bohr2Ang = unit_convert("Bohr", "Ang")
-_J2eV = unit_convert("J", "eV")
-_cm1_eV = constant.h * constant.c * _J2eV
+_cm1_eV = unit_convert("cm^-1", "eV")
 
 
 @set_module("sisl.io.siesta")

From abf3818d9d363343e1e2c8cb7762826d728b1397 Mon Sep 17 00:00:00 2001
From: Nils Wittemeier <nils.wittemeier@gmail.com>
Date: Wed, 22 May 2024 12:35:15 +0200
Subject: [PATCH 05/10] Add test

---
 src/sisl/io/siesta/tests/test_vibra.py | 41 ++++++++++++++++++++++++++
 1 file changed, 41 insertions(+)
 create mode 100644 src/sisl/io/siesta/tests/test_vibra.py

diff --git a/src/sisl/io/siesta/tests/test_vibra.py b/src/sisl/io/siesta/tests/test_vibra.py
new file mode 100644
index 000000000..a6287ce5e
--- /dev/null
+++ b/src/sisl/io/siesta/tests/test_vibra.py
@@ -0,0 +1,41 @@
+# 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 __future__ import annotations
+
+import os.path as osp
+
+import numpy as np
+import pytest
+
+import sisl
+
+pytestmark = [pytest.mark.io, pytest.mark.siesta]
+_dir = osp.join("sisl", "io", "siesta")
+
+
+def test_eigenmode_read(sisl_files):
+    fdf = sisl.get_sile(sisl_files(_dir, "si_vibra.fdf"))
+    geometry = fdf.read_geometry()
+    vectors = sisl.io.siesta.vectorsSileSiesta(
+        sisl_files(_dir, "si_vibra.vectors"),
+        geometry=geometry,
+    )
+
+    blines = fdf.get("BandLines")
+    nk = sum(int(l.split()[0]) for l in blines)
+    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
+    for state in vectors.yield_eigenmode():
+        nmodes += 1
+        nmodes_total += len(state)
+
+    print(nmodes, nmodes_total)
+    assert nmodes == nk
+    assert nmodes_total == geometry.na * 3 * nk

From f7463d2f78bf1f90019f189ddd5cb1468dcd4cf5 Mon Sep 17 00:00:00 2001
From: Nils Wittemeier <nils.wittemeier@gmail.com>
Date: Wed, 22 May 2024 12:43:39 +0200
Subject: [PATCH 06/10] code QL

---
 src/sisl/io/siesta/tests/test_vibra.py | 1 -
 src/sisl/io/siesta/vibra.py            | 2 +-
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/sisl/io/siesta/tests/test_vibra.py b/src/sisl/io/siesta/tests/test_vibra.py
index a6287ce5e..c9092042a 100644
--- a/src/sisl/io/siesta/tests/test_vibra.py
+++ b/src/sisl/io/siesta/tests/test_vibra.py
@@ -5,7 +5,6 @@
 
 import os.path as osp
 
-import numpy as np
 import pytest
 
 import sisl
diff --git a/src/sisl/io/siesta/vibra.py b/src/sisl/io/siesta/vibra.py
index d4ce827c7..685f1fbd0 100644
--- a/src/sisl/io/siesta/vibra.py
+++ b/src/sisl/io/siesta/vibra.py
@@ -6,7 +6,7 @@
 import numpy as np
 
 import sisl._array as _a
-from sisl import Geometry, Lattice, constant
+from sisl import Geometry, Lattice
 from sisl._internal import set_module
 from sisl.messages import warn
 from sisl.physics.phonon import EigenmodePhonon

From e0e183abfc8a3bb5530bf5b04ced69b9cf3ee206 Mon Sep 17 00:00:00 2001
From: Nils Wittemeier <nils.wittemeier@gmail.com>
Date: Wed, 22 May 2024 14:25:50 +0200
Subject: [PATCH 07/10] Remove all mentions of electrons or state in favor of
 modes and phonons

---
 src/sisl/io/siesta/vibra.py | 26 +++++++++++++-------------
 1 file changed, 13 insertions(+), 13 deletions(-)

diff --git a/src/sisl/io/siesta/vibra.py b/src/sisl/io/siesta/vibra.py
index 685f1fbd0..87726e173 100644
--- a/src/sisl/io/siesta/vibra.py
+++ b/src/sisl/io/siesta/vibra.py
@@ -176,7 +176,7 @@ def _r_next_eigenmode(self):
         # Read first eigenvector index
 
         # Determine number of atoms (i.e. rows per mode)
-        state = np.empty((self._nmodes, 3, self._natoms), dtype=np.complex128)
+        mode = np.empty((self._nmodes, 3, self._natoms), dtype=np.complex128)
         c = np.empty(self._nmodes, dtype=np.float64)
         for imode in range(self._nmodes):
             self.readline()
@@ -189,20 +189,20 @@ def _r_next_eigenmode(self):
             self.readline()
             for iatom in range(self._natoms):
                 line = self.readline()
-                state[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()
-                state[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="r")
-        return EigenmodePhonon(state.reshape(self._nmodes, -1), c * _cm1_eV, **info)
+        return EigenmodePhonon(mode.reshape(self._nmodes, -1), c * _cm1_eV, **info)
 
     def yield_eigenmode(self):
-        """Iterates over the states in the vectors file
+        """Iterates over the modes in the vectors file
 
         Yields
         ------
@@ -230,10 +230,10 @@ def read_eigenmode(self, k=(0, 0, 0), ktol: float = 1e-4) -> EigenmodePhonon:
         Parameters
         ----------
         k: array-like of shape (3,), optional
-            The k point of the state you want to find.
+            The k point of the mode you want to find.
         ktol:
             The threshold value for considering two k-points the same (i.e. to match
-            the query k point with the states k point).
+            the query k point with the modes k point).
 
         See Also
         --------
@@ -241,15 +241,15 @@ def read_eigenmode(self, k=(0, 0, 0), ktol: float = 1e-4) -> EigenmodePhonon:
 
         Returns
         -------
-        eigenmodeElectron or None:
-            If found, the state that was queried.
+        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
         """
         # Iterate over all eigenmodes in the file
-        for state in self.yield_eigenmode():
-            if np.allclose(state.info["k"], k, atol=ktol):
-                # This is the state that the user requested
-                return state
+        for mode in self.yield_eigenmode():
+            if np.allclose(mode.info["k"], k, atol=ktol):
+                # This is the mode that the user requested
+                return mode
         return None
 
 

From 5abd747d1a551794cc7d0f57a5b789501a536886 Mon Sep 17 00:00:00 2001
From: Nils Wittemeier <nils.wittemeier@gmail.com>
Date: Wed, 22 May 2024 14:50:02 +0200
Subject: [PATCH 08/10] Update test

---
 src/sisl/io/siesta/tests/test_vibra.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/sisl/io/siesta/tests/test_vibra.py b/src/sisl/io/siesta/tests/test_vibra.py
index c9092042a..fba6da026 100644
--- a/src/sisl/io/siesta/tests/test_vibra.py
+++ b/src/sisl/io/siesta/tests/test_vibra.py
@@ -14,10 +14,10 @@
 
 
 def test_eigenmode_read(sisl_files):
-    fdf = sisl.get_sile(sisl_files(_dir, "si_vibra.fdf"))
+    fdf = sisl.get_sile(sisl_files(_dir, "h_chain_vibra.fdf"))
     geometry = fdf.read_geometry()
     vectors = sisl.io.siesta.vectorsSileSiesta(
-        sisl_files(_dir, "si_vibra.vectors"),
+        sisl_files(_dir, "h_chain_vibra.vectors"),
         geometry=geometry,
     )
 

From 9707908a96ffe8a0f9694086449f82e4abf72cd7 Mon Sep 17 00:00:00 2001
From: Nils Wittemeier <nils.wittemeier@gmail.com>
Date: Wed, 22 May 2024 14:54:46 +0200
Subject: [PATCH 09/10] use regular expression to match 'k *='

---
 src/sisl/io/siesta/vibra.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/sisl/io/siesta/vibra.py b/src/sisl/io/siesta/vibra.py
index 87726e173..942f82ffb 100644
--- a/src/sisl/io/siesta/vibra.py
+++ b/src/sisl/io/siesta/vibra.py
@@ -3,6 +3,8 @@
 # file, You can obtain one at https://mozilla.org/MPL/2.0/.
 from __future__ import annotations
 
+import re
+
 import numpy as np
 
 import sisl._array as _a
@@ -150,7 +152,7 @@ def _r_next_sizes(self):
         nmodes = 1
         while True:
             line = self.readline()
-            if line.startswith("k            ="):
+            if re.match("k *=", line):
                 break
             elif line.startswith("Eigenvector"):
                 nmodes += 1

From 6c94a270fa46ea2c62448e2891d48414dbbde4e0 Mon Sep 17 00:00:00 2001
From: Nils Wittemeier <nils.wittemeier@gmail.com>
Date: Mon, 27 May 2024 12:26:38 +0200
Subject: [PATCH 10/10] Update files submodule

---
 files                       | 2 +-
 src/sisl/io/siesta/vibra.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/files b/files
index b8d5ced4a..cf0032203 160000
--- a/files
+++ b/files
@@ -1 +1 @@
-Subproject commit b8d5ced4aac711a7e00e9939ae0849bc2001af0b
+Subproject commit cf00322030f1c502aff697f9cd1dd136ffa436db
diff --git a/src/sisl/io/siesta/vibra.py b/src/sisl/io/siesta/vibra.py
index 942f82ffb..56fc7fd29 100644
--- a/src/sisl/io/siesta/vibra.py
+++ b/src/sisl/io/siesta/vibra.py
@@ -200,7 +200,7 @@ def _r_next_eigenmode(self):
                 line = self.readline()
                 mode[imode, :, iatom].imag = list(map(float, line.split()))
 
-        info = dict(k=self._convert_k(k), parent=self._parent, gauge="r")
+        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):