added read_brillouinzone to multiple files

The hsx + TSHS files, plus siesta.stdout. Also fixed some typing issues. Signed-off-by: Nick Papior <[email protected]>
zerothi · Jan 22, 2025 · b1adf07 · b1adf07
1 parent 39aab90
commit b1adf07
Show file tree

Hide file tree

Showing 10 changed files with 208 additions and 23 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -9,6 +9,7 @@ we hit release version 1.0.0.
 ## [0.15.3] - YYYY-MM-DD
 
 ### Added
+- `read_brillouinzone` to multiple files
 - `hsxSileSiesta.write_hamiltonian`, allowed writing HSX files
   A long required functionality. This allows one to use the HSX file
   format in intrinsic behavior. TSHS file format will be deprecated in

diff --git a/src/sisl/io/siesta/_src/hsx_read.f90 b/src/sisl/io/siesta/_src/hsx_read.f90
@@ -294,6 +294,90 @@ subroutine read_hsx_ef1_2(fname, Ef)
 
 end subroutine read_hsx_ef1_2
 
+subroutine read_hsx_k(fname, kcell, kdispl)
+
+  implicit none
+
+  integer, parameter :: dp = selected_real_kind(p=15)
+
+  character(len=*), intent(in) :: fname
+  integer, intent(out) :: kcell(3,3)
+  real(dp), intent(out) :: kdispl(3)
+
+! Define f2py intents
+!f2py intent(in)  :: fname
+!f2py intent(out) :: kcell, kdispl
+
+  ! Internal variables and arrays
+  integer :: version
+
+  call read_hsx_version(fname, version)
+
+  if ( version < 2 ) then
+
+    kcell = 0
+    kcell(1,1) = 1
+    kcell(2,2) = 1
+    kcell(3,3) = 1
+
+  else
+
+    call read_hsx_k2(fname, kcell, kdispl)
+
+  end if
+
+end subroutine read_hsx_k
+
+subroutine read_hsx_k2(fname, kcell, kdispl)
+
+  implicit none
+
+  integer, parameter :: dp = selected_real_kind(p=15)
+
+  character(len=*), intent(in) :: fname
+  integer, intent(out) :: kcell(3,3)
+  real(dp), intent(out) :: kdispl(3)
+
+! Define f2py intents
+!f2py intent(in)  :: fname
+!f2py intent(out) :: kcell, kdispl
+
+  integer :: iu, version
+  integer :: na_u, no_u, nspin, nspecies
+
+  call open_file(fname, 'read', 'old', 'unformatted', iu)
+
+  read(iu, iostat=ierr) version
+  call iostat_update(ierr)
+  if ( version /= 2 ) then
+    call iostat_update(-3)
+    return
+  end if
+  read(iu, iostat=ierr) !is_dp
+  call iostat_update(ierr)
+  read(iu, iostat=ierr) na_u, no_u, nspin, nspecies!, nsc
+  call iostat_update(ierr)
+
+  read(iu, iostat=ierr) !ucell, Ef, qtot, temp
+  call iostat_update(ierr)
+  read(iu, iostat=ierr) !isc_off, xa, isa, lasto(1:na_u)
+  call iostat_update(ierr)
+
+  read(iu, iostat=ierr) !(label(is), zval(is), no(is), is=1,nspecies)
+  call iostat_update(ierr)
+  do is = 1, nspecies
+    read(iu, iostat=ierr) !(nquant(is,io), lquant(is,io), zeta(is,io), io=1,no(is))
+    call iostat_update(ierr)
+  end do
+  if ( version == 2 ) then
+     read(iu, iostat=ierr) kcell, kdispl
+     call iostat_update(ierr)
+  end if
+
+  call close_file(iu)
+
+end subroutine read_hsx_k2
+
 subroutine read_hsx_is_dp(fname, is_dp)
 
   implicit none

diff --git a/src/sisl/io/siesta/_src/precision.f90 b/src/sisl/io/siesta/_src/precision.f90
@@ -1,13 +1,10 @@
 module precision
 
   implicit none
-  
+
   integer, parameter :: r4 = selected_real_kind(6)
   integer, parameter :: r8 = selected_real_kind(15)
   integer, parameter :: i4 = selected_int_kind(6)
   integer, parameter :: i8 = selected_int_kind(18)
 
 end module precision
-
-
-
diff --git a/src/sisl/io/siesta/_src/tshs_read.f90 b/src/sisl/io/siesta/_src/tshs_read.f90
@@ -139,6 +139,41 @@ subroutine read_tshs_ef(fname, Ef)
 
 end subroutine read_tshs_ef
 
