diff --git a/files b/files index ac7aa1da50..df3ccb6f2e 160000 --- a/files +++ b/files @@ -1 +1 @@ -Subproject commit ac7aa1da5077f720e795a4b76c0f38fd3fd5a8a3 +Subproject commit df3ccb6f2e843ccafcb24ca5c416d7ac71e41e5a diff --git a/sisl/io/siesta/_src/wfsx_read.f90 b/sisl/io/siesta/_src/wfsx_read.f90 index 9b4c54d1e9..aec2f940a4 100644 --- a/sisl/io/siesta/_src/wfsx_read.f90 +++ b/sisl/io/siesta/_src/wfsx_read.f90 @@ -1,28 +1,31 @@ ! 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_wfsx_sizes(fname, nspin, no_u, nk, Gamma) - use io_m, only: open_file, close_file + +! All routines read_wfsx_next_* and skip_wfsx_* expect that the file unit is handled externally. +! The file unit should already be open and in the correct position. +! On the other hand, routines with another name manage the opening and closing +! themselves. + +! -------------------------------------------------------------- +! Routines to read header information +! -------------------------------------------------------------- +! These routines read the information written at the beggining +! of the WFSX file. + +subroutine read_wfsx_next_sizes(iu, skip_basis, nspin, no_u, nk, Gamma) use io_m, only: iostat_update implicit none ! Input parameters - character(len=*), intent(in) :: fname + integer, intent(in) :: iu + logical, intent(in) :: skip_basis integer, intent(out) :: nspin, no_u, nk logical, intent(out) :: Gamma -! Define f2py intents -!f2py intent(in) :: fname -!f2py intent(out) :: nspin -!f2py intent(out) :: no_u -!f2py intent(out) :: nk -!f2py intent(out) :: Gamma - -! Internal variables and arrays - integer :: iu, ierr - - call open_file(fname, 'read', 'old', 'unformatted', iu) + ! Internal variables and arrays + integer :: ierr read(iu, iostat=ierr) nk, Gamma call iostat_update(ierr) @@ -33,130 +36,117 @@ subroutine read_wfsx_sizes(fname, nspin, no_u, nk, Gamma) read(iu, iostat=ierr) no_u call iostat_update(ierr) - call close_file(iu) + if (skip_basis) then + read(iu, iostat=ierr) ! This contains basis information + call iostat_update(ierr) + endif -end subroutine read_wfsx_sizes +end subroutine read_wfsx_next_sizes -subroutine read_wfsx_index_info(fname, ispin, ik, k, kw, nwf) +subroutine read_wfsx_sizes(fname, nspin, no_u, nk, Gamma) use io_m, only: open_file, close_file - use io_m, only: iostat_update implicit none - integer, parameter :: dp = selected_real_kind(p=15) - ! Input parameters character(len=*), intent(in) :: fname - integer, intent(in) :: ispin, ik - real(dp), intent(out) :: k(3), kw - integer, intent(out) :: nwf - -! Define f2py intents -!f2py intent(in) :: fname -!f2py intent(in) :: ispin -!f2py intent(in) :: ik -!f2py intent(out) :: k -!f2py intent(out) :: kw -!f2py intent(out) :: nwf + integer, intent(out) :: nspin, no_u, nk + logical, intent(out) :: Gamma - integer :: i -! Internal variables and arrays - integer :: iu, ierr + ! Internal variables and arrays + integer :: iu call open_file(fname, 'read', 'old', 'unformatted', iu) - call skip_wfsx_index(iu, ispin, ik) - - ! read information here - read(iu, iostat=ierr) i, k, kw - call iostat_update(ierr) - read(iu, iostat=ierr) ! ispin - call iostat_update(ierr) - read(iu, iostat=ierr) nwf - call iostat_update(ierr) + call read_wfsx_next_sizes(iu, .false., nspin, no_u, nk, Gamma) call close_file(iu) -end subroutine read_wfsx_index_info - -subroutine skip_wfsx_index(iu, ispin, ik) - use io_m, only: open_file, close_file - use io_m, only: iostat_update +end subroutine read_wfsx_sizes - implicit none +subroutine read_wfsx_next_basis(iu, no_u, atom_indices, atom_labels, & + orb_index_atom, orb_n, orb_simmetry) + use io_m, only: iostat_update - integer, parameter :: sp = selected_real_kind(p=6) - integer, parameter :: dp = selected_real_kind(p=15) + implicit none - ! Input parameters - integer, intent(in) :: iu - integer, intent(in) :: ispin, ik + ! Input parameters + integer, intent(in) :: iu, no_u + integer, intent(out) :: atom_indices(no_u), orb_index_atom(no_u) + character, intent(out) :: atom_labels(no_u, 20), orb_simmetry(no_u, 20) + integer, intent(out) :: orb_n(no_u) -! Define f2py intents -!f2py intent(in) :: iu -!f2py intent(in) :: ispin, ik + ! Internal variables and arrays + integer :: j, ierr - integer :: no_u - integer :: nk, spin, i, j -! Internal variables and arrays - integer :: ierr + read(iu, iostat=ierr) (atom_indices(j), atom_labels(j, :), orb_index_atom(j), & + orb_n(j), orb_simmetry(j, :), j=1,no_u) - read(iu, iostat=ierr) nk !, is_gamma - call iostat_update(ierr) - if ( ik > nk ) then - call iostat_update(-1) - end if + call iostat_update(ierr) - read(iu, iostat=ierr) spin - call iostat_update(ierr) - if ( ispin > spin ) then - call iostat_update(-2) - end if +end subroutine read_wfsx_next_basis - read(iu, iostat=ierr) no_u - call iostat_update(ierr) +! -------------------------------------------------------------- +! Routines to read k point information +! -------------------------------------------------------------- +! These routines read the information for a given k point in the +! WFSX file. - read(iu, iostat=ierr) ! basis-information - call iostat_update(ierr) +subroutine read_wfsx_next_info(iu, ispin, ik, k, kw, nwf) + use io_m, only: iostat_update - ! Read pass the spin - do j = 1 , ispin - 1 - do i = 1 , nk - call skip_k() - end do - end do + implicit none - ! skip the indices not requested - do i = 1 , ik - 1 - call skip_k() - end do + integer, parameter :: dp = selected_real_kind(p=15) -contains + ! Input parameters + integer, intent(in) :: iu + integer, intent(out) :: ispin, ik + real(dp), intent(out) :: k(3), kw + integer, intent(out) :: nwf - subroutine skip_k() - integer :: nwf, iwf + ! Internal variables and arrays + integer :: ierr - read(iu, iostat=ierr) ! ik, k, kw + ! read information here + read(iu, iostat=ierr) ik, k, kw call iostat_update(ierr) - read(iu, iostat=ierr) ! ispin + read(iu, iostat=ierr) ispin call iostat_update(ierr) read(iu, iostat=ierr) nwf call iostat_update(ierr) - do iwf = 1, nwf - read(iu, iostat=ierr) ! indwf - call iostat_update(ierr) - read(iu, iostat=ierr) ! eig [eV] - call iostat_update(ierr) - read(iu, iostat=ierr) ! state - call iostat_update(ierr) - end do +end subroutine read_wfsx_next_info - end subroutine skip_k +subroutine read_wfsx_index_info(fname, ispin, ik, k, kw, nwf) + use io_m, only: open_file, close_file + use io_m, only: iostat_update -end subroutine skip_wfsx_index + implicit none + + integer, parameter :: dp = selected_real_kind(p=15) + + ! Input parameters + character(len=*), intent(in) :: fname + integer, intent(in) :: ispin, ik + real(dp), intent(out) :: k(3), kw + integer, intent(out) :: nwf + + integer :: file_ispin, file_ik + ! Internal variables and arrays + integer :: iu + + call open_file(fname, 'read', 'old', 'unformatted', iu) + + call skip_wfsx_index(iu, ispin, ik) -subroutine read_wfsx_index_check(iu, ispin, ik, nwf, fail) + call read_wfsx_next_info(iu, file_ispin, file_ik, k, kw, nwf) + + call close_file(iu) + +end subroutine read_wfsx_index_info + +subroutine read_wfsx_next_index_check(iu, ispin, ik, nwf, fail) use io_m, only: iostat_update, iostat_query implicit none @@ -186,7 +176,184 @@ subroutine read_wfsx_index_check(iu, ispin, ik, nwf, fail) fail = -3 end if -end subroutine read_wfsx_index_check +end subroutine read_wfsx_next_index_check + +! -------------------------------------------------------------- +! Routines that skip over things +! -------------------------------------------------------------- +! Useful if we want to ignore contents of the file. This routines +! ALWAYS depend on the file unit being handled externally. + +subroutine skip_wfsx_next_eigenstate(iu) + use io_m, only: iostat_update + + integer, intent(in) :: iu + + integer :: nwf + + read(iu, iostat=ierr) ! ik, k, kw + call iostat_update(ierr) + read(iu, iostat=ierr) ! ispin + call iostat_update(ierr) + read(iu, iostat=ierr) nwf + call iostat_update(ierr) + + call skip_wfsx_next_vals(iu, nwf) + +end subroutine skip_wfsx_next_eigenstate + +subroutine skip_wfsx_next_vals(iu, nwf) + use io_m, only: iostat_update + + integer, intent(in) :: iu, nwf + + integer :: iwf + + do iwf = 1, nwf + read(iu, iostat=ierr) ! indwf + call iostat_update(ierr) + read(iu, iostat=ierr) ! eig [eV] + call iostat_update(ierr) + read(iu, iostat=ierr) ! state + call iostat_update(ierr) + end do + +end subroutine skip_wfsx_next_vals + +subroutine skip_wfsx_index(iu, ispin, ik) + use io_m, only: iostat_update + + implicit none + + ! Input parameters + integer, intent(in) :: iu + integer, intent(in) :: ispin, ik + + ! Internal variables + integer :: no_u + integer :: nk, nspin, i, j + logical :: Gamma + + call read_wfsx_next_sizes(iu, .true., nspin, no_u, nk, Gamma) + + ! Check that the requested indices make sense + if ( ik > nk ) call iostat_update(-1) + if ( ispin > nspin ) call iostat_update(-2) + + ! Skip until the k index that we need + do i = 1 , ik - 1 + do j = 1 , nspin + call skip_wfsx_next_eigenstate(iu) + end do + end do + + ! Skip until the desired spin index + do i = 1 , ispin - 1 + call skip_wfsx_next_eigenstate(iu) + end do + +end subroutine skip_wfsx_index + +! -------------------------------------------------------------- +! Routines that read the next values +! -------------------------------------------------------------- +! There are three different cases: +! - 1: When only the gamma point is in the WFSX file. The state +! is in a real array. +! - 2: All other cases where spin is not non-colinear. The state +! is in a complex array. +! - 4: For non-colinear spins. The state's first dimension is two +! times larger than the number of actual orbitals in the system. + +subroutine read_wfsx_next_1(iu, no_u, nwf, idx, eig, state) + use io_m, only: iostat_update + + implicit none + + integer, parameter :: sp = selected_real_kind(p=6) + integer, parameter :: dp = selected_real_kind(p=15) + + ! Input parameters + integer, intent(in) :: iu + integer, intent(in) :: no_u, nwf + integer, intent(out) :: idx(nwf) + real(dp), intent(out) :: eig(nwf) + real(sp), intent(out) :: state(no_u, nwf) + + integer :: iwf + ! Internal variables and arrays + integer :: ierr + + do iwf = 1, nwf + read(iu, iostat=ierr) idx(iwf) + call iostat_update(ierr) + read(iu, iostat=ierr) eig(iwf) + call iostat_update(ierr) + read(iu, iostat=ierr) state(:,iwf) + call iostat_update(ierr) + end do + +end subroutine read_wfsx_next_1 + +subroutine read_wfsx_next_2(iu, no_u, nwf, idx, eig, state) + use io_m, only: iostat_update + + implicit none + + integer, parameter :: sp = selected_real_kind(p=6) + integer, parameter :: dp = selected_real_kind(p=15) + + ! Input parameters + integer, intent(in) :: iu + integer, intent(in) :: no_u, nwf + integer, intent(out) :: idx(nwf) + real(dp), intent(out) :: eig(nwf) + complex(sp), intent(out) :: state(no_u, nwf) + + integer :: iwf + ! Internal variables and arrays + integer :: ierr + + do iwf = 1, nwf + read(iu, iostat=ierr) idx(iwf) + call iostat_update(ierr) + read(iu, iostat=ierr) eig(iwf) + call iostat_update(ierr) + read(iu, iostat=ierr) state(:,iwf) + call iostat_update(ierr) + end do + +end subroutine read_wfsx_next_2 + +subroutine read_wfsx_next_4(iu, no_u, nwf, idx, eig, state) + use io_m, only: iostat_update + + implicit none + + integer, parameter :: sp = selected_real_kind(p=6) + integer, parameter :: dp = selected_real_kind(p=15) + + ! Input parameters + integer, intent(in) :: iu + integer, intent(in) :: no_u, nwf + integer, intent(out) :: idx(nwf) + real(dp), intent(out) :: eig(nwf) + complex(sp), intent(out) :: state(2*no_u, nwf) + + integer :: iwf + ! Internal variables and arrays + integer :: ierr + + do iwf = 1, nwf + read(iu, iostat=ierr) idx(iwf) + call iostat_update(ierr) + read(iu, iostat=ierr) eig(iwf) + call iostat_update(ierr) + read(iu, iostat=ierr) state(:,iwf) + call iostat_update(ierr) + end do + +end subroutine read_wfsx_next_4 subroutine read_wfsx_index_1(fname, ispin, ik, no_u, nwf, idx, eig, state) use io_m, only: open_file, close_file @@ -204,34 +371,21 @@ subroutine read_wfsx_index_1(fname, ispin, ik, no_u, nwf, idx, eig, state) real(dp), intent(out) :: eig(nwf) real(sp), intent(out) :: state(no_u, nwf) -! Define f2py intents -!f2py intent(in) :: fname -!f2py intent(in) :: ispin, ik, no_u, nwf -!f2py intent(out) :: idx, eig, state - - integer :: iwf -! Internal variables and arrays + ! Internal variables and arrays integer :: iu, ierr call open_file(fname, 'read', 'old', 'unformatted', iu) call skip_wfsx_index(iu, ispin, ik) - call read_wfsx_index_check(iu, ispin, ik, nwf, ierr) + call read_wfsx_next_index_check(iu, ispin, ik, nwf, ierr) if ( ierr /= 0 ) then call iostat_update(ierr) call close_file(iu) return end if - do iwf = 1, nwf - read(iu, iostat=ierr) idx(iwf) - call iostat_update(ierr) - read(iu, iostat=ierr) eig(iwf) - call iostat_update(ierr) - read(iu, iostat=ierr) state(:,iwf) - call iostat_update(ierr) - end do + call read_wfsx_next_1(iu, no_u, nwf, idx, eig, state) call close_file(iu) @@ -253,34 +407,21 @@ subroutine read_wfsx_index_2(fname, ispin, ik, no_u, nwf, idx, eig, state) real(dp), intent(out) :: eig(nwf) complex(sp), intent(out) :: state(no_u, nwf) -! Define f2py intents -!f2py intent(in) :: fname -!f2py intent(in) :: ispin, ik, no_u, nwf -!f2py intent(out) :: idx, eig, state - - integer :: iwf -! Internal variables and arrays + ! Internal variables and arrays integer :: iu, ierr call open_file(fname, 'read', 'old', 'unformatted', iu) call skip_wfsx_index(iu, ispin, ik) - call read_wfsx_index_check(iu, ispin, ik, nwf, ierr) + call read_wfsx_next_index_check(iu, ispin, ik, nwf, ierr) if ( ierr /= 0 ) then call iostat_update(ierr) call close_file(iu) return end if - do iwf = 1, nwf - read(iu, iostat=ierr) idx(iwf) - call iostat_update(ierr) - read(iu, iostat=ierr) eig(iwf) - call iostat_update(ierr) - read(iu, iostat=ierr) state(:,iwf) - call iostat_update(ierr) - end do + call read_wfsx_next_2(iu, no_u, nwf, idx, eig, state) call close_file(iu) @@ -302,35 +443,62 @@ subroutine read_wfsx_index_4(fname, ispin, ik, no_u, nwf, idx, eig, state) real(dp), intent(out) :: eig(nwf) complex(sp), intent(out) :: state(2*no_u, nwf) -! Define f2py intents -!f2py intent(in) :: fname -!f2py intent(in) :: ispin, ik, no_u, nwf -!f2py intent(out) :: idx, eig, state - - integer :: iwf -! Internal variables and arrays + ! Internal variables and arrays integer :: iu, ierr call open_file(fname, 'read', 'old', 'unformatted', iu) call skip_wfsx_index(iu, ispin, ik) - call read_wfsx_index_check(iu, ispin, ik, nwf, ierr) + call read_wfsx_next_index_check(iu, ispin, ik, nwf, ierr) if ( ierr /= 0 ) then call iostat_update(ierr) call close_file(iu) return end if - do iwf = 1, nwf - read(iu, iostat=ierr) idx(iwf) - call iostat_update(ierr) - read(iu, iostat=ierr) eig(iwf) - call iostat_update(ierr) - read(iu, iostat=ierr) state(:,iwf) - call iostat_update(ierr) - end do + call read_wfsx_next_4(iu, no_u, nwf, idx, eig, state) call close_file(iu) end subroutine read_wfsx_index_4 + +! -------------------------------------------------------------- +! Routines that read the information of all k points +! -------------------------------------------------------------- +! They skip the actual values of the eigenstates and might be +! useful to check the contents of the file. + +subroutine read_wfsx_next_all_info(iu, nspin, nk, ks, kw, nwf) + ! This function should be called after reading sizes. + ! It reads the k value, weight and number of wavefunctions + ! of all the k points in the file. + use io_m, only: iostat_update + + implicit none + + integer, parameter :: dp = selected_real_kind(p=15) + + ! Input parameters + integer, intent(in) :: iu + integer, intent(in) :: nspin, nk + real(dp), intent(out) :: ks(nk, 3), kw(nk) + integer, intent(out) :: nwf(nspin, nk) + + ! Internal variables and arrays + integer :: ik, ispin, file_ik, file_ispin + integer :: l_nspin + + l_nspin = size(nwf, 1) + + do ik = 1 , nk + do ispin = 1, l_nspin + ! Notice that if there is more than one spin, we write more than once + ! to the same position. It doesn't matter, since the k value and the k weight + ! is the same for all spin indices. + call read_wfsx_next_info(iu, file_ispin, file_ik, ks(ik, :), kw(ik), nwf(ispin, ik)) + call skip_wfsx_next_vals(iu, nwf(ispin, ik)) + end do + end do + +end subroutine read_wfsx_next_all_info diff --git a/sisl/io/siesta/binaries.py b/sisl/io/siesta/binaries.py index ccdf95b97e..7d9d233018 100644 --- a/sisl/io/siesta/binaries.py +++ b/sisl/io/siesta/binaries.py @@ -2,8 +2,8 @@ # 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 numbers import Integral -from itertools import product, groupby -from collections import deque +from itertools import product +from collections import deque, namedtuple import numpy as np from numpy import pi @@ -14,10 +14,10 @@ print(e) found_module = False -from ..sile import add_sile, SileError +from ..sile import add_sile, SileError, SileWarning from .sile import SileBinSiesta from sisl._internal import set_module -from sisl.messages import warn, SislError +from sisl.messages import info, warn, SislError from ._help import * import sisl._array as _a @@ -27,6 +27,7 @@ from sisl.unit.siesta import unit_convert from sisl.physics.sparse import SparseOrbitalBZ from sisl.physics import Hamiltonian, DensityMatrix, EnergyDensityMatrix +from sisl.physics import BrillouinZone from sisl.physics.overlap import Overlap from sisl.physics.electron import EigenstateElectron @@ -43,12 +44,6 @@ _eV2Ry = unit_convert('eV', 'Ry') -def _bin_check(obj, method, message): - ierr = _siesta.io_m.iostat_query() - if ierr != 0: - raise SileError(f'{obj!s}.{method} {message} (ierr={ierr})') - - def _toF(array, dtype, scale=None): if scale is None: return array.astype(dtype, order='F', copy=False) @@ -142,9 +137,9 @@ class onlysSileSiesta(SileBinSiesta): def read_supercell(self): """ Returns a SuperCell object from a TranSiesta file """ n_s = _siesta.read_tshs_sizes(self.file)[3] - _bin_check(self, 'read_supercell', 'could not read sizes.') + self._fortran_check('read_supercell', 'could not read sizes.') arr = _siesta.read_tshs_cell(self.file, n_s) - _bin_check(self, 'read_supercell', 'could not read cell.') + self._fortran_check('read_supercell', 'could not read cell.') nsc = np.array(arr[0], np.int32) # We have to transpose since the data is read *as-is* # The cell in fortran files are (:, A1) @@ -162,9 +157,9 @@ def read_geometry(self, geometry=None): sc = self.read_supercell() na = _siesta.read_tshs_sizes(self.file)[1] - _bin_check(self, 'read_geometry', 'could not read sizes.') + self._fortran_check('read_geometry', 'could not read sizes.') arr = _siesta.read_tshs_geom(self.file, na) - _bin_check(self, 'read_geometry', 'could not read geometry.') + self._fortran_check('read_geometry', 'could not read geometry.') # see onlysSileSiesta.read_supercell for .T xyz = arr[0].T * _Bohr2Ang lasto = arr[1] @@ -216,14 +211,14 @@ def read_overlap(self, **kwargs): # read the sizes used... sizes = _siesta.read_tshs_sizes(self.file) - _bin_check(self, 'read_overlap', 'could not read sizes.') + self._fortran_check('read_overlap', 'could not read sizes.') # see onlysSileSiesta.read_supercell for .T isc = _siesta.read_tshs_cell(self.file, sizes[3])[2].T - _bin_check(self, 'read_overlap', 'could not read cell.') + self._fortran_check('read_overlap', 'could not read cell.') no = sizes[2] nnz = sizes[4] ncol, col, dS = _siesta.read_tshs_s(self.file, no, nnz) - _bin_check(self, 'read_overlap', 'could not read overlap matrix.') + self._fortran_check('read_overlap', 'could not read overlap matrix.') # Create the Hamiltonian container S = Overlap(geom, nnzpr=1) @@ -255,7 +250,7 @@ def read_fermi_level(self): Ef : fermi-level of the system """ Ef = _siesta.read_tshs_ef(self.file) * _Ry2eV - _bin_check(self, 'read_fermi_level', 'could not read fermi-level.') + self._fortran_check('read_fermi_level', 'could not read fermi-level.') return Ef @@ -270,15 +265,15 @@ def read_hamiltonian(self, **kwargs): # read the sizes used... sizes = _siesta.read_tshs_sizes(self.file) - _bin_check(self, 'read_hamiltonian', 'could not read sizes.') + self._fortran_check('read_hamiltonian', 'could not read sizes.') # see onlysSileSiesta.read_supercell for .T isc = _siesta.read_tshs_cell(self.file, sizes[3])[2].T - _bin_check(self, 'read_hamiltonian', 'could not read cell.') + self._fortran_check('read_hamiltonian', 'could not read cell.') spin = sizes[0] no = sizes[2] nnz = sizes[4] ncol, col, dH, dS = _siesta.read_tshs_hs(self.file, spin, no, nnz) - _bin_check(self, 'read_hamiltonian', 'could not read Hamiltonian and overlap matrix.') + self._fortran_check('read_hamiltonian', 'could not read Hamiltonian and overlap matrix.') # Check whether it is an orthogonal basis set orthogonal = np.abs(dS).sum() == geom.no @@ -362,7 +357,7 @@ def write_hamiltonian(self, H, **kwargs): csr.ncol, csr.col + 1, _toF(h, np.float64, _eV2Ry), _toF(s, np.float64), isc) - _bin_check(self, 'write_hamiltonian', 'could not write Hamiltonian and overlap matrix.') + self._fortran_check('write_hamiltonian', 'could not write Hamiltonian and overlap matrix.') @set_module("sisl.io.siesta") @@ -374,10 +369,10 @@ def read_density_matrix(self, **kwargs): # Now read the sizes used... spin, no, nsc, nnz = _siesta.read_dm_sizes(self.file) - _bin_check(self, 'read_density_matrix', 'could not read density matrix sizes.') + self._fortran_check('read_density_matrix', 'could not read density matrix sizes.') ncol, col, dDM = _siesta.read_dm(self.file, spin, no, nsc, nnz) - _bin_check(self, 'read_density_matrix', 'could not read density matrix.') + self._fortran_check('read_density_matrix', 'could not read density matrix.') # Try and immediately attach a geometry geom = kwargs.get('geometry', kwargs.get('geom', None)) @@ -448,7 +443,7 @@ def write_density_matrix(self, DM, **kwargs): nsc = DM.geometry.sc.nsc.astype(np.int32) _siesta.write_dm(self.file, nsc, csr.ncol, csr.col + 1, _toF(dm, np.float64)) - _bin_check(self, 'write_density_matrix', 'could not write density matrix.') + self._fortran_check('write_density_matrix', 'could not write density matrix.') @set_module("sisl.io.siesta") @@ -460,9 +455,9 @@ def read_energy_density_matrix(self, **kwargs): # Now read the sizes used... spin, no, nsc, nnz = _siesta.read_tsde_sizes(self.file) - _bin_check(self, 'read_energy_density_matrix', 'could not read energy density matrix sizes.') + self._fortran_check('read_energy_density_matrix', 'could not read energy density matrix sizes.') ncol, col, dEDM = _siesta.read_tsde_edm(self.file, spin, no, nsc, nnz) - _bin_check(self, 'read_energy_density_matrix', 'could not read energy density matrix.') + self._fortran_check('read_energy_density_matrix', 'could not read energy density matrix.') # Try and immediately attach a geometry geom = kwargs.get('geometry', kwargs.get('geom', None)) @@ -515,7 +510,7 @@ def read_fermi_level(self): Ef : fermi-level of the system """ Ef = _siesta.read_tsde_ef(self.file) * _Ry2eV - _bin_check(self, 'read_fermi_level', 'could not read fermi-level.') + self._fortran_check('read_fermi_level', 'could not read fermi-level.') return Ef def write_density_matrices(self, DM, EDM, Ef=0., **kwargs): @@ -567,7 +562,7 @@ def write_density_matrices(self, DM, EDM, Ef=0., **kwargs): _siesta.write_tsde_dm_edm(self.file, nsc, DMcsr.ncol, DMcsr.col + 1, _toF(dm, np.float64), _toF(edm, np.float64, _eV2Ry), Ef * _eV2Ry) - _bin_check(self, 'write_density_matrices', 'could not write DM + EDM matrices.') + self._fortran_check('write_density_matrices', 'could not write DM + EDM matrices.') @set_module("sisl.io.siesta") @@ -875,7 +870,7 @@ def _read_atoms(self, **kwargs): """ Reads basis set and geometry information from the HSX file """ # Now read the sizes used... no, na, nspecies = _siesta.read_hsx_specie_sizes(self.file) - _bin_check(self, 'read_geometry', 'could not read specie sizes.') + self._fortran_check('read_geometry', 'could not read specie sizes.') # Read specie information labels, val_q, norbs, isa = _siesta.read_hsx_species(self.file, nspecies, no, na) # convert to proper string @@ -884,7 +879,7 @@ def _read_atoms(self, **kwargs): labels = list(map(lambda s: b''.join(s).decode('utf-8').strip(), labels.tolist()) ) - _bin_check(self, 'read_geometry', 'could not read species.') + self._fortran_check('read_geometry', 'could not read species.') # to python index isa -= 1 @@ -923,7 +918,7 @@ def _read_atoms(self, **kwargs): atoms = [] for ispecie in range(nspecies): n_l_zeta = _siesta.read_hsx_specie(self.file, ispecie+1, norbs[ispecie]) - _bin_check(self, 'read_geometry', f'could not read specie {ispecie}.') + self._fortran_check('read_geometry', f'could not read specie {ispecie}.') # create orbital # no shell will have l>5, so m=10 should be more than enough m = 10 @@ -940,11 +935,11 @@ def _read_atoms(self, **kwargs): # now read in xij to retrieve atomic positions Gamma, spin, no, no_s, nnz = _siesta.read_hsx_sizes(self.file) - _bin_check(self, 'read_geometry', 'could not read matrix sizes.') + self._fortran_check('read_geometry', 'could not read matrix sizes.') ncol, col, _, dxij = _siesta.read_hsx_sx(self.file, Gamma, spin, no, no_s, nnz) dxij = dxij.T * _Bohr2Ang col -= 1 - _bin_check(self, 'read_geometry', 'could not read xij matrix.') + self._fortran_check('read_geometry', 'could not read xij matrix.') # now create atoms object atoms = Atoms([atoms[ia] for ia in isa]) @@ -956,11 +951,11 @@ def read_hamiltonian(self, **kwargs): # Now read the sizes used... Gamma, spin, no, no_s, nnz = _siesta.read_hsx_sizes(self.file) - _bin_check(self, 'read_hamiltonian', 'could not read Hamiltonian sizes.') + self._fortran_check('read_hamiltonian', 'could not read Hamiltonian sizes.') ncol, col, dH, dS, dxij = _siesta.read_hsx_hsx(self.file, Gamma, spin, no, no_s, nnz) dxij = dxij.T * _Bohr2Ang col -= 1 - _bin_check(self, 'read_hamiltonian', 'could not read Hamiltonian.') + self._fortran_check('read_hamiltonian', 'could not read Hamiltonian.') ptr = _ncol_to_indptr(ncol) xij = SparseCSR((dxij, col, ptr), shape=(no, no_s)) @@ -997,11 +992,11 @@ def read_overlap(self, **kwargs): """ Returns the overlap matrix from the siesta.HSX file """ # Now read the sizes used... Gamma, spin, no, no_s, nnz = _siesta.read_hsx_sizes(self.file) - _bin_check(self, 'read_overlap', 'could not read overlap matrix sizes.') + self._fortran_check('read_overlap', 'could not read overlap matrix sizes.') ncol, col, dS, dxij = _siesta.read_hsx_sx(self.file, Gamma, spin, no, no_s, nnz) dxij = dxij.T * _Bohr2Ang col -= 1 - _bin_check(self, 'read_overlap', 'could not read overlap matrix.') + self._fortran_check('read_overlap', 'could not read overlap matrix.') ptr = _ncol_to_indptr(ncol) xij = SparseCSR((dxij, col, ptr), shape=(no, no_s)) @@ -1035,56 +1030,474 @@ def read_overlap(self, **kwargs): @set_module("sisl.io.siesta") class wfsxSileSiesta(SileBinSiesta): - r""" Binary WFSX file reader for Siesta """ + r""" Binary WFSX file reader for Siesta - def yield_eigenstate(self, parent=None): - r""" Reads eigenstates from the WFSX file + The WFSX file assumes that users initialize the object with + a `parent` argument (or one of the other geometry related objects as + shown below). - Returns - ------- - state: EigenstateElectron + The `parent` argument is necessary to convert WFSX k-points from 1/Ang to + reduced coordinates. + When returning `EigenstateElectron` objects the parent of these objects + are the equivalent of the `parent` argument upon initialization. + Therefore please remember to pass a correct `parent`. + + 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] + sc : SuperCell, optional + a supercell contains the lattice vectors to convert k + """ + + def _setup(self, *args, **kwargs): + """ Simple setup that needs to be overwritten + + All _read_next_* methods expect the fortran file unit to be handled + and that the position in the file is correct. """ - # First query information - nspin, nou, nk, Gamma = _siesta.read_wfsx_sizes(self.file) - _bin_check(self, 'yield_eigenstate', 'could not read sizes.') - - if nspin in [4, 8]: - nspin = 1 # only 1 spin - func = _siesta.read_wfsx_index_4 - elif Gamma: - func = _siesta.read_wfsx_index_1 + super()._setup(*args, **kwargs) + + # default sc + sc = 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, SuperCell): + sc = parent else: - func = _siesta.read_wfsx_index_2 + geometry = parent.geometry + + geometry = kwargs.get("geometry", geometry) + if geometry is not None and sc is None: + sc = geometry.sc + sc = kwargs.get("sc", sc) + if sc is None and geometry is not None: + raise ValueError(f"{self.__class__.__name__}(geometry=Geometry, sc=None) is not an allowed argument combination.") + + if parent is None: + parent = geometry if parent is None: - def convert_k(k): + parent = sc + + self._parent = parent + self._geometry = geometry + self._sc = sc + if self._parent is None and self._geometry is None and self._sc is None: + def conv(k): if not np.allclose(k, 0.): - warn(f"{self.__class__.__name__}.yield_eigenstate returns a k-point in 1/Ang (not in reduced format), please pass 'parent' to ensure reduced k") - return k + warn(f"{self.__class__.__name__} cannot convert stored k-points from 1/Ang to reduced coordinates. Please ensure 'parent=Hamiltonian', 'geometry=Geometry', or 'sc=SuperCell' to ensure reduced k") + return k / _Bohr2Ang else: - # We can succesfully convert to proper reduced k-points - if isinstance(parent, SuperCell): - def convert_k(k): - return np.dot(k, parent.cell.T) / (2 * pi) - else: - def convert_k(k): - return np.dot(k, parent.sc.cell.T) / (2 * pi) + def conv(k): + return (k @ sc.cell.T) / (2 * pi * _Bohr2Ang) + self._convert_k = conv - for ispin, ik in product(range(1, nspin + 1), range(1, nk + 1)): - k, _, nwf = _siesta.read_wfsx_index_info(self.file, ispin, ik) - # Convert to 1/Ang - k /= _Bohr2Ang - _bin_check(self, 'yield_eigenstate', f"could not read index info [{ispin}, {ik}]") + def _open_wfsx(self, mode, rewind=False): + """Open the file unit for the WFSX file. - idx, eig, state = func(self.file, ispin, ik, nou, nwf) - _bin_check(self, 'yield_eigenstate', f"could not read state information [{ispin}, {ik}, {nwf}]") + Here we also initialize some variables to keep track of the state of the read. + """ + self._fortran_open(mode, rewind=rewind) + # 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 + # _ispin is the current (Python) index for the spin-index to be read (only has meaning for a spin-polarized + # WFSX files) + # _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 basis information + # 1 : it means that the file-descriptor is in front of k point information + # 2 : it means that the file-descriptor is in front of k point WFSX values + # + self._state = -1 + self._ispin = 0 + self._ik = 0 - # eig is already in eV - # we probably need to add spin - # see onlysSileSiesta.read_supercell for .T - es = EigenstateElectron(state.T, eig, parent=parent, - k=convert_k(k), gauge="r", index=idx - 1) - yield es + def _close_wfsx(self): + """Close the file unit for the WFSX file. + + We clean the variables used to keep track of read state. + """ + self._fortran_close() + + # Clean variables + del self._state + del self._ik + del self._ispin + try: + del self._sizes + except: + pass + try: + del self._basis + except: + pass + try: + del self._funcs + except: + pass + + def _setup_parsing(self, close=True, skip_basis=True): + """Gets all the things needed to parse the wfsx file. + + Parameters + ----------- + close: bool, optional + Whether the file unit should be closed afterwards. + skip_basis : bool, optional + whether to also read the basis or not + """ + self._open_wfsx('r') + # Read the sizes relevant to the file. + # We also read whether there's only gamma point information or there are multiple points + self._sizes = self._read_next_sizes(skip_basis=skip_basis) + if not skip_basis: + self._basis = self._read_next_basis() + + # Get the functions that should be used to parse state values. + if self._sizes.nspin in (4, 8): + # We will have twice as many coefficients. + func_index = 4 + elif self._sizes.Gamma: + # State values will be in double precision floats + func_index = 1 + else: + # State values will be in double precision complex + func_index = 2 + + Funcs = namedtuple("WFSXReads", ["read_index", "read_next"]) + self._funcs = Funcs( + getattr(_siesta, f"read_wfsx_index_{func_index}"), + getattr(_siesta, f"read_wfsx_next_{func_index}") + ) + + if close: + self._close_wfsx() + + def _read_next_sizes(self, skip_basis=False): + """Reads the sizes if they are the next thing to be read. + + Parameters + ----------- + skip_basis: boolean, optional + Whether this method should also skip over the basis information. + + Returns + ------- + namedtuple : + - 'nspin': int. Number of spin components. + - 'no_u': int. Number of orbitals in the unit cell. + - 'nk': int. Number of k points in the file. + - 'Gamma': bool. Whether the file contains only the gamma point. + """ + # Check that we are in the right position in the file + if self._state != -1: + raise SileError(f"We are not in a position to read the sizes. State is: {self._state}") + # Read the sizes that we can find in the WFSX file + Sizes = namedtuple("Sizes", ["nspin", "no_u", "nk", "Gamma"]) + sizes = _siesta.read_wfsx_next_sizes(self._iu, skip_basis) + # Inform that we are now in front of k point information + self._state = 1 if skip_basis else 0 + # Check that everything went fine + self._fortran_check('_read_sizes', "could not read sizes") + + return Sizes(*sizes) + + def _read_next_basis(self): + """Reads the basis if it is the next thing to be read. + + Returns + ------- + Atoms: + the basis read. + """ + # Check that we are in the right position in the file + if self._state != 0: + raise SileError(f"We are not in a position to read the basis. State is: {self._state}") + # Read the basis information that we can find in the WFSX file + basis_info = _siesta.read_wfsx_next_basis(self._iu, self._sizes.no_u) + # Inform that we are now in front of k point information + self._state = 1 + # Check that everything went fine + self._fortran_check('_read_basis', "could not read basis information") + + # Convert the information to a dict so that code is easier to follow. + basis_info = dict(zip(('atom_indices', 'atom_labels', 'orb_index_atom', 'orb_n', 'orb_symmetry'), basis_info)) + + # Sanitize the string information + for char_key in ("atom_labels", "orb_symmetry"): + basis_info[char_key] = np.array(["".join(label).rstrip() for label in basis_info[char_key].astype(str)]) + + # Find out the unique atom indices + unique_ats = np.unique(basis_info["atom_indices"]) + + def _get_atom_object(at): + """Given an atom index, generates an Atom object with all the information we have about it""" + atom_orbs = np.where(basis_info["atom_indices"] == at)[0] + at_label = basis_info["atom_labels"][atom_orbs[0]] + + orbitals = [ + AtomicOrbital(f"{n}{symmetry}") + for n, symmetry in zip(basis_info["orb_n"][atom_orbs], basis_info["orb_symmetry"][atom_orbs]) + ] + + return Atom(at_label, orbitals=orbitals) + + # Generate the Atoms oject. + return Atoms([_get_atom_object(at) for at in unique_ats]) + + def _read_next_info(self, ispin, ik): + """Reads the next eigenstate information. + + This function should only be called after reading the sizes + or reading the previous state's values. + + Parameters + ---------- + ispin: integer + the (python) spin index of the next eigenstate. + ik: integer + the (python) k index of the next eigenstate. + + Returns + ------- + array of shape (3,): + The k point of the state. + float: + The weight of the k point. + int: + Number of wavefunctions that the state contains. It is needed by the + function that reads the eigenstates values. + """ + # Store the indices of the current state + self._ik = ik + self._ispin = ispin + # Check that we are in a position where we will read state information + if self._state != 1: + raise SileError(f"We are not in a position to read k point information. State is: {self._state}") + + # Read the state information + file_ispin, file_ik, k, weight, nwf = _siesta.read_wfsx_next_info(self._iu) + # Inform that we are now in front of state values + self._state = 2 + # Check that the read went fine + self._fortran_check('_read_next_info', f"could not read next eigenstate info [{ispin + 1}, {ik + 1}]") + + # Check that the read indices match the indices that we were expecting. + if file_ispin != ispin + 1 or file_ik != ik + 1: + self._ik = file_ik - 1 + self._ispin = file_ispin - 1 + raise SileError(f"WFSX indices do not match the expected ones. Expected: [{ispin + 1}, {ik + 1}], found [{file_ispin}, {file_ik}]") + + return k, weight, nwf + + def _read_next_values(self, ispin, ik, nwf): + """Reads the next eigenstate values. + + This function should only be called after reading the state's information + + Parameters + ---------- + ispin: integer + the (python) spin index of the next eigenstate. + ik: integer + the (python) k index of the next eigenstate. + nwf: integer + The number of wavefunctions that the next eigenstate contains. + Should have been obtained by reading the state's info with `_read_next_info`. + + Returns + ------- + array of shape (nwf,): + The indices of each wavefunction that the state contains. + array of shape (nwf,): + The eigenvalues (in eV) of each wavefunction that the state contains. + array of shape (norbitals, nwf): + The coefficients for each wavefunction that the state contains. + """ + # Check that we are in the right position in the file + if self._state != 2: + raise SileError(f"We are not in a position to read k point WFSX values. State is: {self._state}") + + # Read the state values + idx, eig, state = self._funcs.read_next(self._iu, self._sizes.no_u, nwf) + # Inform that we are now in front of the next state info + self._state = 1 + # Check that everything went fine + self._fortran_check('_read_next_values', f"could not read next eigenstate values [{ispin + 1}, {ik + 1}]") + + return idx, eig, state + + def _read_next_eigenstate(self, ispin, ik): + """Reads the next eigenstate in the WFSX file. + + Parameters + ---------- + ispin: integer + the (python) spin index of the next eigenstate. + ik: integer + the (python) k index of the next eigenstate. + + Returns + -------- + EigenstateElectron: + The next eigenstate. + """ + # Get the information of this eigenstate + k, weight, nwf = self._read_next_info(ispin, ik) + # Now that we have the information, we can read the values because + # we know the number of wavefunctions stored in the k point (`nwf`) + idx, eig, state = self._read_next_values(ispin, ik, nwf) + + # Build the info dictionary for the eigenstate to know how it was calculated + # We include the spin index if needed. + info = dict(k=self._convert_k(k), weight=weight, gauge="r", index=idx - 1) + if self._sizes.nspin == 2: + info["spin"] = ispin + + # `eig` is already in eV + # See onlysSileSiesta.read_supercell to understand why we transpose `state` + return EigenstateElectron(state.T, eig, parent=self._parent, **info) + + def read_sizes(self): + """Reads the sizes related to this WFSX file + + Returns + ------- + int : number of spin components + int : number of orbitals in the unit-cell + int : number of k-points + bool : True if the file only contains the Gamma-point + """ + self._open_wfsx("r") + sizes = self._read_next_sizes() + self._close_wfsx() + return sizes + + def read_basis(self): + """Reads the basis contained in the WFSX file. + + The WFSX file only contains information about the atom labels, which atom + each orbital belongs to and the orbital quantum numbers. It is thus not + complete in every sense. + + Returns + ------- + Atoms: + the basis read. + """ + self._open_wfsx("r") + self._sizes = self._read_next_sizes(skip_basis=False) + basis = self._read_next_basis() + self._close_wfsx() + return basis + + def yield_eigenstate(self): + r""" Iterates over the states in the WFSX file + + Yields + ------ + EigenstateElectron + """ + # Open file and get parsing information + self._setup_parsing(close=False) + + if self._sizes.nspin == 2: + itt_spin = range(2) + else: + itt_spin = range(1) + + try: + # Iterate over all eigenstates in the WFSX file, yielding control to the caller at + # each iteration. + for ik, ispin in product(range(self._sizes.nk), itt_spin): + yield self._read_next_eigenstate(ispin, ik) + # We ran out of eigenstates + self._close_wfsx() + except GeneratorExit: + # The loop in which the generator was used has been broken. + self._close_wfsx() + + def read_eigenstate(self, k=(0, 0, 0), spin=0, ktol=1e-4): + """Reads a specific eigenstate from the file. + + This method iterates over the states until it finds a match. Don't 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. + spin: integer, optional + The spin index of the state you want to find. Only meaningful for polarized + calculations. + ktol: float, optional + The threshold value for considering two k-points the same (i.e. to match + the query k point with the state's 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 state.info.get("spin", 0) == spin and np.allclose(state.info['k'], k, atol=ktol): + # This is the state that the user requested + return state + return None + + def read_info(self): + """Reads the information for all the k points contained in the file + + Returns + ------- + k: array of shape (nk, 3) + k values of the k points contained in the file. + weight: array of shape (nk,) + weight of each k point + nwf: array of shape (nspin, nk) + number of wavefunctions that each kpoint(-spin) contains. + """ + # Open file and get parsing information + self._setup_parsing(close=False, skip_basis=True) + + # Check if we are in the correct position in the file (we should be just after the header) + if self._state != 1: + raise ValueError(f"We are not in a position to read eigenstate info in the file. State: {self._state}") + + # Read all the information. Parse here the k values obtained. + # Store the information that should be returned under `returns`. + if self._sizes.nspin == 2: + nspin = 2 + else: + nspin = 1 + k, kw, nwf = _siesta.read_wfsx_next_all_info(self._iu, 1, self._sizes.nk) + + # Close the file unit + self._close_wfsx() + # Check for errors in the read. + self._fortran_check("read_info", "could not read file information.") + return self._convert_k(k), kw, nwf + + def read_brillouinzone(self): + """ Read the brillouin zone object """ + k, weight, _ = self.read_info() + return BrillouinZone(self._parent, k=k, weight=weight) @set_module("sisl.io.siesta") @@ -1099,7 +1512,7 @@ def read_supercell(self, *args, **kwargs): r""" Return the cell contained in the file """ cell = _siesta.read_grid_cell(self.file).T * _Bohr2Ang - _bin_check(self, 'read_supercell', 'could not read cell.') + self._fortran_check('read_supercell', 'could not read cell.') return SuperCell(cell) @@ -1114,7 +1527,7 @@ def read_grid_size(self): # Read the sizes nspin, mesh = _siesta.read_grid_sizes(self.file) - _bin_check(self, 'read_grid_size', 'could not read grid sizes.') + self._fortran_check('read_grid_size', 'could not read grid sizes.') return nspin, mesh def read_grid(self, index=0, dtype=np.float64, *args, **kwargs): @@ -1135,7 +1548,7 @@ def read_grid(self, index=0, dtype=np.float64, *args, **kwargs): nspin, mesh = self.read_grid_size() sc = self.read_supercell() grid = _siesta.read_grid(self.file, nspin, mesh[0], mesh[1], mesh[2]) - _bin_check(self, 'read_grid', 'could not read grid.') + self._fortran_check('read_grid', 'could not read grid.') if isinstance(index, Integral): grid = grid[:, :, :, index] @@ -1184,34 +1597,22 @@ class _gfSileSiesta(SileBinSiesta): ... if new_k: ... f.write_hamiltonian(H, S) ... f.write_self_energy(SeHSE) + """ def _setup(self, *args, **kwargs): """ Simple setup that needs to be overwritten """ - self._iu = -1 + super()._setup(*args, **kwargs) + # The unit convention used for energy-points # This is necessary until Siesta uses CODATA values - if kwargs.get("unit", "old").lower() in ("old", "4.1"): + if kwargs.get("version", "old").lower() in ("old", "4.1"): self._E_Ry2eV = 13.60580 else: self._E_Ry2eV = _Ry2eV - def _is_open(self): - return self._iu != -1 - def _open_gf(self, mode, rewind=False): - if self._is_open() and mode == self._mode: - if rewind: - _siesta.io_m.rewind_file(self._iu) - else: - # retain indices - return - else: - if mode == 'r': - self._iu = _siesta.io_m.open_file_read(self.file) - elif mode == 'w': - self._iu = _siesta.io_m.open_file_write(self.file) - _bin_check(self, '_open_gf', 'could not open for {}.'.format({'r': 'reading', 'w': 'writing'}[mode])) + self._fortran_open(mode, rewind=rewind) # They will at any given time # correspond to the current Python indices that is to be read @@ -1236,11 +1637,9 @@ def _open_gf(self, mode, rewind=False): self._iE = 0 def _close_gf(self): - if not self._is_open(): + if not self._fortran_is_open(): return - # Close it - _siesta.io_m.close_file(self._iu) - self._iu = -1 + self._fortran_close() # Clean variables del self._state @@ -1401,12 +1800,12 @@ def read_header(self): E : energy points in the GF file """ # Ensure it is open (in read-mode) - if self._is_open(): + if self._fortran_is_open(): _siesta.io_m.rewind_file(self._iu) else: self._open_gf('r') nspin, no_u, nkpt, NE = _siesta.read_gf_sizes(self._iu) - _bin_check(self, 'read_header', 'could not read sizes.') + self._fortran_check('read_header', 'could not read sizes.') self._nspin = nspin self._nk = nkpt self._nE = NE @@ -1415,7 +1814,7 @@ def read_header(self): _siesta.io_m.rewind_file(self._iu) self._step_counter('read_header', header=True, read=True) k, E = _siesta.read_gf_header(self._iu, nkpt, NE) - _bin_check(self, 'read_header', 'could not read header information.') + self._fortran_check('read_header', 'could not read header information.') if self._nspin > 2: # non-colinear self._no_u = no_u * 2 @@ -1433,7 +1832,7 @@ def disk_usage(self): ------- estimated disk-space used in GB """ - is_open = self._is_open() + is_open = self._fortran_is_open() if not is_open: self.read_header() @@ -1462,7 +1861,7 @@ def read_hamiltonian(self): """ self._step_counter('read_hamiltonian', HS=True, read=True) H, S = _siesta.read_gf_hs(self._iu, self._no_u) - _bin_check(self, 'read_hamiltonian', 'could not read Hamiltonian and overlap matrices.') + self._fortran_check('read_hamiltonian', 'could not read Hamiltonian and overlap matrices.') # we don't convert to C order! return H * _Ry2eV, S @@ -1480,7 +1879,7 @@ def read_self_energy(self): """ self._step_counter('read_self_energy', read=True) SE = _siesta.read_gf_se(self._iu, self._no_u, self._iE) - _bin_check(self, 'read_self_energy', 'could not read self-energy.') + self._fortran_check('read_self_energy', 'could not read self-energy.') # we don't convert to C order! return SE * _Ry2eV @@ -1496,7 +1895,7 @@ def HkSk(self, k=(0, 0, 0), spin=0): spin : int, optional spin-index for the Hamiltonian and overlap matrices """ - if not self._is_open(): + if not self._fortran_is_open(): self.read_header() # find k-index that is requested @@ -1504,7 +1903,7 @@ def HkSk(self, k=(0, 0, 0), spin=0): _siesta.read_gf_find(self._iu, self._nspin, self._nk, self._nE, self._state, self._ispin, self._ik, self._iE, self._is_read, 0, spin, ik, 0) - _bin_check(self, 'HkSk', 'could not find Hamiltonian and overlap matrix.') + self._fortran_check('HkSk', 'could not find Hamiltonian and overlap matrix.') self._state = 0 self._ispin = spin @@ -1525,7 +1924,7 @@ def self_energy(self, E, k=0, spin=0): spin : int, optional spin-index to retrieve self-energy at """ - if not self._is_open(): + if not self._fortran_is_open(): self.read_header() ik = self.kindex(k) @@ -1533,7 +1932,7 @@ def self_energy(self, E, k=0, spin=0): _siesta.read_gf_find(self._iu, self._nspin, self._nk, self._nE, self._state, self._ispin, self._ik, self._iE, self._is_read, 1, spin, ik, iE) - _bin_check(self, 'self_energy', 'could not find requested self-energy.') + self._fortran_check('self_energy', 'could not find requested self-energy.') self._state = 1 self._ispin = spin @@ -1606,7 +2005,7 @@ def write_header(self, bz, E, mu=0., obj=None): lasto, bloch, 0, mu * _eV2Ry, _toF(k.T, np.float64), w, self._E / self._E_Ry2eV, **sizes) - _bin_check(self, 'write_header', 'could not write header information.') + self._fortran_check('write_header', 'could not write header information.') def write_hamiltonian(self, H, S=None): """ Write the current energy, k-point and H and S to the file @@ -1626,7 +2025,7 @@ def write_hamiltonian(self, H, S=None): _siesta.write_gf_hs(self._iu, self._ik, self._E[self._iE] / self._E_Ry2eV, _toF(H, np.complex128, _eV2Ry), _toF(S, np.complex128), no_u=no) - _bin_check(self, 'write_hamiltonian', 'could not write Hamiltonian and overlap matrices.') + self._fortran_check('write_hamiltonian', 'could not write Hamiltonian and overlap matrices.') def write_self_energy(self, SE): r""" Write the current self energy, k-point and H and S to the file @@ -1645,7 +2044,7 @@ def write_self_energy(self, SE): self._step_counter('write_self_energy', read=True) _siesta.write_gf_se(self._iu, self._ik, self._iE, self._E[self._iE] / self._E_Ry2eV, _toF(SE, np.complex128, _eV2Ry), no_u=no) - _bin_check(self, 'write_self_energy', 'could not write self-energy.') + self._fortran_check('write_self_energy', 'could not write self-energy.') def __len__(self): return self._nE * self._nk * self._nspin @@ -1660,17 +2059,19 @@ def __iter__(self): # get everything e = self._E if self._nspin in [1, 2]: + GFStep = namedtuple("GFStep", ["spin", "do_HS", "k", "E"]) for ispin in range(self._nspin): for k in self._k: - yield ispin, True, k, e[0] + yield GFStep(ispin, True, k, e[0]) for E in e[1:]: - yield ispin, False, k, E + yield GFStep(ispin, False, k, E) else: + GFStep = namedtuple("GFStep", ["do_HS", "k", "E"]) for k in self._k: - yield True, k, e[0] + yield GFStep(True, k, e[0]) for E in e[1:]: - yield False, k, E + yield GFStep(False, k, E) # We will automatically close once we hit the end self._close_gf() diff --git a/sisl/io/siesta/sile.py b/sisl/io/siesta/sile.py index 1ab768b5a2..be8a97d8a6 100644 --- a/sisl/io/siesta/sile.py +++ b/sisl/io/siesta/sile.py @@ -1,8 +1,16 @@ # 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/. + +try: + from . import _siesta + found_bin_module = True +except Exception as e: + print(e) + found_bin_module = False + from sisl._internal import set_module -from ..sile import Sile, SileCDF, SileBin +from ..sile import Sile, SileCDF, SileBin, SileError __all__ = ['SileSiesta', 'SileCDFSiesta', 'SileBinSiesta'] @@ -17,6 +25,8 @@ class SileCDFSiesta(SileCDF): # all netcdf output should not be masked def _setup(self, *args, **kwargs): + super()._setup(*args, **kwargs) + # all NetCDF routines actually returns masked arrays # this is to prevent Siesta CDF files from doing this. if hasattr(self, "fh"): @@ -25,4 +35,39 @@ def _setup(self, *args, **kwargs): @set_module("sisl.io.siesta") class SileBinSiesta(SileBin): - pass + + def _setup(self, *args, **kwargs): + """We set up everything to handle the fortran I/O unit""" + super()._setup(*args, **kwargs) + self._iu = -1 + + def _fortran_check(self, method, message): + ierr = _siesta.io_m.iostat_query() + if ierr != 0: + raise SileError(f'{self!s}.{method} {message} (ierr={ierr})') + + def _fortran_is_open(self): + return self._iu != -1 + + def _fortran_open(self, mode, rewind=False): + if self._fortran_is_open() and mode == self._mode: + if rewind: + _siesta.io_m.rewind_file(self._iu) + else: + # retain indices + return + else: + if mode == 'r': + self._iu = _siesta.io_m.open_file_read(self.file) + elif mode == 'w': + self._iu = _siesta.io_m.open_file_write(self.file) + else: + raise SileError(f"Mode '{mode}' is not an accepted mode to open a fortran file unit. Use only 'r' or 'w'") + self._fortran_check('_fortran_open', 'could not open for {}.'.format({'r': 'reading', 'w': 'writing'}[mode])) + + def _fortran_close(self): + if not self._fortran_is_open(): + return + # Close it + _siesta.io_m.close_file(self._iu) + self._iu = -1 diff --git a/sisl/io/siesta/tests/test_gf.py b/sisl/io/siesta/tests/test_gf.py index 007579952b..964bd3c4c8 100644 --- a/sisl/io/siesta/tests/test_gf.py +++ b/sisl/io/siesta/tests/test_gf.py @@ -90,7 +90,7 @@ def test_gf_write_read_spin(sisl_tmp, sisl_system): gf.write_self_energy(S * e - Hk) # Check it isn't opened - assert not gf._is_open() + assert not gf._fortran_is_open() nspin, no_u, k, E_file = gf.read_header() assert nspin == 2 @@ -132,7 +132,7 @@ def test_gf_write_read_direct(sisl_tmp, sisl_system): gf.write_self_energy(S * e - Hk) # ensure it is not opened - assert not gf._is_open() + assert not gf._fortran_is_open() # First try from beginning for e in [0, 1, E[1], 0, E[0]]: diff --git a/sisl/io/siesta/tests/test_wfsx.py b/sisl/io/siesta/tests/test_wfsx.py new file mode 100644 index 0000000000..e695482546 --- /dev/null +++ b/sisl/io/siesta/tests/test_wfsx.py @@ -0,0 +1,32 @@ +# 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/. +import pytest +import os.path as osp +import sisl +import numpy as np + +pytestmark = [pytest.mark.io, pytest.mark.siesta] +_dir = osp.join('sisl', 'io', 'siesta') + + +def test_wfsx_read(sisl_files): + fdf = sisl.get_sile(sisl_files(_dir, 'bi2se3_3ql.fdf')) + wfsx = sisl.get_sile(sisl_files(_dir, 'bi2se3_3ql.bands.WFSX'), parent=fdf.read_geometry()) + + info = wfsx.read_info() + sizes = wfsx.read_sizes() + basis = wfsx.read_basis() + + # yield states + nstates = 0 + nstates_total = 0 + for state in wfsx.yield_eigenstate(): + nstates += 1 + nstates_total += len(state) + # in this example we a SOC calculation, so we should not muliply by spin + assert nstates == len(info[0]) + assert nstates_total == info[2].sum() + + bz = wfsx.read_brillouinzone() + assert len(bz) == 16 diff --git a/sisl/viz/plots/bands.py b/sisl/viz/plots/bands.py index 672d09a724..3b986666a8 100644 --- a/sisl/viz/plots/bands.py +++ b/sisl/viz/plots/bands.py @@ -6,8 +6,12 @@ import itertools import numpy as np +from numpy.core.defchararray import array +import xarray as xr import sisl +from sisl.physics.brillouinzone import BrillouinZone +from sisl.physics.spin import Spin from ..plot import Plot, entry_point from ..plotutils import find_files from ..input_fields import ( @@ -41,6 +45,11 @@ class BandsPlot(Plot): of the dict: { 'x': 'y': 'z': 'divisions': 'name': Tick that should be displayed at this corner of the path. } + wfsx_file: wfsxSileSiesta, optional + The WFSX file to get the eigenstates. In standard SIESTA + nomenclature, this should probably be the *.bands.WFSX file, as it is + the one that contains the eigenstates for the band + structure. aiida_bands: optional An aiida BandsData node. add_band_data: optional @@ -122,6 +131,15 @@ class BandsPlot(Plot): BandStructureInput(key="band_structure", name="Band structure"), + SileInput(key='wfsx_file', name='Path to WFSX file', + dtype=sisl.io.siesta.wfsxSileSiesta, + default=None, + help="""The WFSX file to get the eigenstates. + In standard SIESTA nomenclature, this should probably be the *.bands.WFSX file, as it is the one + that contains the eigenstates for the band structure. + """ + ), + AiidaNodeInput(key="aiida_bands", name="Aiida BandsData node", default=None, help="""An aiida BandsData node.""" @@ -285,13 +303,25 @@ def _after_init(self): self.add_shortcut("g", "Toggle gap", self.toggle_gap) + @entry_point('bands file', 0) + def _read_siesta_output(self, bands_file, band_structure): + """ + Reads the bands information from a SIESTA bands file. + """ + if band_structure: + raise ValueError("A path was provided, therefore we can not use the .bands file even if there is one") + + self.bands_data = self.get_sile(bands_file or "bands_file").read_data(as_dataarray=True) + + # Define the spin class of the results we have retrieved + if len(self.bands_data.spin.values) == 2: + self.spin = sisl.Spin("p") + @entry_point('aiida bands', 1) def _read_aiida_bands(self, aiida_bands): """ Creates the bands plot reading from an aiida BandsData node. """ - import xarray as xr - plot_data = aiida_bands._get_bandplot_data(cartesian=True) bands = plot_data["y"] @@ -317,12 +347,178 @@ def _read_aiida_bands(self, aiida_bands): attrs={**tick_info} ) - @entry_point('band structure', 2) + def _get_eigenstate_wrapper(self, k_vals, extra_vars=()): + """Helper function to build the function to call on each eigenstate. + + Parameters + ---------- + k_vals: array_like of shape (nk,) + The (linear) values of the k points. This will be used for plotting + the bands. + extra_vars: array-like of dict, optional + This argument determines the extra quantities that should be included + in the final dataset of the bands. Energy and spin moments (if available) + are already included, so no need to pass them here. + Each item of the array defines a new quantity and should contain a dictionary + with the following keys: + - 'name', str: The name of the quantity. + - 'getter', callable: A function that gets 3 arguments: eigenstate, plot and + spin index, and returns the values of the quantity in a numpy array. This + function will be called for each eigenstate object separately. That is, once + for each (k-point, spin) combination. + - 'coords', tuple of str: The names of the dimensions of the returned array. + The number of coordinates should match the number of dimensions. + of + - 'coords_values', dict: If this variable introduces a new coordinate, you should + pass the values for that coordinate here. If the coordinates were already defined + by another variable, they will already have values. If you are unsure that the + coordinates are new, just pass the values for them, they will get overwritten. + + Returns + -------- + function: + The function that should be called for each eigenstate and will return a tuple of size + n_vars with the values for each variable. + tuple of dicts: + A tuple containing the dictionaries that define all variables. Exactly the same as + the passed `extra_vars`, but with the added Energy and spin moment (if available) variables. + dict: + Dictionary containing the values for each coordinate involved in the dataset. + """ + # In case it is a non_colinear or spin-orbit calculation we will get the spin moments + if not self.spin.is_diagonal: + def _spin_moment_getter(eigenstate, plot, spin): + return eigenstate.spin_moment().real + + extra_vars = ({ + "coords": ("band", "axis"), "coords_values": dict(axis=["x", "y", "z"]), + "name": "spin_moments", "getter": _spin_moment_getter}, + *extra_vars) + + # Define the available spin indices. Notice that at the end the spin dimension + # is removed from the dataset unless the calculation is spin polarized. So having + # spin_indices = [0] is just for convenience. + spin_indices = [0] + if self.spin.is_polarized: + spin_indices = [0, 1] + + # Add a variable to get the eigenvalues. + all_vars = ({ + "coords": ("band",), "coords_values": {"spin": spin_indices, "k": k_vals}, + "name": "E", "getter": lambda eigenstate, self, spin: eigenstate.eig}, + *extra_vars + ) + + # Now build the function that will be called for each eigenstate and will + # return the values for each variable. + def bands_wrapper(eigenstate, spin_index): + return tuple(var["getter"](eigenstate, self, spin_index) for var in all_vars) + + # Finally get the values for all coordinates involved. + coords_values = {} + for var in all_vars: + coords_values.update(var.get("coords_values", {})) + + return bands_wrapper, all_vars, coords_values + + @entry_point('wfsx file', 2) + def _read_from_wfsx(self, root_fdf, wfsx_file, extra_vars=()): + """Plots bands from the eigenvalues contained in a WFSX file. + + It also needs to get a geometry. + """ + # Get the fdf sile + fdf = self.get_sile(root_fdf or "root_fdf") + # Read the geometry from the fdf sile + self.geometry = fdf.read_geometry(output=True) + + # Get the wfsx file + wfsx_sile = self.get_sile(wfsx_file or "wfsx_file") + + # Now read all the information of the k points from the WFSX file + wfsx_info = wfsx_sile.read_info(parent=self.geometry) + # Get the number of wavefunctions in the file while performing a quick check + nwf = np.unique(wfsx_info['nwf']) + if len(nwf) > 1: + raise ValueError(f"File {wfsx_sile.file} contains different number of wavefunctions in some k points") + nwf = nwf[0] + # From the k values read in the file, build a brillouin zone object. + # We will use it just to get the linear k values for plotting. + bz = BrillouinZone(self.geometry, k=wfsx_info['k'], weight=wfsx_info['kw']) + + # Read the sizes of the file, which contain the number of spin channels + # and the number of orbitals and the number of k points. + wfsx_sizes = wfsx_sile.read_sizes() + + # Find out the spin class of the calculation. + self.spin = Spin({ + 1: Spin.UNPOLARIZED, 2: Spin.POLARIZED, + 4: Spin.NONCOLINEAR, 8: Spin.SPINORBIT + }[wfsx_sizes['nspin']]) + # Now find out how many spin channels we need. Note that if there is only + # one spin channel there will be no "spin" dimension on the final dataset. + nspin = 2 if self.spin.is_polarized else 1 + + # Get the wrapper function that we should call on each eigenstate. + # This also returns the coordinates and names to build the final dataset. + bands_wrapper, all_vars, coords_values = self._get_eigenstate_wrapper( + sisl.physics.linspace_bz(bz), extra_vars=extra_vars + ) + # Make sure all coordinates have values so that we can assume the shape + # of arrays below. + coords_values['band'] = np.arange(0, nwf) + coords_values['orb'] = np.arange(0, wfsx_sizes['nou']) + + self.ticks = None + + # Initialize all the arrays. For each quantity we will initialize + # an array of the needed shape. + arrays = {} + for var in all_vars: + # These are all the extra dimensions of the quantity. Note that a + # quantity does not need to have extra dimensions. + extra_shape = [len(coords_values[coord]) for coord in var['coords']] + # First two dimensions will always be the spin channel and the k index. + # Then add potential extra dimensions. + shape = (nspin, len(bz), *extra_shape) + # Initialize the array. + arrays[var['name']] = np.empty(shape, dtype=var.get('dtype', np.float64)) + + # Loop through eigenstates in the WFSX file and add their contribution to the bands + ik = -1 + for eigenstate in wfsx_sile.yield_eigenstate(self.geometry): + spin = eigenstate.info.get("spin", 0) + # Every time we encounter spin 0, we are in a new k point. + if spin == 0: + ik +=1 + if ik == 0: + # If this is the first eigenstate we read, get the wavefunction + # indices. We will assume that ALL EIGENSTATES have the same indices. + # Note that we already checked previously that they all have the same + # number of wfs, so this is a fair assumption. + coords_values['band'] = eigenstate.info['index'] + + # Get all the values for this eigenstate. + returns = bands_wrapper(eigenstate, spin_index=spin) + # And store them in the respective arrays. + for var, vals in zip(all_vars, returns): + arrays[var['name']][spin, ik] = vals + + # Now that we have all the values, just build the dataset. + self.bands_data = xr.Dataset( + data_vars={ + var['name']: (("spin", "k", *var['coords']), arrays[var['name']]) + for var in all_vars + } + ).assign_coords(coords_values) + + self.bands_data.attrs = {"ticks": None, "ticklabels": None, "parent": bz} + + @entry_point('band structure', 3) def _read_from_H(self, band_structure, extra_vars=()): """ Uses a sisl's `BandStructure` object to calculate the bands. """ - import xarray as xr if band_structure is None: raise ValueError("No band structure (k points path) was provided") @@ -337,38 +533,17 @@ def _read_from_H(self, band_structure, extra_vars=()): self.ticks = band_structure.lineartick() - # In case it is a non_colinear or spin-orbit calculation we will get the spin moments - if not self.spin.is_diagonal: - def _spin_moment_getter(eigenstate, plot, spin): - return eigenstate.spin_moment().real - - extra_vars = ({ - "coords": ("band", "axis"), "coords_values": dict(axis=["x", "y", "z"]), - "name": "spin_moments", "getter": _spin_moment_getter}, - *extra_vars) - - def bands_wrapper(eigenstate, spin_index): - returns = [] - for extra_var in extra_vars: - returns.append(extra_var["getter"](eigenstate, self, spin_index)) - return (eigenstate.eig, *returns) - - # Define the available spins - spin_indices = [0] - if self.spin.is_polarized: - spin_indices = [0, 1] - - # Get the eigenstates for all the available spin components - bands_arrays = [] - name = ["E"] - coords = [('band'), ] - coords_values = {"spin": spin_indices, "k": band_structure.lineark()} - for extra_var in extra_vars: - name.append(extra_var["name"]) - coords.append(extra_var["coords"]) - coords_values.update(extra_var.get("coords_values", {})) + # Get the wrapper function that we should call on each eigenstate. + # This also returns the coordinates and names to build the final dataset. + bands_wrapper, all_vars, coords_values= self._get_eigenstate_wrapper( + band_structure.lineark(), extra_vars=extra_vars + ) - for spin_index in spin_indices: + # Get a dataset with all values for all spin indices + spin_datasets = [] + coords = [var['coords'] for var in all_vars] + name = [var['name'] for var in all_vars] + for spin_index in coords_values['spin']: # Non collinear routines don't accept the keyword argument "spin" spin_kwarg = {"spin": spin_index} @@ -382,31 +557,15 @@ def bands_wrapper(eigenstate, spin_index): coords=coords, name=name, ) - bands_arrays.append(spin_bands) + spin_datasets.append(spin_bands) # Merge everything into a single dataset with a spin dimension - self.bands_data = xr.concat(bands_arrays, "spin").assign_coords(coords_values) + self.bands_data = xr.concat(spin_datasets, "spin").assign_coords(coords_values) - self.bands['k'] = band_structure.lineark() # Inform of where to place the ticks - self.bands_data.attrs = {"ticks": self.ticks[0], "ticklabels": self.ticks[1], **bands_arrays[0].attrs} - - @entry_point('bands file', 0) - def _read_siesta_output(self, bands_file, band_structure): - """ - Reads the bands information from a SIESTA bands file. - """ - if band_structure: - raise ValueError("A path was provided, therefore we can not use the .bands file even if there is one") - - self.bands_data = self.get_sile(bands_file or "bands_file").read_data(as_dataarray=True) - - # Define the spin class of the results we have retrieved - if len(self.bands_data.spin.values) == 2: - self.spin = sisl.Spin("p") + self.bands_data.attrs = {"ticks": self.ticks[0], "ticklabels": self.ticks[1], **spin_datasets[0].attrs} def _after_read(self): - import xarray as xr if isinstance(self.bands_data, xr.DataArray): attrs = self.bands_data.attrs self.bands_data = xr.Dataset({"E": self.bands_data}) diff --git a/sisl/viz/plots/fatbands.py b/sisl/viz/plots/fatbands.py index 9e0365beb0..8991267cb9 100644 --- a/sisl/viz/plots/fatbands.py +++ b/sisl/viz/plots/fatbands.py @@ -4,6 +4,7 @@ import numpy as np import sisl +from sisl.physics.spin import Spin from ..plot import entry_point from .bands import BandsPlot from ..plotutils import random_color @@ -126,20 +127,6 @@ class FatbandsPlot(BandsPlot): _parameters = ( - SileInput(key='wfsx_file', name='Path to WFSX file', - dtype=sisl.io.siesta.wfsxSileSiesta, - default=None, - help="""The WFSX file to get the weights of the different orbitals in the bands. - In standard SIESTA nomenclature, this should be the *.bands.WFSX file, as it is the one - that contains the weights that correspond to the bands. - - This file is only meaningful (and required) if fatbands are plotted from the .bands file. - Otherwise, the bands and weights will be generated from the hamiltonian by sisl. - If the *.bands file is provided but the wfsx one isn't, we will try to find it. - If `bands_file` is SystemLabel.bands, we will look for SystemLabel.bands.WFSX - """ - ), - FloatInput(key='scale', name='Scale factor', default=None, help="""The factor by which the width of all fatbands should be multiplied. @@ -190,83 +177,21 @@ class FatbandsPlot(BandsPlot): def weights(self): return self.bands_data["weight"] - @entry_point("siesta output", 0) - def _read_siesta_output(self, wfsx_file, bands_file, root_fdf): + @entry_point("wfsx file", 0) + def _read_from_wfsx(self, root_fdf, wfsx_file): """Generates fatbands from SIESTA output. - Uses the `.bands` file to read the bands and a `.wfsx` file to - retrieve the wavefunctions coefficients. + Uses the `.wfsx` file to retrieve the eigenstates. From them, it computes + all the needed quantities (eigenvalues, orbital contribution, ...). """ - from xarray import DataArray, Dataset - - # Try to get the wfsx file either by user input or by guessing it - # from bands_file - bands_file = self.get_sile(bands_file or "bands_file").file - if wfsx_file is None: - wfsx_file = bands_file.with_suffix(bands_file.suffix + ".WFSX") - - # We will need the overlap matrix from the hamiltonian to get the correct - # weights. - # If there is no root_fdf we will try to guess it from bands_file - if root_fdf is None and not hasattr(self, "H"): - possible_fdf = bands_file.with_suffix(".fdf") - print(f"We are assuming that the fdf associated to {bands_file} is {possible_fdf}."+ - ' If it is not, please provide a "root_fdf" by using the update_settings method.') - self.update_settings(root_fdf=root_fdf, run_updates=False) - - self.setup_hamiltonian() - - # If the wfsx doesn't exist, we will not even bother to read the bands - if not wfsx_file.exists(): - raise ValueError(f"We did not find a WFSX file in the location {wfsx_file}") - - # Otherwise we will make BandsPlot read the bands - super()._read_siesta_output() - - # And then read the weights from the wfsx file - wfsx_sile = self.get_sile(wfsx_file) - - weights = [] - for i, state in enumerate(wfsx_sile.yield_eigenstate(self.H)): - # Each eigenstate represents all the states for a given k-point - - # Get the band indices to which these states correspond - if i == 0: - bands = state.info["indices"] - - # Get the weights for this eigenstate - weights.append(state.norm2(sum=False)) - - weights = np.array(weights).real - - # Finally, build the weights dataarray so that it can be used by _set_data - weights = DataArray( - weights, - coords={ - "k": self.bands.k, - "band": bands, - "orb": np.arange(0, weights.shape[2]), - }, - dims=("k", "band", "orb") - ) - - # Add the spin dimension so that the weights array is normalized, - # even though spin is not yet supported by this entrypoint - weights = weights.expand_dims("spin") - - # Merge everything into a dataset - attrs = self.bands_data.attrs - self.bands_data = Dataset({"E": self.bands_data, "weight": weights}) - self.bands_data.attrs = attrs - - # Set up the options for the 'groups' setting based on the plot's associated geometry - self._set_group_options() + self._entry_point_with_extra_vars(super()._read_from_wfsx) @entry_point("hamiltonian", 1) def _read_from_H(self): - """ - Calculates the fatbands from a sisl hamiltonian. - """ + """Calculates the fatbands from a sisl hamiltonian.""" + self._entry_point_with_extra_vars(super()._read_from_H) + + def _entry_point_with_extra_vars(self, entry_point): # Define the function that will "catch" each eigenstate and # build the weights array. See BandsPlot._read_from_H to understand where # this will go exactly @@ -286,7 +211,7 @@ def _weights_from_eigenstate(eigenstate, plot, spin_index): # thanks to the above step bands_read = False; err = None try: - super()._read_from_H(extra_vars=[{"coords": ("band", "orb"), "name": "weight", "getter": _weights_from_eigenstate}]) + entry_point(extra_vars=[{"coords": ("band", "orb"), "name": "weight", "getter": _weights_from_eigenstate}]) bands_read = True except Exception as e: # Let's keep this error, we are going to at least set the group options so that the @@ -307,13 +232,7 @@ def _set_group_options(self): if band_struct is not None: self.geometry = band_struct.parent.geometry - if getattr(self, "H", None) is not None: - spin = self.H.spin - else: - # There is yet no spin support reading from bands.WFSX - spin = sisl.Spin.UNPOLARIZED - - self.get_param('groups').update_options(self.geometry, spin) + self.get_param('groups').update_options(self.geometry, self.spin) def _set_data(self): # We get the information that the Bandsplot wants to send to the drawer diff --git a/sisl/viz/plots/grid.py b/sisl/viz/plots/grid.py index c399381d9f..f513baa3f1 100644 --- a/sisl/viz/plots/grid.py +++ b/sisl/viz/plots/grid.py @@ -1264,6 +1264,10 @@ class WavefunctionPlot(GridPlot): eigenstate: EigenstateElectron, optional The eigenstate that contains the coefficients of the wavefunction. Note that an eigenstate can contain coefficients for multiple states. + wfsx_file: wfsxSileSiesta, optional + Siesta WFSX file to directly read the coefficients from. + If the root_fdf file is provided but the wfsx one isn't, we will try + to find it as SystemLabel.WFSX. geometry: Geometry, optional Necessary to generate the grid and to plot the wavefunctions, since the basis orbitals are needed. If you provide a @@ -1398,6 +1402,15 @@ class WavefunctionPlot(GridPlot): """ ), + SileInput(key='wfsx_file', name='Path to WFSX file', + dtype=sisl.io.siesta.wfsxSileSiesta, + default=None, + help="""Siesta WFSX file to directly read the coefficients from. + If the root_fdf file is provided but the wfsx one isn't, we will try to find it + as SystemLabel.WFSX. + """ + ), + SislObjectInput(key='geometry', name='Geometry', default=None, dtype=sisl.Geometry, @@ -1450,7 +1463,27 @@ def _read_nosource(self, eigenstate): self.eigenstate = eigenstate - @entry_point('hamiltonian', 1) + @entry_point('Siesta WFSX file', 1) + def _read_from_WFSX_file(self, wfsx_file, k, spin, root_fdf): + """Reads the wavefunction coefficients from a SIESTA WFSX file""" + # Try to read the geometry + fdf = self.get_sile(root_fdf or "root_fdf") + if fdf is None: + raise ValueError("The setting 'root_fdf' needs to point to an fdf file with a geometry") + geometry = fdf.read_geometry(output=True) + + # Get the WFSX file. If not provided, it is inferred from the fdf. + wfsx = self.get_sile(wfsx_file or "wfsx_file") + if not wfsx.file.exists(): + raise ValueError(f"File '{wfsx.file}' does not exist.") + + # Try to find the eigenstate that we need + self.eigenstate = wfsx.read_eigenstate(k=k, spin=spin[0], parent=geometry) + if self.eigenstate is None: + # We have not found it. + raise ValueError(f"A state with k={k} was not found in file {wfsx.file}.") + + @entry_point('hamiltonian', 2) def _read_from_H(self, k, spin): """ Calculates the eigenstates from a Hamiltonian and then generates the wavefunctions. @@ -1464,6 +1497,23 @@ def _after_read(self): # calling it later in _set_data pass + def _get_eigenstate(self, i): + + if "index" in self.eigenstate.info: + wf_i = np.nonzero(self.eigenstate.info["index"] == i)[0] + if len(wf_i) == 0: + raise ValueError(f"Wavefunction with index {i} is not present in the eigenstate. Available indices: {self.eigenstate.info['index']}." + f"Entry point used: {self.source._name}") + wf_i = wf_i[0] + else: + max_index = len(self.eigenstate) + if i > max_index: + raise ValueError(f"Wavefunction with index {i} is not present in the eigenstate. Available range: [0, {max_index}]." + f"Entry point used: {self.source._name}") + wf_i = i + + return self.eigenstate[wf_i] + def _set_data(self, i, geometry, grid, k, grid_prec, nsc): if geometry is not None: @@ -1498,16 +1548,18 @@ def _set_data(self, i, geometry, grid, k, grid_prec, nsc): tiled_geometry = tiled_geometry.tile(sc_i, ax) nsc[ax] = 1 + is_gamma = (np.array(k) == 0).all() if grid is None: - dtype = np.float64 if (np.array(k) == 0).all() else np.complex128 + dtype = np.float64 if is_gamma else np.complex128 self.grid = sisl.Grid(grid_prec, geometry=tiled_geometry, dtype=dtype) # GridPlot's after_read basically sets the x_range, y_range and z_range options # which need to know what the grid is, that's why we are calling it here super()._after_read() - self.eigenstate[i].wavefunction(self.grid) + state = self._get_eigenstate(i) + state.wavefunction(self.grid) - return super()._set_data(nsc=nsc) + return super()._set_data(nsc=nsc, trace_name=f"WF {i} ({state.eig[0]:.2f} eV)") GridPlot.backends.register_child(WavefunctionPlot.backends) diff --git a/sisl/viz/plots/pdos.py b/sisl/viz/plots/pdos.py index 81aaf916bc..05524b09fd 100644 --- a/sisl/viz/plots/pdos.py +++ b/sisl/viz/plots/pdos.py @@ -34,6 +34,11 @@ class PdosPlot(Plot): tbt_nc: tbtncSileTBtrans, optional This parameter explicitly sets a .TBT.nc file. Otherwise, the PDOS file is attempted to read from the fdf file + wfsx_file: wfsxSileSiesta, optional + The WFSX file to get the eigenstates. In standard SIESTA + nomenclature, this should probably be the *.fullBZ.WFSX file, as it + is the one that contains the eigenstates from the full + brillouin zone. geometry: Geometry or sile (or path to file) that contains a geometry, optional If this is passed, the geometry that has been read is ignored and this one is used instead. @@ -122,6 +127,15 @@ class PdosPlot(Plot): help = """This parameter explicitly sets a .TBT.nc file. Otherwise, the PDOS file is attempted to read from the fdf file """ ), + SileInput(key='wfsx_file', name='Path to WFSX file', + dtype=sisl.io.siesta.wfsxSileSiesta, + default=None, + help="""The WFSX file to get the eigenstates. + In standard SIESTA nomenclature, this should probably be the *.fullBZ.WFSX file, as it is the one + that contains the eigenstates from the full brillouin zone. + """ + ), + GeometryInput( key = "geometry", name = "Geometry to force on the plot", group="dataread", @@ -300,7 +314,7 @@ def _read_siesta_output(self, pdos_file): #Get the info from the .PDOS file self.geometry, self.E, self.PDOS = self.get_sile(pdos_file or "pdos_file").read_data() - @entry_point("TB trans", 1) + @entry_point("TB trans", 2) def _read_TBtrans(self, root_fdf, tbt_nc): """ Reads the PDOS from a *.TBT.nc file coming from a TBtrans run. @@ -320,7 +334,43 @@ def _read_TBtrans(self, root_fdf, tbt_nc): # Read the geometry from the TBT.nc file and get only the device part self.geometry = tbt_sile.read_geometry(**read_geometry_kwargs).sub(tbt_sile.a_dev) - @entry_point('hamiltonian', 2) + @entry_point('wfsx file', 3) + def _read_from_wfsx(self, root_fdf, wfsx_file, Erange, nE, E0, distribution): + """Generates the PDOS values from a file containing eigenstates.""" + # Read the hamiltonian. We need it because we need the overlap matrix. + if not hasattr(self, "H"): + self.setup_hamiltonian() + + # Get the wfsx file + wfsx_sile = self.get_sile(wfsx_file or "wfsx_file") + + # Read the sizes of the file, which contain the number of spin channels + # and the number of orbitals and the number of k points. + sizes = wfsx_sile.read_sizes() + # Check that spin sizes of hamiltonian and wfsx file match + assert self.H.spin.size == sizes['nspin'], \ + f"Hamiltonian has spin size {self.H.spin.size} while file has spin size {sizes['nspin']}" + # Get the size of the spin channel. The size returned might be 8 if it is a spin-orbit + # calculation, but we need only 4 spin channels (total, x, y and z), same as with non-colinear + nspin = min(4, sizes['nspin']) + + # Get the energies for which we need to calculate the PDOS. + self.E = np.linspace(Erange[0], Erange[-1], nE) + E0 + + # Initialize the PDOS array + self.PDOS = np.zeros((nspin, sizes["nou"], self.E.shape[0]), dtype=np.float64) + + # Loop through eigenstates in the WFSX file and add their contribution to the PDOS. + # Note that we pass the hamiltonian as the parent here so that the overlap matrix + # for each point can be calculated by eigenstate.PDOS() + for eigenstate in wfsx_sile.yield_eigenstate(self.H): + spin = eigenstate.info.get("spin", 0) + if nspin == 4: + spin = slice(None) + + self.PDOS[spin] += eigenstate.PDOS(self.E, distribution=distribution) * eigenstate.info.get("weight", 1) + + @entry_point('hamiltonian', 4) def _read_from_H(self, kgrid, kgrid_displ, Erange, nE, E0, distribution): """ Calculates the PDOS from a sisl Hamiltonian. @@ -338,7 +388,7 @@ def _read_from_H(self, kgrid, kgrid_displ, Erange, nE, E0, distribution): self.E = np.linspace(Erange[0], Erange[-1], nE) + E0 - self.mp = sisl.MonkhorstPack(self.H, kgrid, kgrid_displ) + self.bz = sisl.MonkhorstPack(self.H, kgrid, kgrid_displ) # Define the available spins spin_indices = [0] @@ -348,7 +398,7 @@ def _read_from_H(self, kgrid, kgrid_displ, Erange, nE, E0, distribution): # Calculate the PDOS for all available spins PDOS = [] for spin in spin_indices: - with self.mp.apply(pool=_do_parallel_calc) as parallel: + with self.bz.apply(pool=_do_parallel_calc) as parallel: spin_PDOS = parallel.average.eigenstate( spin=spin, wrap=lambda eig: eig.PDOS(self.E, distribution=distribution) diff --git a/sisl/viz/plots/tests/test_bands.py b/sisl/viz/plots/tests/test_bands.py index 5cc120fab6..86dae1d9d7 100644 --- a/sisl/viz/plots/tests/test_bands.py +++ b/sisl/viz/plots/tests/test_bands.py @@ -76,8 +76,11 @@ def init_func_and_attrs(self, request, siesta_test_files): # from two different prespectives if name.startswith("sisl_H_path"): # Passing a list of points (as if we were interacting from a GUI) + # We want 6 points in total. This is the path that we want to get: + # [0,0,0] --2-- [2/3, 1/3, 0] --1-- [1/2, 0, 0] path = [{"active": True, "x": x, "y": y, "z": z, "divisions": 3, "name": tick} for tick, (x, y, z) in zip(["Gamma", "M", "K"], [[0, 0, 0], [2/3, 1/3, 0], [1/2, 0, 0]])] + path[-1]['divisions'] = 2 init_func = partial(H.plot.bands, band_structure=path) else: