From b030984e17bc4c6f2342a68d867626329e7e1151 Mon Sep 17 00:00:00 2001
From: Nick Papior <nickpapior@gmail.com>
Date: Wed, 15 Jan 2025 15:17:57 +0100
Subject: [PATCH 1/2] enabled writing HSX files

This not only allows writing HSX files, but also
completes the version == 1, version == 2 specification.

The later version is going to be the new default
for transiesta runs, given that it contains all information
that the TSHS file contains + more.

Tests are added to ensure they are consistent.

However, the charges for the orbitals are not correctly
stored, that should be solved later on.

Signed-off-by: Nick Papior <nickpapior@gmail.com>
---
 CHANGELOG.md                          |   5 +
 src/sisl/_core/orbital.py             |   3 +
 src/sisl/io/siesta/CMakeLists.txt     |   4 -
 src/sisl/io/siesta/_src/hs_read.f90   | 144 -----------------
 src/sisl/io/siesta/_src/hsx_read.f90  | 177 +++++++++++++++------
 src/sisl/io/siesta/_src/hsx_write.f90 | 201 +++++++++++++++---------
 src/sisl/io/siesta/binaries.py        | 216 +++++++++++++++++++++++---
 src/sisl/io/siesta/tests/test_hsx.py  |  39 +++++
 src/sisl/io/sile.py                   |   7 +-
 9 files changed, 496 insertions(+), 300 deletions(-)
 delete mode 100644 src/sisl/io/siesta/_src/hs_read.f90

diff --git a/CHANGELOG.md b/CHANGELOG.md
index f1985eb43d..62c0305114 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,11 @@ we hit release version 1.0.0.
 ## [0.15.3] - YYYY-MM-DD
 
 ### Added
