changed HSX file reads

When reading the Hamiltonian, it now shifts
the Fermi-level to 0. This is to be compatible
with other tools, and if users reads directly
from the file.

The fdf file reads still performs normally.

Added Sile._log to perform loggings. It
is a simple wrapper for the logging module.

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

CHANGELOG.md

@@ -95,6 +95,7 @@ we hit release version 1.0.0.
 - removed `Selector` and `TimeSelector`, they were never used internally
 
 ### Changed
+- `hsxSileSiesta.read_hamiltonian` now implicitly shifts Fermi-level to 0 (for newer HSX versions)
 - deprecated `periodic` to `axes` argument in `BrillouinZone.volume`
 - changed `Eigenmode.displacement` shape, please read the documentation
 - bumped minimal Python version to 3.9, #640

src/sisl/io/siesta/_src/hsx_read.f90

@@ -222,11 +222,7 @@ subroutine read_hsx_ef(fname, Ef)
 
   call read_hsx_version(fname, version)
 
-  if ( version == 0 ) then ! old
-
-    call read_hsx_ef0(fname, Ef)
-
-  else if ( version == 1 ) then
+  if ( version == 1 ) then
 
     call read_hsx_ef1(fname, Ef)
 

src/sisl/io/siesta/binaries.py

@@ -409,9 +409,11 @@ def read_hamiltonian(self, geometry=None, **kwargs) -> Hamiltonian:
         H._csr._nnz = len(col)
 
         if orthogonal:
+            self._log("read orthogonal hamiltonian")
             H._csr._D = _a.emptyd([nnz, spin])
             H._csr._D[:, :] = dH[:, :] * _Ry2eV
         else:
+            self._log("read non-orthogonal hamiltonian")
             H._csr._D = _a.emptyd([nnz, spin + 1])
             H._csr._D[:, :spin] = dH[:, :] * _Ry2eV
             H._csr._D[:, spin] = dS[:]
@@ -428,7 +430,7 @@ def read_hamiltonian(self, geometry=None, **kwargs) -> Hamiltonian:
             print(f"Number of orbitals: {no}")
             print(idx)
             raise SileError(
-                f"{self!s}.read_hamiltonian could not assert "
+                f"{self!r}.read_hamiltonian could not assert "
                 "the supercell connections in the primary unit-cell."
             )
 
@@ -443,7 +445,7 @@ def write_hamiltonian(self, H, **kwargs):
         csr = H.transpose(spin=False, sort=False)._csr
         if csr.nnz == 0:
             raise SileError(
-                f"{self!s}.write_hamiltonian cannot write "
+                f"{self!r}.write_hamiltonian cannot write "
                 "a zero element sparse matrix!"
             )
 
@@ -460,6 +462,7 @@ def write_hamiltonian(self, H, **kwargs):
 
         # 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)
@@ -477,6 +480,7 @@ def write_hamiltonian(self, H, **kwargs):
             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]
 
@@ -542,7 +546,7 @@ def read_density_matrix(self, **kwargs) -> DensityMatrix:
 
         if geom.no != no:
             raise SileError(
-                f"{self!s}.read_density_matrix could not use the "
+                f"{self!r}.read_density_matrix could not use the "
                 "passed geometry as the number of atoms or orbitals is "
                 "inconsistent with DM file."
             )
@@ -568,7 +572,7 @@ def read_density_matrix(self, **kwargs) -> DensityMatrix:
         if nsc[0] != 0 or geom.no_s >= col.max():
             _csr_from_siesta(geom, DM._csr)
         else:
