plot_read_unk Subroutine

  subroutine plot_read_unk(wfnname, formatted, spinors, inc_band, num_inc, num_bands, &
                           kpt, ngx, ngy, ngz, wvfn, wvfn_nc, stdout, error, comm)
    !! Read a single UNK file into `wvfn` (or `wvfn_nc` for spinors), keeping
    !! the first `num_inc` bands flagged in `inc_band` (bands not flagged are
    !! read and discarded).
    !!
    !! Failures set the error and return immediately. This is safe even though
    !! the caller distributes k-points over ranks: set_error_file synchronises
    !! the failing rank via comms_sync_error, and the other ranks pick the
    !! error up in the matching handshake at their next synchronised
    !! collective (e.g. the comms_reduce after the caller's k-point loop).

    use w90_constants, only: dp
    use w90_comms, only: w90_comm_type
    use w90_error, only: w90_error_type, set_error_file

    implicit none

    character(len=*), intent(in) :: wfnname
    !! Name of the UNK file to read
    logical, intent(in) :: formatted
    !! Whether the UNK file is formatted
    logical, intent(in) :: spinors
    !! Whether the wavefunctions are spinors
    logical, intent(in) :: inc_band(:)
    !! Bands to keep
    integer, intent(in) :: num_inc
    !! Number of bands to keep
    integer, intent(in) :: num_bands
    !! Total number of bands in the file
    integer, intent(in) :: kpt
    !! Expected k-point index (checked against the file header)
    integer, intent(in) :: ngx, ngy, ngz
    !! Expected grid dimensions (checked against the file header)
    complex(kind=dp), allocatable, intent(inout) :: wvfn(:, :)
    !! Destination, band as second index (caller-allocated; if .not. spinors)
    complex(kind=dp), intent(inout), allocatable :: wvfn_nc(:, :, :)
    !! Spinor destination, spin as third index (caller-allocated; if spinors)
    integer, intent(in) :: stdout
    type(w90_error_type), allocatable, intent(out) :: error
    type(w90_comm_type), intent(in) :: comm

    integer :: file_unit, ierr, ix, iy, iz, ik, nbnd, counter, loop_b, ngpts

    logical :: have_file

    ngpts = ngx*ngy*ngz

    ! A missing UNK file must be detected here rather than silently
    ! (re)created by the open below (status='old' forbids creation).
    inquire (file=wfnname, exist=have_file)
    if (.not. have_file) then
      call set_error_file(error, 'plot_wannier: file '//wfnname//' not found', comm)
      return
    end if

    if (formatted) then
      open (newunit=file_unit, file=wfnname, form='formatted', status='old', iostat=ierr)
    else
      open (newunit=file_unit, file=wfnname, form='unformatted', status='old', iostat=ierr)
    end if
    if (ierr /= 0) then
      call set_error_file(error, 'plot_wannier: could not open file '//wfnname, comm)
      return
    end if

    if (formatted) then
      read (file_unit, *, iostat=ierr) ix, iy, iz, ik, nbnd
    else
      read (file_unit, iostat=ierr) ix, iy, iz, ik, nbnd
    end if
    if (ierr /= 0) then
      call set_error_file(error, 'plot_wannier: error reading file '//wfnname, comm)
      close (file_unit)
      return
    end if

    if ((ix /= ngx) .or. (iy /= ngy) .or. (iz /= ngz) .or. (ik /= kpt)) then
      write (stdout, '(1x,a,a)') 'WARNING: mismatch in file', trim(wfnname)
      write (stdout, '(1x,5(a6,I5))') '   ix=', ix, '   iy=', iy, '   iz=', iz, '   ik=', ik, ' nbnd=', nbnd
      write (stdout, '(1x,5(a6,I5))') '  ngx=', ngx, '  ngy=', ngy, '  ngz=', ngz, '  kpt=', kpt, 'bands=', num_bands
      call set_error_file(error, 'plot_wannier: mismatch in file '//wfnname, comm)
      close (file_unit)
      return
    end if

    counter = 1
    do loop_b = 1, num_bands
      if (counter > num_inc) exit
      if (.not. spinors) then
        call plot_read_unk_band(wvfn(:, counter))
      else
        call plot_read_unk_band(wvfn_nc(:, counter, 1)) ! up-spinor
        if (ierr == 0) call plot_read_unk_band(wvfn_nc(:, counter, 2)) ! down-spinor
      end if
      if (ierr /= 0) then
        call set_error_file(error, 'plot_wannier: error reading file '//wfnname, comm)
        close (file_unit)
        return
      end if
      if (inc_band(loop_b)) counter = counter + 1
    end do

    close (file_unit)

  contains

    subroutine plot_read_unk_band(band)
      !! Read one band's worth of grid values, leaving any failure in `ierr`
      complex(kind=dp), intent(out) :: band(:)

      real(kind=dp) :: w_real, w_imag
      integer :: nx

      if (formatted) then
        do nx = 1, ngpts
          read (file_unit, *, iostat=ierr) w_real, w_imag
          if (ierr /= 0) return
          band(nx) = cmplx(w_real, w_imag, kind=dp)
        end do
      else
        read (file_unit, iostat=ierr) (band(nx), nx=1, ngpts)
      end if
    end subroutine plot_read_unk_band

  end subroutine plot_read_unk