+- `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
+  later Siesta releases, see https://gitlab.com/siesta-project/siesta/-/issues/183
+  This allows one to write HSX files in version 1, or 2, 2 is default.
 - `astype` for sparse matrices, enables one to change data-types, #865
   This should be preferred over transform which can't do real->complex
   of spin matrices.
diff --git a/src/sisl/_core/orbital.py b/src/sisl/_core/orbital.py
index 7f95889db2..a73945c733 100644
--- a/src/sisl/_core/orbital.py
+++ b/src/sisl/_core/orbital.py
@@ -1041,6 +1041,9 @@ def __init__(self, *args, **kwargs):
         if n <= 0:
             raise ValueError(f"{self.__class__.__name__} n must be >= 1")
 
+        if zeta <= 0:
+            raise ValueError(f"{self.__class__.__name__} zeta must be >= 1")
+
         if self.l >= len(_rspher_harm_fact):
             raise ValueError(
                 f"{self.__class__.__name__} does not implement shells l>={len(_rspher_harm_fact)}!"
diff --git a/src/sisl/io/siesta/CMakeLists.txt b/src/sisl/io/siesta/CMakeLists.txt
index fc0bdd34dd..5b4fccfd54 100644
--- a/src/sisl/io/siesta/CMakeLists.txt
+++ b/src/sisl/io/siesta/CMakeLists.txt
@@ -14,7 +14,6 @@ set(sig_sources
     grid_read grid_write
     gf_read gf_write
     tsde_read tsde_write
-    hs_read
     wfsx_read
   )
 # The REGEX allows Option 2 to pass through correctly
@@ -30,9 +29,6 @@ set_source_files_properties("${name}"
   )
 
 get_filename_component(name _src/hsx_write.f90 NAME)
-set_source_files_properties("${name}"
-  PROPERTIES SISL_SIGNATURE_SKIP "write_hsx1"
-  )
 
 # Additional sources which we do not need interfaces
 # to in the Python front
diff --git a/src/sisl/io/siesta/_src/hs_read.f90 b/src/sisl/io/siesta/_src/hs_read.f90
deleted file mode 100644
index 9ecf2a8524..0000000000
--- a/src/sisl/io/siesta/_src/hs_read.f90
+++ /dev/null
@@ -1,144 +0,0 @@
-! 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/.
-subroutine read_hs_header(fname, Gamma, nspin, no_u, no_s, maxnh)
-  use io_m, only: open_file, close_file
-  use io_m, only: iostat_update
-
-  implicit none
-
-  ! Input parameters
-  character(len=*), intent(in) :: fname
-  logical, intent(out) :: Gamma
-  integer, intent(out) :: no_u, no_s, nspin, maxnh
-
-! Define f2py intents
-!f2py intent(in)  :: fname
-!f2py intent(out) :: Gamma, no_u, no_s, nspin, maxnh
-
-! Internal variables and arrays
-  integer :: iu, ierr
-
-  call open_file(fname, 'read', 'old', 'unformatted', iu)
-
-  read(iu, iostat=ierr) no_u, no_s, nspin, maxnh
-  call iostat_update(ierr)
-
-  read(iu, iostat=ierr) Gamma
-  call iostat_update(ierr)
-
-  call close_file(iu)
-
-end subroutine read_hs_header
-
-subroutine read_hs(fname, nspin, no_u,no_s,maxnh, &
-    numh,listhptr,listh,H,S,xij)
-  use io_m, only: open_file, close_file
-  use io_m, only: iostat_update
-
-  implicit none
-
-  ! Precision 
-  integer, parameter :: sp = selected_real_kind(p=6)
-  integer, parameter :: dp = selected_real_kind(p=15)
-
-  ! Input parameters
-  character(len=*), intent(in) :: fname
-  integer, intent(in) :: no_u, no_s, nspin, maxnh
-  integer, intent(out) :: listh(maxnh), numh(no_u), listhptr(no_u)
-  real(dp), intent(out) :: H(maxnh,nspin), S(maxnh), xij(3,maxnh)
-
-! Define f2py intents
-!f2py intent(in)  :: fname
-!f2py intent(in) :: no_u, no_s, nspin, maxnh
-!f2py intent(out) :: numh, listhptr, listh
-!f2py intent(out) :: H, S, xij
-
-! Internal variables and arrays
-  integer :: iu, ierr
-  integer :: is, ih, im, ip
-
-  ! Local readables
-  logical :: Gamma
-  integer :: lno_s, lno_u, lnspin, lmaxnh
-
-  call open_file(fname, 'read', 'old', 'unformatted', iu)
-
-  read(iu, iostat=ierr) lno_u, lno_s, lnspin, lmaxnh
-  call iostat_update(ierr)
-  if ( lno_u /= no_u ) stop 'Error in reading data, not allocated, no_u'
-  if ( lno_s /= no_s ) stop 'Error in reading data, not allocated, no_s'
-  if ( lnspin /= nspin ) stop 'Error in reading data, not allocated, nspin'
-  if ( lmaxnh /= maxnh ) stop 'Error in reading data, not allocated, maxnh'
-
-  read(iu, iostat=ierr) Gamma
-  call iostat_update(ierr)
-  if ( Gamma ) then
-    if ( no_u /= no_s ) then
-      stop 'Error in reading data, not allocated, Gamma'
-    end if
-  else if ( no_u == no_s ) then
-    stop 'Error in reading data, not allocated, Gamma'
-  end if
-
-! Read out indxuo
-  if (.not. Gamma) then
-    read(iu, iostat=ierr) ! indxuo
-    call iostat_update(ierr)
-  end if
-
-  read(iu, iostat=ierr) numh
-  call iostat_update(ierr)
-
-! Create listhptr
-  listhptr(1) = 0
-  do ih = 2 , no_u
-    listhptr(ih) = listhptr(ih-1) + numh(ih-1)
-  end do
-
-! Read listh
-  do ih = 1 , no_u
-    ip = listhptr(ih)
-    do im = 1 , numh(ih)
-      read(iu, iostat=ierr) listh(ip+im)
-      call iostat_update(ierr)
-    end do
-  end do
-
-! Read Hamiltonian
-  do is = 1 , nspin
-    do ih = 1 , no_u
-      ip = listhptr(ih)
-      do im = 1 , numh(ih)
-        read(iu, iostat=ierr) H(ip+im,is)
-        call iostat_update(ierr)
-      end do
-    end do
-  end do
-
-! Read Hamiltonian
-  do ih = 1 , no_u
-    ip = listhptr(ih)
-    do im = 1 , numh(ih)
-      read(iu, iostat=ierr) S(ip+im)
-      call iostat_update(ierr)
-    end do
-  end do
-
-  read(iu, iostat=ierr) !qtot, temp 
-  call iostat_update(ierr)
-
-! Read xij
-  if ( .not. Gamma ) then
-    do ih = 1 , no_u
-      ip = listhptr(ih)
-      do im = 1 , numh(ih)
-        read(iu, iostat=ierr) xij(1:3,ip + im)
-        call iostat_update(ierr)
-      end do
-    end do
-  end if
-
-  call close_file(iu)
-
-end subroutine read_hs
diff --git a/src/sisl/io/siesta/_src/hsx_read.f90 b/src/sisl/io/siesta/_src/hsx_read.f90
index 33fbadbb46..17baf94f20 100644
--- a/src/sisl/io/siesta/_src/hsx_read.f90
+++ b/src/sisl/io/siesta/_src/hsx_read.f90
@@ -34,7 +34,7 @@ subroutine read_hsx_version(fname, version)
     rewind(iu)
     read(iu, iostat=ierr) version
 
-    if ( version /= 1 ) then
+    if ( all(version /= [1, 2]) ) then
       ! Signal we do not know about this file
       call iostat_update(-3)
     end if
@@ -68,9 +68,9 @@ subroutine read_hsx_sizes(fname, nspin, na_u, no_u, no_s, maxnh)
 
     call read_hsx_sizes0(fname, nspin, na_u, no_u, no_s, maxnh)
 
-  else if ( version == 1 ) then
+  else if ( any(version == [1, 2]) ) then
 
-    call read_hsx_sizes1(fname, nspin, na_u, no_u, no_s, maxnh)
+    call read_hsx_sizes1_2(fname, nspin, na_u, no_u, no_s, maxnh)
 
   end if
 
@@ -132,7 +132,7 @@ subroutine read_hsx_sizes0(fname, nspin, na_u, no_u, no_s, maxnh)
 
   do is = 1, nspecies
     do io = 1 , no(is)
-      read(iu, iostat=ierr) !n(is,io), l(is,io), zeta(is,io)
+      read(iu, iostat=ierr) !n(io,is), l(io,is), zeta(io,is)
       call iostat_update(ierr)
     end do
   end do
@@ -145,7 +145,7 @@ subroutine read_hsx_sizes0(fname, nspin, na_u, no_u, no_s, maxnh)
 
 end subroutine read_hsx_sizes0
 
-subroutine read_hsx_sizes1(fname, nspin, na_u, no_u, no_s, maxnh)
+subroutine read_hsx_sizes1_2(fname, nspin, na_u, no_u, no_s, maxnh)
   use io_m, only: open_file, close_file
   use io_m, only: iostat_update
 
@@ -160,7 +160,6 @@ subroutine read_hsx_sizes1(fname, nspin, na_u, no_u, no_s, maxnh)
 !f2py intent(out) :: nspin, na_u, no_u, no_s, maxnh
 
   ! Internal variables and arrays
-  logical :: Gamma
   integer :: version, nspecies, nsc(3)
   integer, allocatable :: numh(:)
   integer :: is
@@ -170,7 +169,7 @@ subroutine read_hsx_sizes1(fname, nspin, na_u, no_u, no_s, maxnh)
 
   read(iu, iostat=ierr) version
   call iostat_update(ierr)
-  if ( version /= 1 ) then
+  if ( all(version /= [1, 2]) ) then
     call iostat_update(-3)
     return
   end if
@@ -179,19 +178,22 @@ subroutine read_hsx_sizes1(fname, nspin, na_u, no_u, no_s, maxnh)
   read(iu, iostat=ierr) na_u, no_u, nspin, nspecies, nsc
   call iostat_update(ierr)
   no_s = product(nsc) * no_u
-  Gamma = no_s == no_u
 
   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
 
   allocate(numh(no_u))
   read(iu, iostat=ierr) numh
@@ -201,7 +203,7 @@ subroutine read_hsx_sizes1(fname, nspin, na_u, no_u, no_s, maxnh)
 
   call close_file(iu)
 
-end subroutine read_hsx_sizes1
+end subroutine read_hsx_sizes1_2
 
 subroutine read_hsx_ef(fname, Ef)
 
@@ -222,9 +224,9 @@ subroutine read_hsx_ef(fname, Ef)
 
   call read_hsx_version(fname, version)
 
-  if ( version == 1 ) then
+  if ( any(version == [1, 2]) ) then
 
-    call read_hsx_ef1(fname, Ef)
+    call read_hsx_ef1_2(fname, Ef)
 
   end if
 
@@ -251,7 +253,7 @@ subroutine read_hsx_ef0(fname, Ef)
 
 end subroutine read_hsx_ef0
 
-subroutine read_hsx_ef1(fname, Ef)
+subroutine read_hsx_ef1_2(fname, Ef)
   use io_m, only: open_file, close_file
   use io_m, only: iostat_update
 
@@ -276,7 +278,7 @@ subroutine read_hsx_ef1(fname, Ef)
 
   read(iu, iostat=ierr) version
   call iostat_update(ierr)
-  if ( version /= 1 ) then
+  if ( all(version /= [1, 2]) ) then
     call iostat_update(-3)
     return
   end if
@@ -290,14 +292,76 @@ subroutine read_hsx_ef1(fname, Ef)
 
   call close_file(iu)
 
-end subroutine read_hsx_ef1
+end subroutine read_hsx_ef1_2
+
+subroutine read_hsx_is_dp(fname, is_dp)
+
+  implicit none
+
+  ! Input parameters
+  character(len=*), intent(in) :: fname
+  logical, intent(out) :: is_dp
+
+! Define f2py intents
+!f2py intent(in)  :: fname
+!f2py intent(out) :: is_dp
+
+  ! Internal variables and arrays
+  integer :: version
+
+  call read_hsx_version(fname, version)
+
+  if ( version == 0 ) then
+     is_dp = .false.
+
+  else if ( any(version == [1, 2]) ) then
+
+    call read_hsx_is_dp1_2(fname, is_dp)
+
+  end if
+
+end subroutine read_hsx_is_dp
+
+
+subroutine read_hsx_is_dp1_2(fname, is_dp)
+  use io_m, only: open_file, close_file
+  use io_m, only: iostat_update
+
+  implicit none
+
+  ! Input parameters
+  character(len=*), intent(in) :: fname
+  logical, intent(out) :: is_dp
+
+! Define f2py intents
+!f2py intent(in)  :: fname
+!f2py intent(out) :: is_dp
+
+  ! Internal variables and arrays
+  integer :: version
+  integer :: iu, ierr
+
+  call open_file(fname, 'read', 'old', 'unformatted', iu)
+
+  read(iu, iostat=ierr) version
+  call iostat_update(ierr)
+  if ( all(version /= [1, 2]) ) then
+    call iostat_update(-3)
+    return
+  end if
+  read(iu, iostat=ierr) is_dp
+  call iostat_update(ierr)
+
+  call close_file(iu)
+
+end subroutine read_hsx_is_dp1_2
 
 
 !< Internal method for skipping species information in version 1 and later
 !< The unit *must* be located just after
 !<     read(iu) isc_off, xa, isa, lasto(1:na_u)
 
-subroutine internal_read_hsx_skip_species1(iu, nspecies)
+subroutine internal_read_hsx_skip_species1_2(iu, nspecies)
   use io_m, only: iostat_update
 
   implicit none
@@ -316,7 +380,7 @@ subroutine internal_read_hsx_skip_species1(iu, nspecies)
     call iostat_update(ierr)
   end do
 
-end subroutine internal_read_hsx_skip_species1
+end subroutine internal_read_hsx_skip_species1_2
 
 
 subroutine read_hsx_hsx0(fname, nspin, no_u, no_s, maxnh, &
@@ -427,7 +491,7 @@ subroutine read_hsx_hsx0(fname, nspin, no_u, no_s, maxnh, &
 
 end subroutine read_hsx_hsx0
 
-subroutine read_hsx_hsx1(fname, nspin, no_u, no_s, maxnh, &
+subroutine read_hsx_hsx1_2(fname, nspin, no_u, no_s, maxnh, &
     numh, listh, H, S, isc)
   use io_m, only: open_file, close_file
   use io_m, only: iostat_update
@@ -467,7 +531,7 @@ subroutine read_hsx_hsx1(fname, nspin, no_u, no_s, maxnh, &
   ! Read overall data
   read(iu, iostat=ierr) version
   call iostat_update(ierr)
-  if ( version /= 1 ) then
+  if ( all(version /= [1, 2]) ) then
     call iostat_update(-3)
     return
   end if
@@ -487,7 +551,12 @@ subroutine read_hsx_hsx1(fname, nspin, no_u, no_s, maxnh, &
   read(iu, iostat=ierr) isc!, xa, isa, lasto
   call iostat_update(ierr)
 
-  call internal_read_hsx_skip_species1(iu, nspecies)
+  call internal_read_hsx_skip_species1_2(iu, nspecies)
+
+  if ( version == 2 ) then
+    read(iu, iostat=ierr) ! kcell, kdispl
+    call iostat_update(ierr)
+  end if
 
   read(iu, iostat=ierr) numh
   call iostat_update(ierr)
@@ -554,7 +623,7 @@ subroutine read_hsx_hsx1(fname, nspin, no_u, no_s, maxnh, &
 
   call close_file(iu)
 
-end subroutine read_hsx_hsx1
+end subroutine read_hsx_hsx1_2
 
 
 subroutine read_hsx_sx0(fname, nspin, no_u, no_s, maxnh, &
@@ -669,7 +738,7 @@ subroutine read_hsx_sx0(fname, nspin, no_u, no_s, maxnh, &
 
 end subroutine read_hsx_sx0
 
-subroutine read_hsx_sx1(fname, nspin, no_u, no_s, maxnh, &
+subroutine read_hsx_sx1_2(fname, nspin, no_u, no_s, maxnh, &
     numh, listh, S, isc)
   use io_m, only: open_file, close_file
   use io_m, only: iostat_update
@@ -709,7 +778,7 @@ subroutine read_hsx_sx1(fname, nspin, no_u, no_s, maxnh, &
   ! Read overall data
   read(iu, iostat=ierr) version
   call iostat_update(ierr)
-  if ( version /= 1 ) then
+  if ( all(version /= [1, 2]) ) then
     call iostat_update(-3)
     return
   end if
@@ -729,7 +798,12 @@ subroutine read_hsx_sx1(fname, nspin, no_u, no_s, maxnh, &
   read(iu, iostat=ierr) isc!, xa, isa, lasto
   call iostat_update(ierr)
 
-  call internal_read_hsx_skip_species1(iu, nspecies)
+  call internal_read_hsx_skip_species1_2(iu, nspecies)
+
+  if ( version == 2 ) then
+    read(iu, iostat=ierr) ! kcell, kdispl
+    call iostat_update(ierr)
+  end if
 
   read(iu, iostat=ierr) numh
   call iostat_update(ierr)
@@ -783,9 +857,9 @@ subroutine read_hsx_sx1(fname, nspin, no_u, no_s, maxnh, &
 
   call close_file(iu)
 
-end subroutine read_hsx_sx1
+end subroutine read_hsx_sx1_2
 
-subroutine read_hsx_geom1(fname, na_u, cell, nsc, xa, lasto)
+subroutine read_hsx_geom1_2(fname, na_u, cell, nsc, xa, lasto)
   use io_m, only: open_file, close_file
   use io_m, only: iostat_update
 
@@ -812,7 +886,7 @@ subroutine read_hsx_geom1(fname, na_u, cell, nsc, xa, lasto)
   ! Read overall data
   read(iu, iostat=ierr) version
   call iostat_update(ierr)
-  if ( version /= 1 ) then
+  if ( all(version /= [1, 2]) ) then
     call iostat_update(-3)
     return
   end if
@@ -834,7 +908,7 @@ subroutine read_hsx_geom1(fname, na_u, cell, nsc, xa, lasto)
 
   call close_file(iu)
 
-end subroutine read_hsx_geom1
+end subroutine read_hsx_geom1_2
 
 subroutine read_hsx_species_sizes(fname, no_u, na_u, nspecies)
 
@@ -857,8 +931,8 @@ subroutine read_hsx_species_sizes(fname, no_u, na_u, nspecies)
 
   if ( version == 0 ) then
     call read_hsx_species_sizes0(fname, no_u, na_u, nspecies)
-  else if ( version == 1 ) then
-    call read_hsx_species_sizes1(fname, no_u, na_u, nspecies)
+  else if ( any(version == [1, 2]) ) then
+    call read_hsx_species_sizes1_2(fname, no_u, na_u, nspecies)
   end if
 
 end subroutine read_hsx_species_sizes
@@ -915,7 +989,7 @@ subroutine read_hsx_species_sizes0(fname, no_u, na_u, nspecies)
     call iostat_update(ierr)
   end do
 
-  ! Read Hamiltonian
+  ! Read Hamiltonian + overlap
   do is = 1, no_u * (lnspin + 1)
     read(iu, iostat=ierr) ! H and S
     call iostat_update(ierr)
@@ -950,7 +1024,7 @@ subroutine read_hsx_species_sizes0(fname, no_u, na_u, nspecies)
 
 end subroutine read_hsx_species_sizes0
 
-subroutine read_hsx_species_sizes1(fname, no_u, na_u, nspecies)
+subroutine read_hsx_species_sizes1_2(fname, no_u, na_u, nspecies)
   use io_m, only: open_file, close_file
   use io_m, only: iostat_update
 
@@ -977,7 +1051,7 @@ subroutine read_hsx_species_sizes1(fname, no_u, na_u, nspecies)
   ! Read overall data
   read(iu, iostat=ierr) version
   call iostat_update(ierr)
-  if ( version /= 1 ) then
+  if ( all(version /= [1, 2]) ) then
     call iostat_update(-3)
     return
   end if
@@ -991,7 +1065,7 @@ subroutine read_hsx_species_sizes1(fname, no_u, na_u, nspecies)
 
   call close_file(iu)
 
-end subroutine read_hsx_species_sizes1
+end subroutine read_hsx_species_sizes1_2
 
 subroutine read_hsx_species_info(fname, nspecies, no_u, na_u, label, zval, no, isa)
 
@@ -1020,9 +1094,9 @@ subroutine read_hsx_species_info(fname, nspecies, no_u, na_u, label, zval, no, i
 
     call read_hsx_species_info0(fname, nspecies, no_u, na_u, label, zval, no, isa)
 
-  else if ( version == 1 ) then
+  else if ( any(version == [1, 2]) ) then
 
-    call read_hsx_species_info1(fname, nspecies, no_u, na_u, label, zval, no, isa)
+    call read_hsx_species_info1_2(fname, nspecies, no_u, na_u, label, zval, no, isa)
 
   end if
 
@@ -1040,7 +1114,7 @@ subroutine read_hsx_species_info0(fname, nspecies, no_u, na_u, label, zval, no,
   ! Input parameters
   character(len=*), intent(in) :: fname
   integer, intent(in) :: nspecies, no_u, na_u
-  character(len=1), intent(out) :: label(20,nspecies)
+  character(len=20), intent(out) :: label(nspecies)
   real(dp), intent(out) :: zval(nspecies)
   integer, intent(out) :: no(nspecies)
   integer, intent(out) :: isa(na_u)
@@ -1112,7 +1186,7 @@ subroutine read_hsx_species_info0(fname, nspecies, no_u, na_u, label, zval, no,
   if ( lnspecies /= nspecies ) stop 'Error in reading data, not allocated, nspecies'
   call iostat_update(ierr)
 
-  read(iu, iostat=ierr) (label(1:20,is), zval(is),no(is),is=1,nspecies)
+  read(iu, iostat=ierr) (label(is), zval(is),no(is),is=1,nspecies)
   call iostat_update(ierr)
   do is = 1, nspecies
     do io = 1 , no(is)
@@ -1131,7 +1205,7 @@ subroutine read_hsx_species_info0(fname, nspecies, no_u, na_u, label, zval, no,
 
 end subroutine read_hsx_species_info0
 
-subroutine read_hsx_species_info1(fname, nspecies, no_u, na_u, label, zval, no, isa)
+subroutine read_hsx_species_info1_2(fname, nspecies, no_u, na_u, label, zval, no, isa)
   use io_m, only: open_file, close_file
   use io_m, only: iostat_update
 
@@ -1142,7 +1216,7 @@ subroutine read_hsx_species_info1(fname, nspecies, no_u, na_u, label, zval, no,
   ! Input parameters
   character(len=*), intent(in) :: fname
   integer, intent(in) :: nspecies, no_u, na_u
-  character(len=1), intent(out) :: label(20,nspecies)
+  character(len=20), intent(out) :: label(nspecies)
   real(dp), intent(out) :: zval(nspecies)
   integer, intent(out) :: no(nspecies)
   integer, intent(out) :: isa(na_u)
@@ -1167,7 +1241,7 @@ subroutine read_hsx_species_info1(fname, nspecies, no_u, na_u, label, zval, no,
   ! Read overall data
   read(iu, iostat=ierr) version
   call iostat_update(ierr)
-  if ( version /= 1 ) then
+  if ( all(version /= [1, 2]) ) then
     call iostat_update(-3)
     return
   end if
@@ -1177,6 +1251,7 @@ subroutine read_hsx_species_info1(fname, nspecies, no_u, na_u, label, zval, no,
 
   read(iu, iostat=ierr) lna_u, lno_u, lnspin, lnspecies, nsc
   call iostat_update(ierr)
+  if ( lna_u /= na_u ) call iostat_update(-6)
   if ( lnspecies /= nspecies ) call iostat_update(-6)
   if ( lno_u /= no_u ) call iostat_update(-6)
 
@@ -1189,12 +1264,12 @@ subroutine read_hsx_species_info1(fname, nspecies, no_u, na_u, label, zval, no,
   call iostat_update(ierr)
   deallocate(isc, xa)
 
-  read(iu, iostat=ierr) (label(1:20,is), zval(is),no(is),is=1,nspecies)
+  read(iu, iostat=ierr) (label(is), zval(is), no(is),is=1,nspecies)
   call iostat_update(ierr)
 
   call close_file(iu)
 
-end subroutine read_hsx_species_info1
+end subroutine read_hsx_species_info1_2
 
 subroutine read_hsx_species(fname, ispecies, no_species, n_species, l_species, zeta_species)
 
@@ -1219,9 +1294,9 @@ subroutine read_hsx_species(fname, ispecies, no_species, n_species, l_species, z
 
     call read_hsx_species0(fname, ispecies, no_species, n_species, l_species, zeta_species)
 
-  else if ( version == 1 ) then
+  else if ( any(version == [1, 2]) ) then
 
-    call read_hsx_species1(fname, ispecies, no_species, n_species, l_species, zeta_species)
+    call read_hsx_species1_2(fname, ispecies, no_species, n_species, l_species, zeta_species)
 
   end if
 
@@ -1313,7 +1388,7 @@ subroutine read_hsx_species0(fname, ispecies, no_species, n_species, l_species,
   allocate(zval(lnspecies))
   allocate(no(lnspecies))
 
-  read(iu, iostat=ierr) (label(is), zval(is),no(is),is=1,lnspecies)
+  read(iu, iostat=ierr) (label(is), zval(is), no(is),is=1,lnspecies)
   call iostat_update(ierr)
   do is = 1, lnspecies
     if ( is == ispecies ) then
@@ -1323,7 +1398,7 @@ subroutine read_hsx_species0(fname, ispecies, no_species, n_species, l_species,
       end do
     else
       do io = 1 , no(is)
-        read(iu, iostat=ierr) !n(is,io), l(is,io), zeta(is,io)
+        read(iu, iostat=ierr) !n(io,is), l(io,is), zeta(io,is)
         call iostat_update(ierr)
       end do
     end if
@@ -1335,7 +1410,7 @@ subroutine read_hsx_species0(fname, ispecies, no_species, n_species, l_species,
 
 end subroutine read_hsx_species0
 
-subroutine read_hsx_species1(fname, ispecies, no_species, n_species, l_species, zeta_species)
+subroutine read_hsx_species1_2(fname, ispecies, no_species, n_species, l_species, zeta_species)
   use io_m, only: open_file, close_file
   use io_m, only: iostat_update
 
@@ -1365,7 +1440,7 @@ subroutine read_hsx_species1(fname, ispecies, no_species, n_species, l_species,
   ! Read overall data
   read(iu, iostat=ierr) version
   call iostat_update(ierr)
-  if ( version /= 1 ) then
+  if ( all(version /= [1, 2]) ) then
     call iostat_update(-3)
     return
   end if
@@ -1394,11 +1469,11 @@ subroutine read_hsx_species1(fname, ispecies, no_species, n_species, l_species,
       read(iu, iostat=ierr) (n_species(io), l_species(io), zeta_species(io), io=1,no_species)
       call iostat_update(ierr)
     else
-      read(iu, iostat=ierr) !n(is,io), l(is,io), zeta(is,io)
+      read(iu, iostat=ierr) !n(io,is), l(io,is), zeta(io,is)
       call iostat_update(ierr)
     end if
   end do
 
   call close_file(iu)
 
-end subroutine read_hsx_species1
+end subroutine read_hsx_species1_2
diff --git a/src/sisl/io/siesta/_src/hsx_write.f90 b/src/sisl/io/siesta/_src/hsx_write.f90
index 59b5beefc2..060a5ad242 100644
--- a/src/sisl/io/siesta/_src/hsx_write.f90
+++ b/src/sisl/io/siesta/_src/hsx_write.f90
@@ -1,8 +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/.
-subroutine write_hsx0(fname, Gamma, no_u, no_s, nspin, maxnh, &
-    numh, listhptr, listh, H, S, xij, Qtot, temp)
+subroutine write_hsx1(fname, &
+    label, Z, no, n, l, zeta, &
+    xa, isa, lasto, &
+    ucell, nsc, isc_off, &
+    ncol, row, &
+    H, S, is_dp, &
+    Ef, Qtot, temp, &
+    nspecies, na_u, no_max, n_s, nspin, no_u, nnz)
 
   use precision, only: r4, r8
   use io_m, only: open_file, close_file
@@ -11,88 +17,124 @@ subroutine write_hsx0(fname, Gamma, no_u, no_s, nspin, maxnh, &
   implicit none
 
   ! Input parameters
-  character(len=*) :: fname
-  logical :: Gamma
-  integer :: no_u, no_s, nspin, maxnh
-  integer :: listh(maxnh), numh(no_u), listhptr(no_u)
-  real(r8) :: H(maxnh,nspin), S(maxnh), xij(3,maxnh), Qtot, temp
+  character(len=*), intent(in) :: fname
+  integer, intent(in) :: nspecies
+  integer, intent(in) :: na_u, isa(na_u), lasto(na_u)
+  real(r8), intent(in) :: xa(3,na_u)
+  real(r8), intent(in) :: ucell(3,3)
+  character(len=20), intent(in) :: label(nspecies)
+  real(r8), intent(in) :: Z(nspecies)
+  integer, intent(in) :: no(nspecies)
+  integer, intent(in) :: no_max
+  integer, intent(in), dimension(no_max, nspecies) :: n, l, zeta
+  integer, intent(in) :: n_s, nsc(3), isc_off(3, n_s)
+  integer, intent(in) :: nspin, no_u, nnz
+  integer, intent(in) :: ncol(no_u), row(nnz)
+  real(r8), intent(in) :: H(nnz, nspin), S(nnz)
+  logical, intent(in) :: is_dp
+  real(r8), intent(in) :: Ef, Qtot, temp
 
 ! Define f2py intents
-!f2py intent(in) :: fname, Gamma, no_u, no_s, nspin, maxnh
-!f2py intent(in) :: numh, listhptr, listh, H, S, xij, Qtot, temp
+!f2py intent(in) :: fname
+!f2py intent(in) :: na_u, xa, isa, lasto
+!f2py intent(in) :: nspecies, label, Z, no
+!f2py intent(in) :: no_max, n, l, zeta
+!f2py intent(in) :: ucell, n_s, nsc, isc_off
+!f2py intent(in) :: no_u, nnz, nspin
+!f2py intent(in) :: ncol, row
+!f2py intent(in) :: H, S, is_dp, Ef, Qtot, temp
 
 ! Internal variables and arrays
   integer :: iu, ierr
-  integer :: is, ih, im, k
-  integer :: indxuo(no_s)
+  integer :: i, is, idx
 
   ! Open file (ensure we start from a clean slate)!
   call open_file(fname, 'write', 'unknown', 'unformatted', iu)
 
-! Write overall data
-  write(iu, iostat=ierr) no_u, no_s, nspin, maxnh
+  write(iu, iostat=ierr) 1
   call iostat_update(ierr)
 
-! Write logical
-  write(iu, iostat=ierr) Gamma
+! Write header information
+  write(iu, iostat=ierr) is_dp
   call iostat_update(ierr)
 
-! Write out indxuo
-  if (.not. Gamma) then
-    do ih = 1 , no_s
-      im = mod(ih,no_u)
-      if ( im == 0 ) im = no_u
-      indxuo(ih) = im
-    end do
-    write(iu, iostat=ierr) (indxuo(ih),ih=1,no_s)
+! Write overall data
+  write(iu, iostat=ierr) na_u, no_u, nspin, nspecies, nsc
+  call iostat_update(ierr)
+  write(iu, iostat=ierr) ucell, Ef, qtot, temp
+  call iostat_update(ierr)
+  write(iu, iostat=ierr) isc_off, xa, isa, lasto
+  call iostat_update(ierr)
+
+  write(iu, iostat=ierr) (label(is), Z(is), no(is), is=1,nspecies)
+  call iostat_update(ierr)
+  do is = 1, nspecies
+    write(iu, iostat=ierr) (n(i,is), l(i,is), zeta(i,is), i=1,no(is))
     call iostat_update(ierr)
-  end if
+  end do
 
-  write(iu, iostat=ierr) (numh(ih),ih=1,no_u)
+  write(iu, iostat=ierr) ncol
   call iostat_update(ierr)
 
-! Write listh
-  do ih = 1 , no_u
-    write(iu, iostat=ierr) (listh(listhptr(ih)+im),im = 1,numh(ih))
+  idx = 1
+  do i = 1, no_u
+    write(iu, iostat=ierr) row(idx:idx+ncol(i)-1)
     call iostat_update(ierr)
+    idx = idx + ncol(i)
   end do
 
-! Write Hamiltonian
-  do is = 1 , nspin
-    do ih = 1 , no_u
-      write(iu, iostat=ierr) (real(H(listhptr(ih)+im,is),kind=r4),im=1,numh(ih))
+  if ( is_dp ) then
+
+    do is = 1 , nspin
+      idx = 1
+      do i = 1 , no_u
+        write(iu, iostat=ierr) H(idx:idx+ncol(i)-1,is)
+        call iostat_update(ierr)
+        idx = idx + ncol(i)
+      end do
+    end do
+
+    idx = 1
+    do i = 1 , no_u
+      write(iu, iostat=ierr) S(idx:idx+ncol(i)-1)
       call iostat_update(ierr)
+      idx = idx + ncol(i)
     end do
-  end do
 
-! Write Overlap matrix
-  do ih = 1,no_u
-    write(iu, iostat=ierr) (real(S(listhptr(ih)+im),kind=r4),im = 1,numh(ih))
-    call iostat_update(ierr)
-  end do
+  else
 
-  write(iu, iostat=ierr) Qtot,temp
-  call iostat_update(ierr)
+    do is = 1 , nspin
+      idx = 1
+      do i = 1 , no_u
+        write(iu, iostat=ierr) real(H(idx:idx+ncol(i)-1,is), r4)
+        call iostat_update(ierr)
+        idx = idx + ncol(i)
+      end do
+    end do
 
-  do ih = 1 , no_u
-    write(iu, iostat=ierr) ((real(xij(k,listhptr(ih)+im),kind=r4), k=1,3),im =1,numh(ih))
-    call iostat_update(ierr)
-  end do
+    idx = 1
+    do i = 1 , no_u
+      write(iu, iostat=ierr) real(S(idx:idx+ncol(i)-1), r4)
+      call iostat_update(ierr)
+      idx = idx + ncol(i)
+    end do
+
+  end if
 
   call close_file(iu)
 
-end subroutine write_hsx0
+end subroutine write_hsx1
 
-subroutine write_hsx1(fname, &
-    nspecies, &
-    na_u, xa, isa, lasto, &
-    label, Z, no_max, no, n, l, zeta, &
-    ucell, n_s, nsc, isc_off, &
-    nspin, &
-    no_u, no_s, nnz, &
+subroutine write_hsx2(fname, &
+    label, Z, no, n, l, zeta, &
+    xa, isa, lasto, &
+    ucell, nsc, isc_off, &
     ncol, row, &
     H, S, is_dp, &
-    Ef, Qtot, temp)
+    Ef, Qtot, temp, &
+    k_cell, k_displ, &
+    ! Sizes in the end
+    nspecies, no_max, na_u, n_s, nspin, no_u, nnz)
 
   use precision, only: r4, r8
   use io_m, only: open_file, close_file
@@ -101,27 +143,34 @@ subroutine write_hsx1(fname, &
   implicit none
 
   ! Input parameters
-  character(len=*) :: fname
-  integer :: nspecies
-  integer :: na_u, isa(na_u), lasto(na_u)
-  real(r8) :: xa(3,na_u)
-  real(r8) :: ucell(3,3)
-  character(len=*) :: label(nspecies)
-  integer :: Z(nspecies), no_max, no(nspecies)
-  integer, dimension(nspecies, no_max) :: n, l, zeta
-  integer :: n_s, nsc(3), isc_off(3, n_s)
-  integer :: nspin, no_u, no_s, nnz
-  integer :: ncol(no_u), row(nnz)
-  real(r8) :: H(nnz, nspin), S(nnz)
-  logical :: is_dp
-  real(r8) :: Ef, Qtot, temp
+  character(len=*), intent(in) :: fname
+  integer, intent(in) :: nspecies
+  character(len=20), intent(in) :: label(nspecies)
+  real(r8), intent(in) :: Z(nspecies)
+  integer, intent(in) :: no(nspecies)
+  integer, intent(in) :: no_max
+  integer, intent(in), dimension(no_max, nspecies) :: n, l, zeta
+  integer, intent(in) :: na_u, isa(na_u), lasto(na_u)
+  real(r8), intent(in) :: xa(3,na_u), ucell(3,3)
+  integer, intent(in) :: n_s, nsc(3), isc_off(3, n_s)
+  integer, intent(in) :: nspin, no_u, nnz
+  integer, intent(in) :: ncol(no_u), row(nnz)
+  real(r8), intent(in) :: H(nnz, nspin), S(nnz)
+  logical, intent(in) :: is_dp
+  real(r8), intent(in) :: Ef, Qtot, temp
+  integer, intent(in) :: k_cell(3,3)
+  real(r8), intent(in) :: k_displ(3)
 
 ! Define f2py intents
-!f2py intent(in) :: fname, nspecies, na_u, xa, isa, lasto
-!f2py intent(in) :: label, Z, no_max, no, n, l, zeta
-!f2py intent(in) :: ucell, n_s, nsc, isc_off, nspin
-!f2py intent(in) :: no_u, no_s, nnz, ncol, row
+!f2py intent(in) :: fname
+!f2py intent(in) :: na_u, xa, isa, lasto
+!f2py intent(in) :: nspecies, label, Z, no
+!f2py intent(in) :: no_max, n, l, zeta
+!f2py intent(in) :: ucell, n_s, nsc, isc_off
+!f2py intent(in) :: no_u, nnz, nspin
+!f2py intent(in) :: ncol, row
 !f2py intent(in) :: H, S, is_dp, Ef, Qtot, temp
+!f2py intent(in) :: k_cell, k_displ
 
 ! Internal variables and arrays
   integer :: iu, ierr
@@ -130,7 +179,7 @@ subroutine write_hsx1(fname, &
   ! Open file (ensure we start from a clean slate)!
   call open_file(fname, 'write', 'unknown', 'unformatted', iu)
 
-  write(iu, iostat=ierr) 1
+  write(iu, iostat=ierr) 2
   call iostat_update(ierr)
 
 ! Write header information
@@ -148,10 +197,13 @@ subroutine write_hsx1(fname, &
   write(iu, iostat=ierr) (label(is), Z(is), no(is), is=1,nspecies)
   call iostat_update(ierr)
   do is = 1, nspecies
-    write(iu, iostat=ierr) (n(is,i), l(is,i), zeta(is,i), i=1,no(is))
+    write(iu, iostat=ierr) (n(i,is), l(i,is), zeta(i,is), i=1,no(is))
     call iostat_update(ierr)
   end do
 
+  write(iu, iostat=ierr) k_cell, k_displ
+  call iostat_update(ierr)
+
   write(iu, iostat=ierr) ncol
   call iostat_update(ierr)
 
@@ -202,5 +254,4 @@ subroutine write_hsx1(fname, &
 
   call close_file(iu)
 
-end subroutine write_hsx1
-
+end subroutine write_hsx2
diff --git a/src/sisl/io/siesta/binaries.py b/src/sisl/io/siesta/binaries.py
index e50f4644d0..568c33f2f8 100644
--- a/src/sisl/io/siesta/binaries.py
+++ b/src/sisl/io/siesta/binaries.py
@@ -10,11 +10,13 @@
 
 import numpy as np
 
+from sisl.physics.brillouinzone import MonkhorstPack
+
 try:
     from . import _siesta
 
     has_fortran_module = True
-except ImportError:
+except ImportError as _e:
     has_fortran_module = False
 
 import sisl._array as _a
@@ -134,7 +136,7 @@ def get_copy(geom, is_copy):
             f"{cls.__name__}.{method} has non-equal number of supercells, will copy geometry and use file supercell count."
         )
         geom, is_copy = get_copy(geom, is_copy)
-        geom.set_nsc(geom_b.nsc)
+    geom.set_nsc(geom_b.nsc)
 
     # Now for the difficult part.
     # If there is a mismatch in the number of orbitals we will
@@ -788,18 +790,20 @@ def write_density_matrices(self, DM, EDM, Ef: float = 0.0, **kwargs):
 class hsxSileSiesta(SileBinSiesta):
     """Hamiltonian and overlap matrix file
 
-    This file does not contain all information regarding the system.
+    Historically there are several HSX file versions output from Siesta.
 
-    To ensure no errors are being raised one should pass a `Geometry` with
-    correct number of atoms and correct number of supercells.
-    The number of orbitals will be updated in the returned matrices geometry.
+    Since Siesta 5.3, the HSX file format is containing all the information present in
+    the TranSiesta specific TSHS file-format, plus more information.
+    As such, the TSHS file format is deprecated and will not be used anymore.
+    Therefore users should adapt their workflows to be done with HSX
+    files.
 
-    >>> hsx = hsxSileSiesta("siesta.HSX")
-    >>> HS = hsx.read_hamiltonian() # may fail
-    >>> HS = hsx.read_hamiltonian(geometry=<>) # should run correctly if above satisfied
+    Basically, one should never work with HSX files from Siesta <= 4, since they are
+    not information complete.
 
-    Users are adviced to use the `tshsSileSiesta` instead since that correctly contains
-    all information.
+    Notes
+    -----
+    At some point we will deprecate the older HSX file versions (<= 4).
     """
 
     @property
@@ -807,6 +811,11 @@ def version(self) -> int:
         """The version of the file"""
         return int(_siesta.read_hsx_version(self.file))
 
+    @property
+    def _stored_dtype(self):
+        is_dp = _siesta.read_hsx_is_dp(self.file)
+        return {True: np.float64, False: np.float32}[is_dp]
+
     def _xij2system(self, xij, geometry=None, **kwargs):
         """Create a new geometry with *correct* nsc and somewhat correct xyz
 
@@ -1178,7 +1187,7 @@ def read_basis(self, **kwargs) -> Atoms:
             m = 0
             orbs = []
             for n, l, zeta in zip(*n_l_zeta):
-                if old_values != (n, l, zeta):
+                if old_values != (n, l, zeta) or m > l:
                     old_values = (n, l, zeta)
                     m = -l
                 orbs.append(AtomicOrbital(n=n, l=l, m=m, zeta=zeta, R=-1.0))
@@ -1260,17 +1269,19 @@ def _r_lattice_v0(self, **kwargs):
     def _r_lattice_v1(self, **kwargs):
         _, na, _, _, _ = _siesta.read_hsx_sizes(self.file)
 
-        cell, nsc, _, _ = _siesta.read_hsx_geom1(self.file, na)
-
+        cell, nsc, _, _ = _siesta.read_hsx_geom1_2(self.file, na)
+        self._log("lattice v1_2")
         return Lattice(cell.T * _Bohr2Ang, nsc=nsc)
 
+    # Same
+    _r_lattice_v2 = _r_lattice_v1
+
     def read_lattice(self, **kwargs) -> Lattice:
         """Read the lattice from the file
 
         This will always work on new files Siesta >=5, but only sometimes on older
         versions of the HSX file format.
         """
-        self._log("lattice v1")
         return getattr(self, f"_r_lattice_v{self.version}")(**kwargs)
 
     def _r_geometry_v0(self, **kwargs):
@@ -1293,7 +1304,7 @@ def _r_geometry_v0(self, **kwargs):
 
     def _r_geometry_v1(self, **kwargs):
         self._log(
-            "geometry v1 %r, %r",
+            f"geometry v{self.version} %r, %r",
             kwargs.get("atoms", None),
             kwargs.get("geometry", None),
         )
@@ -1303,11 +1314,13 @@ def _r_geometry_v1(self, **kwargs):
         # now read coordinates and cell sizes
         _, na, _, _, _ = _siesta.read_hsx_sizes(self.file)
 
-        cell, nsc, xa, _ = _siesta.read_hsx_geom1(self.file, na)
+        cell, nsc, xa, _ = _siesta.read_hsx_geom1_2(self.file, na)
 
         lattice = Lattice(cell.T * _Bohr2Ang, nsc=nsc)
         return Geometry(xa.T * _Bohr2Ang, atoms, lattice=lattice)
 
+    _r_geometry_v2 = _r_geometry_v1
+
     def read_geometry(self, **kwargs) -> Geometry:
         """Read the geometry from the file
 
@@ -1377,7 +1390,7 @@ def _r_hamiltonian_v0(self, **kwargs):
         H._csr._D[:, spin] = dS[:]
 
         _mat_siesta2sisl(H)
-        H = H.astype(dtype=kwargs.get("dtype"), copy=False)
+        H = H.astype(dtype=kwargs.get("dtype", self._stored_dtype), copy=False)
 
         # Convert the supercells to sisl supercells
         if no_s // no == np.prod(geom.nsc):
@@ -1396,7 +1409,7 @@ def _r_hamiltonian_v1(self, **kwargs):
 
         spin, _, no, no_s, nnz = _siesta.read_hsx_sizes(self.file)
         self._fortran_check("read_hamiltonian", "could not read Hamiltonian sizes.")
-        ncol, col, dH, dS, isc = _siesta.read_hsx_hsx1(self.file, spin, no, no_s, nnz)
+        ncol, col, dH, dS, isc = _siesta.read_hsx_hsx1_2(self.file, spin, no, no_s, nnz)
         col -= 1
         self._fortran_check("read_hamiltonian", "could not read Hamiltonian.")
 
@@ -1423,7 +1436,7 @@ def _r_hamiltonian_v1(self, **kwargs):
         H._csr._D[:, spin] = dS[:]
 
         _mat_siesta2sisl(H)
-        H = H.astype(dtype=kwargs.get("dtype"), copy=False)
+        H = H.astype(dtype=kwargs.get("dtype", self._stored_dtype), copy=False)
 
         # Convert the supercells to sisl supercells
         _csr_from_sc_off(H.geometry, isc.T, H._csr)
@@ -1434,8 +1447,10 @@ def _r_hamiltonian_v1(self, **kwargs):
 
         return H.transpose(spin=False, sort=kwargs.get("sort", True))
 
+    _r_hamiltonian_v2 = _r_hamiltonian_v1
+
     def read_hamiltonian(self, **kwargs) -> Hamiltonian:
-        """Returns the electronic structure from the siesta.TSHS file"""
+        """Returns the electronic structure from the siesta.HSX file"""
         return getattr(self, f"_r_hamiltonian_v{self.version}")(**kwargs)
 
     def _r_overlap_v0(self, **kwargs):
@@ -1456,6 +1471,9 @@ def _r_overlap_v0(self, **kwargs):
                 "inconsistent with HSX file."
             )
 
+        # Get dtype
+        dtype = kwargs.get("dtype", self._stored_dtype)
+
         # Create the Hamiltonian container
         S = Overlap(geom, nnzpr=1)
 
@@ -1466,7 +1484,7 @@ def _r_overlap_v0(self, **kwargs):
         S._csr.col = col.astype(np.int32, copy=False)
         S._csr._nnz = len(col)
 
-        S._csr._D = _a.empty([nnz, 1], dtype=dS.dtype)
+        S._csr._D = _a.empty([nnz, 1], dtype=dtype)
         S._csr._D[:, 0] = dS[:]
 
         # Convert the supercells to sisl supercells
@@ -1483,7 +1501,7 @@ def _r_overlap_v1(self, **kwargs):
         # Now read the sizes used...
         spin, _, no, no_s, nnz = _siesta.read_hsx_sizes(self.file)
         self._fortran_check("read_overlap", "could not read overlap matrix sizes.")
-        ncol, col, dS, isc = _siesta.read_hsx_sx1(self.file, spin, no, no_s, nnz)
+        ncol, col, dS, isc = _siesta.read_hsx_sx1_2(self.file, spin, no, no_s, nnz)
         col -= 1
         self._fortran_check("read_overlap", "could not read overlap matrix.")
 
@@ -1494,6 +1512,9 @@ def _r_overlap_v1(self, **kwargs):
                 "inconsistent with HSX file."
             )
 
+        # Get dtype
+        dtype = kwargs.get("dtype", self._stored_dtype)
+
         # Create the Hamiltonian container
         S = Overlap(geom, nnzpr=1)
 
@@ -1504,7 +1525,7 @@ def _r_overlap_v1(self, **kwargs):
         S._csr.col = col.astype(np.int32, copy=False)
         S._csr._nnz = len(col)
 
-        S._csr._D = _a.empty([nnz, 1], dtype=dS.dtype)
+        S._csr._D = _a.empty([nnz, 1], dtype=dtype)
         S._csr._D[:, 0] = dS[:]
 
         _csr_from_sc_off(S.geometry, isc.T, S._csr)
@@ -1512,10 +1533,155 @@ def _r_overlap_v1(self, **kwargs):
         # not really necessary with Hermitian transposing, but for consistency
         return S.transpose(sort=kwargs.get("sort", True))
 
+    _r_overlap_v2 = _r_overlap_v1
+
     def read_overlap(self, **kwargs) -> Overlap:
-        """Returns the electronic structure from the siesta.TSHS file"""
+        """Returns the electronic structure from the siesta.HSX file"""
         return getattr(self, f"_r_overlap_v{self.version}")(**kwargs)
 
+    def write_hamiltonian(self, H, **kwargs):
+        """Writes the Hamiltonian to a siesta.HSX file"""
+        # we sort below, so no need to do it here
+        # see hsxSileSiesta.read_overlap for .transpose()
+        H = H.transpose(spin=False, sort=False)
+        if H._csr.nnz == 0:
+            raise SileError(
+                f"{self!r}.write_hamiltonian cannot write "
+                "a zero element sparse matrix!"
+            )
+
+        # Default to writing version 2
+        version = kwargs.get("version", 2)
+        if version not in (1, 2):
+            raise NotImplementedError(
+                f"{self!r}.write_hamiltonian only " "supports one of versions [1, 2]."
+            )
+        sort = kwargs.get("sort", True)
+        bz = kwargs.get("bz", None)  # In case the BZ is an argument
+        dtype = kwargs.get("dtype", np.float64)
+        if dtype not in (np.float32, np.float64):
+            raise ValueError(
+                f"{self!r}.write_hamiltonian only "
+                "supports float32 or float64 in the HSX file format."
+            )
+
+        # Convert to siesta CSR
+        _csr_to_siesta(H.geometry, H._csr, diag=True)
+        H.finalize(sort=sort)
+
+        # we have the required dtype above.
+        # However, the fortran interface requires a double, so we have to
+        # use a bool to signal the stored data.
+        H = H.astype(dtype=np.float64, copy=False)
+        _mat_sisl2siesta(H)
+        csr = H._csr
+
+        # Extract the data to pass to the fortran routine
+        cell = H.geometry.cell
+        xyz = H.geometry.xyz
+
+        # Get H and S
+        if H.orthogonal:
+            self._log("writing orthogonal hamiltonian")
+            s = csr.copy(dims=0)
+            s.empty(keep_nnz=True)
+            s += s.diags(1)
+            missing_diags = s.nnz - csr.nnz
+            if missing_diags:
+                # This should never happen as _csr_to_siesta should ensure
+                # that the diagonal entries always exists.
+                # Hence it gets changed before writing.
+                # Not compeletely optimal, but for now this is OK
+                raise SislError(
+                    f"{self.__class__.__name__}.write_hamiltonian "
+                    "The diagonal elements of your orthogonal Hamiltonian "
+                    f"have not been defined. Got {len(csr) - missing_diags} elements, expected {len(csr)}."
+                )
+            h = csr._D
+            s = s._D
+        else:
+            self._log("writing non-orthogonal hamiltonian")
+            h = csr._D[:, : H.S_idx]
+            s = csr._D[:, H.S_idx]
+
+        # Get shorter variants
+        nsc = H.geometry.nsc[:].astype(np.int32)
+        isc = _siesta.siesta_sc_off(*nsc)
+
+        # Correct arguments depending on BZ
+        kcell = _a.zerosi([3, 3])
+        kdispl = _a.zerosd([3])
+        if isinstance(bz, MonkhorstPack):
+            kcell = np.diag(bz._diag)
+            kdispl = bz._displ
+        else:
+            # TODO, perhaps do something based on the type of BZ
+            for i in range(3):
+                if nsc[i] > 1:
+                    kcell[i, i] = 2
+                else:
+                    kcell[i, i] = 1
+
+        geom = H.geometry
+        atoms = geom.atoms
+
+        # Populate arguments
+        args = []
+
+        def func(atom_):
+            atom, _ = atom_
+            return f"{atom.tag:20s}"[:20], atom.Z, len(atom)
+
+        args.extend(list(map(list, zip(*map(func, atoms.iter(species=True))))))
+
+        # Next values
+        no_max = (atoms.lasto - atoms.firsto[:-1]).max() + 1
+        # Atoms require 'n' and 'zeta' to be >=1
+        n = _a.zerosi([no_max, atoms.nspecies], order="F") + 1
+        l = np.zeros_like(n)
+        zeta = np.zeros_like(n) + 1
+        try:
+            for i_s, (atom, _) in enumerate(atoms.iter(species=True)):
+                for io, orb in enumerate(atom):
+                    n[io, i_s] = orb.n
+                    l[io, i_s] = orb.l
+                    zeta[io, i_s] = orb.zeta
+        except:
+            pass  # should we warn?
+
+        args.extend([n, l, zeta])
+
+        args.extend(
+            [
+                _toF(xyz.T, np.float64, 1 / _Bohr2Ang),
+                atoms.species + 1,
+                atoms.firsto[:-1],
+            ]
+        )
+
+        # Now cell etc.
+        args.append(_toF(cell.T, np.float64, 1 / _Bohr2Ang))
+        args.append(nsc)
+        args.append(isc)  # already in correct shape!
+
+        # Now the matrices
+        args.extend([csr.ncol, csr.col + 1])
+        args.append(_toF(h, np.float64, _eV2Ry))
+        args.append(_toF(s, np.float64))
+        args.append(dtype == np.float64)
+        args.extend([0.0, 0.0, 0.0])
+
+        if version == 2:
+            args.append(kcell)
+            args.append(kdispl)
+
+        # see hsxSileSiesta.read_lattice for .T
+        write_hsx = getattr(_siesta, f"write_hsx{version}")
+        write_hsx(self.file, *args)
+        self._fortran_check(
+            "write_hamiltonian", "could not write Hamiltonian and overlap matrix."
+        )
+
 
 @set_module("sisl.io.siesta")
 class wfsxSileSiesta(SileBinSiesta):
diff --git a/src/sisl/io/siesta/tests/test_hsx.py b/src/sisl/io/siesta/tests/test_hsx.py
index ff059f54ea..d282680ee2 100644
--- a/src/sisl/io/siesta/tests/test_hsx.py
+++ b/src/sisl/io/siesta/tests/test_hsx.py
@@ -156,6 +156,45 @@ def test_si_pdos_kgrid_hsx_1_same_tshs(sisl_files, sisl_tmp):
     assert np.allclose(gx.nsc, gt.nsc)
 
 
+def test_hsx_si_pdos_kgrid_tofromnc(sisl_files, sisl_tmp):
+    pytest.importorskip("netCDF4")
+    si = sisl.get_sile(sisl_files("siesta", "Si_pdos_k", "Si_pdos.TSHS"))
+    HS1 = si.read_hamiltonian()
+    f = sisl_tmp("tmp.HSX")
+    fnc = sisl_tmp("tmp.nc")
+
+    HS1.write(f)
+    HS1.write(fnc)
+
+    HS2 = sisl.get_sile(f).read_hamiltonian()
+    HS2nc = sisl.get_sile(fnc).read_hamiltonian()
+    assert HS1._csr.spsame(HS2._csr)
+    assert HS1._csr.spsame(HS2nc._csr)
+    HS1.finalize()
+    HS2.finalize()
+    assert np.allclose(HS1._csr._D, HS2._csr._D)
+    HS2nc.finalize()
+    assert np.allclose(HS1._csr._D, HS2nc._csr._D)
+
+
+@pytest.mark.parametrize("version", [1, 2])
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_hsx_si_pdos_kgrid_dtype(sisl_files, sisl_tmp, version, dtype):
+    pytest.importorskip("netCDF4")
+    si = sisl.get_sile(sisl_files("siesta", "Si_pdos_k", "Si_pdos.TSHS"))
+    HS = si.read_hamiltonian()
+    f = sisl_tmp("tmp.1.HSX")
+
+    HS.write(f, version=version, dtype=dtype)
+
+    HS1 = sisl.get_sile(f).read_hamiltonian()
+    assert HS1.dtype == dtype
+    assert HS._csr.spsame(HS1._csr)
+    HS.finalize()
+    HS1.finalize()
+    assert np.allclose(HS._csr._D, HS1._csr._D)
+
+
 @pytest.mark.parametrize("order", ["HSX", "TSHS"])
 def test_srti03_H_eig(sisl_files, sisl_tmp, order):
     fdf = sisl.get_sile(sisl_files("siesta", "SrTiO3", "non-collinear", "SrTiO3.fdf"))
diff --git a/src/sisl/io/sile.py b/src/sisl/io/sile.py
index de9c41012f..3a73c9c180 100644
--- a/src/sisl/io/sile.py
+++ b/src/sisl/io/sile.py
@@ -363,8 +363,13 @@ def get_eligibles(end, rules):
                     f"Cannot determine the exact Sile requested, multiple hits: {tuple(e.cls.__name__ for e in eligibles)}"
                 )
 
+        # Print-out error on which extensions it tried (and full filename)
+        if len(end_list) == 1:
+            ext_list = end_list
+        else:
+            ext_list = end_list[1:]
         raise NotImplementedError(
-            f"Sile for file '{filename}' could not be found, "
+            f"Sile for file '{filename}' ({ext_list}) could not be found, "
             "possibly the file has not been implemented."
         )
 

From 850896082ed0c06603393e96096a12316748be54 Mon Sep 17 00:00:00 2001
From: Nick Papior <nickpapior@gmail.com>
Date: Wed, 15 Jan 2025 15:29:10 +0100
Subject: [PATCH 2/2] removed BaseException and made BTD coherent

Signed-off-by: Nick Papior <nickpapior@gmail.com>
---
 src/sisl/io/siesta/binaries.py | 2 +-
 src/sisl/io/siesta/stdout.py   | 4 ++--
 src/sisl_toolbox/btd/_btd.py   | 9 ++++++++-
 tools/changelog.py             | 2 +-
 4 files changed, 12 insertions(+), 5 deletions(-)

diff --git a/src/sisl/io/siesta/binaries.py b/src/sisl/io/siesta/binaries.py
index 568c33f2f8..357e3df782 100644
--- a/src/sisl/io/siesta/binaries.py
+++ b/src/sisl/io/siesta/binaries.py
@@ -1646,7 +1646,7 @@ def func(atom_):
                     n[io, i_s] = orb.n
                     l[io, i_s] = orb.l
                     zeta[io, i_s] = orb.zeta
-        except:
+        except Exception:
             pass  # should we warn?
 
         args.extend([n, l, zeta])
diff --git a/src/sisl/io/siesta/stdout.py b/src/sisl/io/siesta/stdout.py
index 27e8852d47..8a7ff45485 100644
--- a/src/sisl/io/siesta/stdout.py
+++ b/src/sisl/io/siesta/stdout.py
@@ -119,7 +119,7 @@ def _parse_version(attr, instance, match):
     version, *spec = opt.split("-", maxsplit=1)
     try:
         version = tuple(int(v) for v in version.split("."))
-    except BaseException:
+    except Exception:
         version = (0, 0, 0)
 
     # Convert version to a tuple
@@ -1431,7 +1431,7 @@ def try_parse_int(s):
                             "Mulliken Atomic Populations",
                             "Mulliken Net Atomic Populations",
                         ]
-                except BaseException:
+                except Exception:
                     pass
 
         else:
diff --git a/src/sisl_toolbox/btd/_btd.py b/src/sisl_toolbox/btd/_btd.py
index f5f40c5886..0eb8ab47d8 100644
--- a/src/sisl_toolbox/btd/_btd.py
+++ b/src/sisl_toolbox/btd/_btd.py
@@ -615,7 +615,14 @@ def from_fdf(cls, fdf, prefix="TBT", use_tbt_se=False, eta=None):
         is_tbtrans = prefix.upper() == "TBT"
 
         # Read the device H, only valid for TBT stuff
-        Hdev = si.get_sile(fdf.get("TBT.HS", f"{slabel}.TSHS")).read_hamiltonian()
+        for hs_ext in ("TS.HSX", "TSHS", "HSX"):
+            if Path(f"{slabel}.{hs_ext}").exists():
+                # choose a sane default (if it exists!)
+                hs_default = f"{slabel}.{hs_ext}"
+                break
+        else:
+            hs_default = f"{slabel}.TSHS"
+        Hdev = si.get_sile(fdf.get("TBT.HS", hs_default)).read_hamiltonian()
 
         def get_line(line):
             """Parse lines in the %block constructs of fdf's"""
diff --git a/tools/changelog.py b/tools/changelog.py
index 2809d24cb9..0fe138421f 100644
--- a/tools/changelog.py
+++ b/tools/changelog.py
@@ -123,7 +123,7 @@ def get_pull_requests(repo, revision_range):
 
         try:
             prs.append(repo.get_pull(n))
-        except BaseException:
+        except Exception:
             pass
 
     return prs