-            warn(f"{self!s}.read_density_matrix may result in a wrong sparse pattern!")
+            warn(f"{self!r}.read_density_matrix may result in a wrong sparse pattern!")
 
         DM = DM.transpose(spin=False, sort=kwargs.get("sort", True))
         _add_overlap(
@@ -584,7 +588,7 @@ def write_density_matrix(self, DM, **kwargs):
         # This ensures that we don"t have any *empty* elements
         if csr.nnz == 0:
             raise SileError(
-                f"{self!s}.write_density_matrix cannot write "
+                f"{self!r}.write_density_matrix cannot write "
                 "a zero element sparse matrix!"
             )
 
@@ -648,7 +652,7 @@ def read_energy_density_matrix(self, **kwargs) -> EnergyDensityMatrix:
 
         if geom.no != no:
             raise SileError(
-                f"{self!s}.read_energy_density_matrix could "
+                f"{self!r}.read_energy_density_matrix could "
                 "not use the passed geometry as the number of atoms or orbitals "
                 "is inconsistent with DM file."
             )
@@ -677,7 +681,7 @@ def read_energy_density_matrix(self, **kwargs) -> EnergyDensityMatrix:
             _csr_from_siesta(geom, EDM._csr)
         else:
             warn(
-                f"{self!s}.read_energy_density_matrix may result in a wrong sparse pattern!"
+                f"{self!r}.read_energy_density_matrix may result in a wrong sparse pattern!"
             )
 
         EDM = EDM.transpose(spin=False, sort=kwargs.get("sort", True))
@@ -719,7 +723,7 @@ def write_density_matrices(self, DM, EDM, Ef=0.0, **kwargs):
 
         if DMcsr.nnz == 0:
             raise SileError(
-                f"{self!s}.write_density_matrices cannot write "
+                f"{self!r}.write_density_matrices cannot write "
                 "a zero element sparse matrix!"
             )
 
@@ -736,7 +740,7 @@ def write_density_matrices(self, DM, EDM, Ef=0.0, **kwargs):
             np.allclose(DMcsr.ncol, EDMcsr.ncol) and np.allclose(DMcsr.col, EDMcsr.col)
         ):
             raise ValueError(
-                f"{self!s}.write_density_matrices got non compatible "
+                f"{self!r}.write_density_matrices got non compatible "
                 "DM and EDM matrices."
             )
 
@@ -1094,6 +1098,9 @@ def conv(orbs, atm):
 
     def read_basis(self, **kwargs) -> Atoms:
         """Reads basis set and geometry information from the HSX file"""
+        self._log(
+            "basis %r %r", kwargs.get("atoms", None), kwargs.get("geometry", None)
+        )
         # Now read the sizes used...
         no, na, nspecies = _siesta.read_hsx_species_sizes(self.file)
         self._fortran_check("read_basis", "could not read specie sizes.")
@@ -1179,25 +1186,42 @@ def read_basis(self, **kwargs) -> Atoms:
         def get_best(aF, aI):
             if len(aF) != len(aI):
                 # the file has the correct number of orbitals
+                self._log("basis file atom=%r better than input %r", aF, aI)
                 return aF
 
             if not isinstance(aF, AtomUnknown):
                 if aF.Z != aI.Z:
+                    self._log(
+                        "basis file atom=%r has different Z than input %r", aF, aI
+                    )
                     return aF
 
             # check for orbitals being atomicorbital
             for orb in aI:
                 if not isinstance(orb, AtomicOrbital):
+                    self._log(
+                        "basis file, input atom=%r does not use AtomicOrbital", aI
+                    )
                     return aF
 
             for oF, oI in zip(aF, aI):
-                for prop in ("n", "l", "m", "zeta", "P"):
+                # we can't check polarization, not stored in the HSX file format
+                for prop in ("n", "l", "m", "zeta"):
                     if getattr(oF, prop) != getattr(oI, prop):
+                        self._log(
+                            "basis file atom, orbital %r,%r property %s cannot be matched with input atom, orbital %r,%r",
+                            aF,
+                            oF,
+                            prop,
+                            aI,
+                            oI,
+                        )
                         return aF
 
             # the atoms are the same, so we select the input atom
             # since it likely contains the spherical functions
             # and the charge
+            self._log("basis file, choose input atom=%r over file atom %r", aI, aF)
             return aI
 
         atoms = Atoms(map(get_best, atoms, base))
@@ -1231,11 +1255,16 @@ def read_lattice(self, **kwargs) -> Lattice:
         This will always work on new files Siesta >=5, but only sometimes on older
         versions of the HSX file format.
         """
-        version = _siesta.read_hsx_version(self.file)
-        return getattr(self, f"_r_lattice_v{version}")(**kwargs)
+        self._log("lattice v1")
+        return getattr(self, f"_r_lattice_v{self.version}")(**kwargs)
 
     def _r_geometry_v0(self, **kwargs):
         """Read the geometry from the old file version"""
+        self._log(
+            "geometry v0 %r, %r",
+            kwargs.get("atoms", None),
+            kwargs.get("geometry", None),
+        )
         spin, _, no, no_s, nnz = _siesta.read_hsx_sizes(self.file)
         self._fortran_check("read_geometry", "could not read geometry sizes.")
         ncol, col, _, _, dxij = _siesta.read_hsx_hsx0(self.file, spin, no, no_s, nnz)
@@ -1248,6 +1277,11 @@ def _r_geometry_v0(self, **kwargs):
         return geom
 
     def _r_geometry_v1(self, **kwargs):
+        self._log(
+            "geometry v1 %r, %r",
+            kwargs.get("atoms", None),
+            kwargs.get("geometry", None),
+        )
         # first read the atoms object
         atoms = self.read_basis(**kwargs)
 
@@ -1265,20 +1299,19 @@ def read_geometry(self, **kwargs) -> Geometry:
         This will always work on new files Siesta >=5, but only sometimes on older
         versions of the HSX file format.
         """