+subroutine read_tshs_k(fname, kcell, kdispl)
+
+  implicit none
+
+  integer, parameter :: dp = selected_real_kind(p=15)
+
+  character(len=*), intent(in) :: fname
+  integer, intent(out) :: kcell(3,3)
+  real(dp), intent(out) :: kdispl(3)
+
+! Define f2py intents
+!f2py intent(in)  :: fname
+!f2py intent(out) :: kcell, kdispl
+
+  integer :: iu, version
+
+  call open_file(fname, 'read', 'old', 'unformatted', iu)
+
+  read(iu, iostat=ierr) ! version
+  call iostat_update(ierr)
+  read(iu, iostat=ierr) ! na_u, no_u, no_s, nspin, n_nzsg
+  call iostat_update(ierr)
+  read(iu, iostat=ierr) ! nsc
+  call iostat_update(ierr)
+  read(iu, iostat=ierr) ! cell, xa
+  call iostat_update(ierr)
+  read(iu, iostat=ierr) ! Gamma, TSGamma, onlyS
+  call iostat_update(ierr)
+  read(iu, iostat=ierr) kcell, kdispl
+  call iostat_update(ierr)
+
+  call close_file(iu)
+
+end subroutine read_tshs_k
+
 subroutine read_tshs_cell(fname, n_s, nsc, cell, isc)
   use io_m, only: open_file, close_file
   use io_m, only: iostat_update
@@ -482,4 +517,3 @@ subroutine read_tshs_s(fname, no_u, nnz, ncol, list_col, S)
   call close_file(iu)
 
 end subroutine read_tshs_s
-
diff --git a/src/sisl/io/siesta/binaries.py b/src/sisl/io/siesta/binaries.py
@@ -342,6 +342,13 @@ def read_fermi_level(self) -> float:
         self._fortran_check("read_fermi_level", "could not read fermi-level.")
         return Ef
 
+    def read_brillouinzone(self, trs: bool = True) -> MonkhorstPack:
+        """Read the Brillouin zone object"""
+        kcell, kdispl = _siesta.read_tshs_k(self.file)
+        self._fortran_check("read_brillouinzone", "could not read the file.")
+        geom = self.read_geometry()
+        return MonkhorstPack(geom, kcell, displacement=kdispl, trs=trs)
+
 
 @set_module("sisl.io.siesta")
 class tshsSileSiesta(onlysSileSiesta):
@@ -1358,6 +1365,13 @@ def read_fermi_level(self, **kwargs) -> float:
 
         return Ef * _Ry2eV
 
+    def read_brillouinzone(self, trs: bool = True) -> MonkhorstPack:
+        """Read the Brillouin zone object"""
+        kcell, kdispl = _siesta.read_hsx_k(self.file)
+        self._fortran_check("read_brillouinzone", "could not read the file.")
+        geom = self.read_geometry()
+        return MonkhorstPack(geom, kcell, displacement=kdispl, trs=trs)
+
     def _r_hamiltonian_v0(self, **kwargs):
         geom = self.read_geometry(**kwargs)
 

diff --git a/src/sisl/io/siesta/kp.py b/src/sisl/io/siesta/kp.py
@@ -3,14 +3,14 @@
 # file, You can obtain one at https://mozilla.org/MPL/2.0/.
 from __future__ import annotations
 
-from typing import Optional, Union
+from typing import Optional
 
 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.typing import LatticeLike
 from sisl.unit.siesta import unit_convert
 
 from ..sile import add_sile, sile_fh_open, sile_raise_write
@@ -20,16 +20,14 @@
 
 Bohr2Ang = unit_convert("Bohr", "Ang")
 
-TLattice = Optional[Union[Lattice, LatticeChild]]
-
 
 @set_module("sisl.io.siesta")
 class kpSileSiesta(SileSiesta):
     """k-points file in 1/Bohr units"""
 
     @sile_fh_open()
     @deprecate_argument("sc", "lattice", "use lattice= instead of sc=", "0.15", "0.16")
-    def read_data(self, lattice: TLattice = None):
+    def read_data(self, lattice: Optional[LatticeLike] = None):
         """Returns K-points from the file (note that these are in reciprocal units)
 
         Parameters
@@ -68,7 +66,7 @@ def write_data(self, k, weight, fmt: str = ".9e"):
            k-points in units 1/Bohr
         weight : array_like
            same length as k, weights of k-points
-        fmt : str, optional
+        fmt :
            format for the k-values
         """
         sile_raise_write(self)
@@ -82,12 +80,12 @@ def write_data(self, k, weight, fmt: str = ".9e"):
 
     @sile_fh_open()
     @deprecate_argument("sc", "lattice", "use lattice= instead of sc=", "0.15", "0.16")
-    def read_brillouinzone(self, lattice: TLattice) -> BrillouinZone:
+    def read_brillouinzone(self, lattice: LatticeLike) -> BrillouinZone:
         """Returns K-points from the file (note that these are in reciprocal units)
 
         Parameters
         ----------
-        lattice : LatticeChild
+        lattice :
            required supercell for the BrillouinZone object
 
         Returns
@@ -97,20 +95,18 @@ def read_brillouinzone(self, lattice: TLattice) -> BrillouinZone:
         k, w = self.read_data(lattice)
         from sisl.physics.brillouinzone import BrillouinZone
 
-        bz = BrillouinZone(lattice)
-        bz._k = k
-        bz._w = w
+        bz = BrillouinZone(lattice, k=k, weight=w)
         return bz
 
     @sile_fh_open()
-    def write_brillouinzone(self, bz: BrillouinZone, fmt: str = ".9e"):
+    def write_brillouinzone(self, bz: BrillouinZone, fmt: str = ".9e") -> None:
         """Writes BrillouinZone-points to file
 
         Parameters
         ----------
