bug: fixed final errors in HSXv1

This should cover most things.

Also the fdf is now defaulting to also reading from
the HSX file

sisl/io/siesta/_src/hsx_read.f90

@@ -274,7 +274,7 @@ subroutine read_hsx_ef1(fname, Ef)
 
   ! Internal variables and arrays
   integer :: version
-  real(dp) :: ucell(3,3), qtot, temp
+  real(dp) :: ucell(3,3)
   integer :: iu, ierr
 
   call open_file(fname, 'read', 'old', 'unformatted', iu)
@@ -290,7 +290,7 @@ subroutine read_hsx_ef1(fname, Ef)
   read(iu, iostat=ierr) !na_u, no_u, nspin, nspecies, nsc
   call iostat_update(ierr)
 
-  read(iu, iostat=ierr) ucell, Ef, qtot, temp
+  read(iu, iostat=ierr) ucell, Ef! , qtot, temp
   call iostat_update(ierr)
 
   call close_file(iu)
@@ -750,11 +750,9 @@ subroutine read_hsx_sx1(fname, nspin, no_u, no_s, maxnh, &
   end do
 
   ! Read Hamiltonian
-  do is = 1 , nspin
-    do ih = 1 , no_u
-      read(iu, iostat=ierr) !
-      call iostat_update(ierr)
-    end do
+  do is = 1 , nspin * no_u
+    read(iu, iostat=ierr) !
+    call iostat_update(ierr)
   end do
 
   if ( is_dp ) then
@@ -807,7 +805,7 @@ subroutine read_hsx_geom1(fname, na_u, cell, nsc, xa, lasto)
 !f2py intent(out) :: cell, nsc, isa, xa, lasto
 
   integer :: iu, ierr, version
-  integer :: lna_u, lno_u, lnspin, nspecies
+  integer :: lna_u, no_u, nspin, nspecies
   integer, allocatable :: isc(:,:), isa(:)
 
   call open_file(fname, 'read', 'old', 'unformatted', iu)
@@ -820,7 +818,10 @@ subroutine read_hsx_geom1(fname, na_u, cell, nsc, xa, lasto)
     return
   end if
 
-  read(iu, iostat=ierr) lna_u, lno_u, lnspin, nspecies, nsc
+  read(iu, iostat=ierr) ! is_dp
+  call iostat_update(ierr)
+
+  read(iu, iostat=ierr) lna_u, no_u, nspin, nspecies, nsc
   call iostat_update(ierr)
   if ( lna_u /= na_u ) call iostat_update(-6)
 

sisl/io/siesta/binaries.py

@@ -1035,8 +1035,8 @@ def _r_geometry_v1(self, **kwargs):
 
         cell, nsc, xa, _ = _siesta.read_hsx_geom1(self.file, na)
 
-        sc = SuperCell(cell, nsc=nsc)
-        return Geometry(xa, atoms, sc=sc)
+        sc = SuperCell(cell.T * _Bohr2Ang, nsc=nsc)
+        return Geometry(xa.T * _Bohr2Ang, atoms, sc=sc)
 
     def read_geometry(self, **kwargs):
         """ Read the geometry from the file
@@ -1072,7 +1072,7 @@ def _r_hamiltonian_v0(self, **kwargs):
         geom = self._xij2system(xij, kwargs.get("geometry", kwargs.get("geom", None)))
 
         if geom.no != no or geom.no_s != no_s:
-            raise SileError(f"{str(self)}.read_hamiltonian could not use the "
+            raise SileError(f"{self!s}.read_hamiltonian could not use the "
                             "passed geometry as the number of atoms or orbitals is "
                             "inconsistent with HSX file.")
 
@@ -1100,14 +1100,16 @@ def _r_hamiltonian_v0(self, **kwargs):
 
     def _r_hamiltonian_v1(self, **kwargs):
         # Now read the sizes used...
+        geom = self.read_geometry(**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)
         col -= 1
         self._fortran_check("read_hamiltonian", "could not read Hamiltonian.")
 
         if geom.no != no or geom.no_s != no_s:
-            raise SileError(f"{str(self)}.read_hamiltonian could not use the "
+            raise SileError(f"{self!s}.read_hamiltonian could not use the "
                             "passed geometry as the number of atoms or orbitals is "
                             "inconsistent with HSX file.")
 
@@ -1116,7 +1118,8 @@ def _r_hamiltonian_v1(self, **kwargs):
 
         # Create the new sparse matrix
         H._csr.ncol = ncol.astype(np.int32, copy=False)
-        H._csr.ptr = ptr
+        H._csr.ptr = _ncol_to_indptr(ncol).astype(np.int32, copy=False)
+
         # Correct fortran indices
         H._csr.col = col.astype(np.int32, copy=False)
         H._csr._nnz = len(col)
@@ -1150,7 +1153,7 @@ def _r_overlap_v0(self, **kwargs):
         self._fortran_check("read_overlap", "could not read overlap matrix.")
 
         if geom.no != no or geom.no_s != no_s:
-            raise SileError(f"{str(self)}.read_overlap could not use the "
+            raise SileError(f"{self!s}.read_overlap could not use the "
                             "passed geometry as the number of atoms or orbitals is "
                             "inconsistent with HSX file.")
 
@@ -1186,7 +1189,7 @@ def _r_overlap_v1(self, **kwargs):
         self._fortran_check("read_overlap", "could not read overlap matrix.")
 
         if geom.no != no or geom.no_s != no_s:
-            raise SileError(f"{str(self)}.read_overlap could not use the "
+            raise SileError(f"{self!s}.read_overlap could not use the "
                             "passed geometry as the number of atoms or orbitals is "
                             "inconsistent with HSX file.")
 

sisl/io/siesta/fdf.py

@@ -1258,7 +1258,7 @@ def read_geometry(self, output=False, *args, **kwargs):
             the fdf file).
         order: list of str, optional
             the order of which to try and read the geometry.
-            By default this is ``["XV", "nc", "fdf", "TSHS", "STRUCT"]`` if `output` is true
+            By default this is ``["XV", "nc", "fdf", "TSHS", "HSX", "STRUCT"]`` if `output` is true
             If `order` is present `output` is disregarded.
 
         Examples
@@ -1275,7 +1275,7 @@ def read_geometry(self, output=False, *args, **kwargs):
         # code to correct.
         ##
         if output:
-            order = _listify_str(kwargs.pop("order", ["XV", "nc", "fdf", "TSHS", "STRUCT"]))
+            order = _listify_str(kwargs.pop("order", ["XV", "nc", "fdf", "TSHS", "HSX", "STRUCT"]))
         else:
             order = _listify_str(kwargs.pop("order", ["fdf"]))
         for f in order:
@@ -1345,6 +1345,15 @@ def _r_geometry_tshs(self):
             return tshsSileSiesta(f).read_geometry(geometry=self.read_geometry(False))
         return None
 
+    def _r_geometry_hsx(self):
+        # Read geometry from <>.TSHS file
+        f = self.dir_file(self.get("SystemLabel", default="siesta") + ".HSX")
+        _track_file(self._r_geometry_hsx, f)
+        if f.is_file():
+            # Default to a geometry with the correct atomic numbers etc.
+            return hsxSileSiesta(f).read_geometry(geometry=self.read_geometry(False))
+        return None
+
     def _r_geometry_fdf(self, *args, **kwargs):
         """ Returns Geometry object from the FDF file
 
@@ -2066,7 +2075,7 @@ def _r_hamiltonian_tshs(self, *args, **kwargs):
         if f.is_file():
             if "geometry" not in kwargs:
                 # to ensure we get the correct orbital count
-                kwargs["geometry"] = self.read_geometry(True, order=["nc", "TSHS"])
+                kwargs["geometry"] = self.read_geometry(True, order=["nc", "TSHS", "HSX"])
             H = tshsSileSiesta(f).read_hamiltonian(*args, **kwargs)
         return H
 
@@ -2076,9 +2085,10 @@ def _r_hamiltonian_hsx(self, *args, **kwargs):
         _track_file(self._r_hamiltonian_hsx, f)
         H = None
         if f.is_file():
-            if "geometry" not in kwargs:
-                # to ensure we get the correct orbital count
-                kwargs["geometry"] = self.read_geometry(True, order=["nc", "TSHS", "fdf"])
+            if hsxSileSiesta(f).version == 0:
+                if "geometry" not in kwargs:
+                    # to ensure we get the correct orbital count
+                    kwargs["geometry"] = self.read_geometry(True, order=["nc", "TSHS", "fdf"])
             H = hsxSileSiesta(f).read_hamiltonian(*args, **kwargs)
             Ef = self.read_fermi_level()
             if Ef is None:
@@ -2110,12 +2120,7 @@ def label_file(suffix):
             return self.dir_file(f_label).with_suffix(suffix)
 
         # The default on all sub-parsers are the retrieval and setting
-
-        d = {
-            "_fdf": self,
-            "_fdf_first": True,
-        }
-        namespace = default_namespace(**d)
+        namespace = default_namespace(_fdf=self, _fdf_first=True)
 
         ep = sp.add_parser("edit",
                            help="Change or read and print data from the fdf file")

sisl/io/siesta/tests/test_hsx.py

@@ -111,5 +111,19 @@ def test_si_pdos_kgrid_hsx_versions_s(sisl_files, sisl_tmp):
     assert HS0._csr.spsame(S1._csr)
     assert S0._csr.spsame(HS1._csr)
 
-    assert np.allclose(HS0._csr._D[:, HS0.S_idx], S0._csr._D)
-    assert np.allclose(HS1._csr._D[:, HS1.S_idx], S1._csr._D)
+    assert np.allclose(HS0._csr._D[:, HS0.S_idx], S0._csr._D[:, 0])
+    assert np.allclose(HS1._csr._D[:, HS1.S_idx], S1._csr._D[:, 0])
+
+
+def test_si_pdos_kgrid_hsx_1_same_tshs(sisl_files, sisl_tmp):
+    HSX = sisl.get_sile(sisl_files(_dir, "si_pdos_kgrid.1.HSX"))
+    TSHS = sisl.get_sile(sisl_files(_dir, "si_pdos_kgrid.TSHS"))
+
+    HSX = HSX.read_hamiltonian()
+    TSHS = TSHS.read_hamiltonian()
+
+    gx = HSX.geometry
+    gt = TSHS.geometry
+    assert np.allclose(gx.sc.cell, gt.sc.cell)
+    assert np.allclose(gx.xyz, gt.xyz)
+    assert np.allclose(gx.nsc, gt.nsc)