-        version = _siesta.read_hsx_version(self.file)
-        return getattr(self, f"_r_geometry_v{version}")(**kwargs)
+        return getattr(self, f"_r_geometry_v{self.version}")(**kwargs)
 
     def read_fermi_level(self, **kwargs) -> float:
         """Reads the fermi level in the file
 
         Only valid for files created by Siesta >=5.
         """
+        if self.version == 0:
+            # It isn't implemented in the old file format.
+            warn(f"{self!r} does not contain the fermi-level (too old version file).")
+            return None
         Ef = _siesta.read_hsx_ef(self.file)
-        msg = self._fortran_check(
-            "read_fermi_level", "could not read Fermi-level", ret_msg=True
-        )
-        if msg:
-            warn(msg)
+        self._fortran_check("read_fermi_level", "could not read Fermi-level")
         return Ef * _Ry2eV
 
     def _r_hamiltonian_v0(self, **kwargs):
@@ -1293,7 +1326,7 @@ def _r_hamiltonian_v0(self, **kwargs):
 
         if geom.no != no or geom.no_s != no_s:
             raise SileError(
-                f"{self!s}.read_hamiltonian could not use the "
+                f"{self!r}.read_hamiltonian could not use the "
                 "passed geometry as the number of atoms or orbitals is "
                 "inconsistent with HSX file."
             )
@@ -1332,7 +1365,7 @@ def _r_hamiltonian_v1(self, **kwargs):
 
         if geom.no != no or geom.no_s != no_s:
             raise SileError(
-                f"{self!s}.read_hamiltonian could not use the "
+                f"{self!r}.read_hamiltonian could not use the "
                 "passed geometry as the number of atoms or orbitals is "
                 "inconsistent with HSX file."
             )
@@ -1357,12 +1390,15 @@ def _r_hamiltonian_v1(self, **kwargs):
         # Convert the supercells to sisl supercells
         _csr_from_sc_off(H.geometry, isc.T, H._csr)
 
+        # Correct the fermi-level here
+        Ef = self.read_fermi_level()
+        H.shift(-Ef)
+
         return H.transpose(spin=False, sort=kwargs.get("sort", True))
 
     def read_hamiltonian(self, **kwargs) -> Hamiltonian:
         """Returns the electronic structure from the siesta.TSHS file"""
-        version = _siesta.read_hsx_version(self.file)
-        return getattr(self, f"_r_hamiltonian_v{version}")(**kwargs)
+        return getattr(self, f"_r_hamiltonian_v{self.version}")(**kwargs)
 
     def _r_overlap_v0(self, **kwargs):
         """Returns the overlap matrix from the siesta.HSX file"""
@@ -1377,7 +1413,7 @@ def _r_overlap_v0(self, **kwargs):
 
         if geom.no != no or geom.no_s != no_s:
             raise SileError(
-                f"{self!s}.read_overlap could not use the "
+                f"{self!r}.read_overlap could not use the "
                 "passed geometry as the number of atoms or orbitals is "
                 "inconsistent with HSX file."
             )
@@ -1415,7 +1451,7 @@ def _r_overlap_v1(self, **kwargs):
 
         if geom.no != no or geom.no_s != no_s:
             raise SileError(
-                f"{self!s}.read_overlap could not use the "
+                f"{self!r}.read_overlap could not use the "
                 "passed geometry as the number of atoms or orbitals is "
                 "inconsistent with HSX file."
             )
@@ -1440,8 +1476,7 @@ def _r_overlap_v1(self, **kwargs):
 
     def read_overlap(self, **kwargs) -> Overlap:
         """Returns the electronic structure from the siesta.TSHS file"""
-        version = _siesta.read_hsx_version(self.file)
-        return getattr(self, f"_r_overlap_v{version}")(**kwargs)
+        return getattr(self, f"_r_overlap_v{self.version}")(**kwargs)
 
 
 @set_module("sisl.io.siesta")