-        bz : BrillouinZone
+        bz :
            object contain all weights and k-points
-        fmt : str, optional
+        fmt :
            format for the k-values
         """
         # And convert to 1/Bohr
@@ -148,12 +144,12 @@ def read_data(self):
 
     @sile_fh_open()
     @deprecate_argument("sc", "lattice", "use lattice= instead of sc=", "0.15", "0.16")
-    def read_brillouinzone(self, lattice: TLattice) -> BrillouinZone:
+    def read_brillouinzone(self, lattice: LatticeLike) -> BrillouinZone:
         """Returns K-points from the file
 
         Parameters
         ----------
-        lattice : LatticeChild
+        lattice :
            required supercell for the BrillouinZone object
 
         Returns
@@ -169,7 +165,7 @@ def read_brillouinzone(self, lattice: TLattice) -> BrillouinZone:
         return bz
 
     @sile_fh_open()
-    def write_brillouinzone(self, bz: BrillouinZone, fmt: str = ".9e"):
+    def write_brillouinzone(self, bz: BrillouinZone, fmt: str = ".9e") -> None:
         """Writes BrillouinZone-points to file
 
         Parameters

diff --git a/src/sisl/io/siesta/siesta_nc.py b/src/sisl/io/siesta/siesta_nc.py
@@ -19,6 +19,7 @@
     EnergyDensityMatrix,
     Hamiltonian,
 )
+from sisl.physics.brillouinzone import MonkhorstPack
 from sisl.physics.overlap import Overlap
 from sisl.unit.siesta import unit_convert
 
@@ -410,6 +411,14 @@ def read_grid(self, name, index=0, **kwargs) -> Grid:
 
         return grid
 
+    def read_brillouinzone(self, trs: bool = True) -> MonkhorstPack:
+        """Read the Brillouin zone object"""
+        settings = self.groups["BASIS"]
+        kcell = settings.variables["BZ"][:, :]
+        kdispl = settings.variables["BZ_displ"][:]
+        geom = self.read_geometry()
+        return MonkhorstPack(geom, kcell, displacement=kdispl, trs=trs)
+
     def write_basis(self, atoms: Atoms):
         """Write the current atoms orbitals as the basis
 

diff --git a/src/sisl/io/siesta/stdout.py b/src/sisl/io/siesta/stdout.py
@@ -17,6 +17,7 @@
 from sisl._internal import set_module
 from sisl.messages import deprecation, warn
 from sisl.physics import Spin
+from sisl.physics.brillouinzone import MonkhorstPack
 from sisl.unit.siesta import unit_convert
 from sisl.utils import PropertyDict
 from sisl.utils.cmd import *
@@ -769,6 +770,45 @@ def assign(out, key, val):
 
         return out
 
+    @SileBinder()
+    @sile_fh_open()
+    def read_brillouinzone(self, trs: bool = True) -> MonkhorstPack:
+        r"""Parses the k-grid section if present, otherwise returns the
+        :math:`\Gamma`-point"""
+
+        # store position, read geometry, then reposition file.
+        tell = self.fh.tell()
+        geom = self.read_geometry()
+        self.fh.seek(tell)
+
+        found, line = self.step_to(
+            "siesta: k-grid: Number of k-points", allow_reread=False
+        )
+
+        if not found:
+            return MonkhorstPack(geom, 1, trs=trs)
+
+        nk = int(line.split("=")[-1])
+
+        # default kcell and kdispl
+        kcell = np.diag([1, 1, 1])
+        kdispl = np.zeros(3)
+
+        # if next line also contains 'siesta: k-grid:' then we can step_to
+        line = self.readline()
+        if line.startswith("siesta: k-grid:"):
+            if self.step_to(
+                "siesta: k-grid: Supercell and displacements", allow_reread=False
+            )[0]:
+                for i in range(3):
+                    line = self.readline().split()
+                    kdispl[i] = float(line[-1])
+                    kcell[i, 2] = int(line[-2])
+                    kcell[i, 1] = int(line[-3])
+                    kcell[i, 0] = int(line[-4])
+
+        return MonkhorstPack(geom, kcell, displacement=kdispl, trs=trs)
+
     def read_data(self, *args, **kwargs) -> Any:
         """Read specific content in the Siesta out file
 

diff --git a/src/sisl/io/sile.py b/src/sisl/io/sile.py
@@ -1199,7 +1199,7 @@ def __iter__(self):
             yield l
             l = self.readline(comment=True)
 
-    def readline(self, comment=False):
+    def readline(self, comment: bool = False) -> str:
         r"""Reads the next line of the file"""
         l = self.fh.readline()
         self._line += 1