src/sisl/io/siesta/fdf.py

@@ -161,7 +161,7 @@ def _pushfile(self, f):
             self.fh = gzip.open(self.dir_file(f"{f}.gz"), mode="rt")
         else:
             warn(
-                f"{self!s} is trying to include file: {f} but the file seems not to exist? Will disregard file!"
+                f"{self!r} is trying to include file: {f} but the file seems not to exist? Will disregard file!"
             )
 
     def _popfile(self):
@@ -666,7 +666,7 @@ def write_block(atoms, append, write_block):
                 idx = list2str(geometry.names[n] + 1).replace("-", " -- ")
                 if len(idx) > 200:
                     info(
-                        f"{self!s}.write_geometry will not write the constraints for {n} (too long line)."
+                        f"{self!r}.write_geometry will not write the constraints for {n} (too long line)."
                     )
                 else:
                     _write_block = write_block(idx, append, _write_block)
@@ -832,7 +832,7 @@ def _r_lattice_xv(self, *args, **kwargs):
 
     def _r_lattice_struct(self, *args, **kwargs):
         """Returns `Lattice` object from the STRUCT files"""
-        for end in ["STRUCT_NEXT_ITER", "STRUCT_OUT", "STRUCT_IN"]:
+        for end in ("STRUCT_NEXT_ITER", "STRUCT_OUT", "STRUCT_IN"):
             f = self.dir_file(self.get("SystemLabel", default="siesta") + f".{end}")
             _track_file(self._r_lattice_struct, f)
             if f.is_file():
@@ -1228,7 +1228,7 @@ def _dynamical_matrix_from_fc(self, geom, FC, FC_atoms, *args, **kwargs):
                 daxis = geom_tile.xyz[:, ax] - geom.xyz[:, ax]
                 if not np.allclose(daxis, daxis[0], rtol=0.0, atol=0.01):
                     raise SislError(
-                        f"{self!s}.read_dynamical_matrix(FC) could "
+                        f"{self!r}.read_dynamical_matrix(FC) could "
                         "not figure out the tiling method for the supercell"
                     )
 
@@ -1858,7 +1858,7 @@ def _r_basis_ion(self):
         if atoms_species:
             return Atoms([atoms[spc] for spc in atoms_species])
         warn(
-            f"{self!s} does not contain the AtomicCoordinatesAndAtomicSpecies block, basis set definition may not contain all atoms."
+            f"{self!r} does not contain the AtomicCoordinatesAndAtomicSpecies block, basis set definition may not contain all atoms."
         )
         return Atoms(atoms)
 
@@ -1924,7 +1924,7 @@ def _r_basis_fdf(self):
             return Atoms([atoms[spc] for spc in atoms_species])
 
         warn(
-            f"{self!s} does not contain the AtomicCoordinatesAndAtomicSpecies block, basis set definition may not contain all atoms."
+            f"{self!r} does not contain the AtomicCoordinatesAndAtomicSpecies block, basis set definition may not contain all atoms."
         )
         return Atoms(atoms)
 
@@ -2077,7 +2077,7 @@ def _r_add_overlap(self, parent_call, M):
                 raise ValueError
         except Exception:
             warn(
-                f"{self!s} could not succesfully read the overlap matrix in {parent_call}."
+                f"{self!r} could not succesfully read the overlap matrix in {parent_call}."
             )
 
     def read_density_matrix(self, *args, **kwargs) -> DensityMatrix:
@@ -2292,13 +2292,14 @@ def _r_hamiltonian_hsx(self, *args, **kwargs):
             if "atoms" not in kwargs:
                 kwargs["atoms"] = self.read_basis(order="^hsx")
             H = hsx.read_hamiltonian(*args, **kwargs)
-            Ef = self.read_fermi_level()
-            if Ef is None:
-                info(
-                    f"{self!s}.read_hamiltonian from HSX file failed shifting to the Fermi-level."
-                )
-            else:
-                H.shift(-Ef)
+            if hsx.version == 0:
+                Ef = self.read_fermi_level()
+                if Ef is None:
+                    info(
+                        f"{self!r}.read_hamiltonian from HSX file failed shifting to the Fermi-level."
+                    )
+                else:
+                    H.shift(-Ef)
         return H
 
     @default_ArgumentParser(description="Manipulate a FDF file.")