conv_read_spn Subroutine

  subroutine conv_read_spn(stdout, seedname)
    !================================================!
    !
    !! Read unformatted spn file
    !
    !================================================!

    use w90_constants, only: eps6, dp
    use w90spn_parameters, only: num_bands, num_kpts

    implicit none

    integer, intent(in) :: stdout
    character(len=50), intent(in)  :: seedname

    integer :: spn_unit, m, n, ik, ierr, s, counter
    complex(kind=dp), allocatable :: spn_temp(:, :)

    write (stdout, '(3a)') 'Reading information from unformatted file ', trim(seedname), '.spn :'

    open (newunit=spn_unit, file=trim(seedname)//'.spn', status='old', form='unformatted', err=109)

    ! Read comment line
    read (spn_unit, err=110, end=110) header
    header = ADJUSTL(header)
    write (stdout, '(1x,a)') trim(header)

    ! Consistency checks
    read (spn_unit, err=110, end=110) num_bands, num_kpts
    write (stdout, '(1x,a,i0)') "Number of bands: ", num_bands
    write (stdout, '(1x,a,i0)') "Number of k-points: ", num_kpts

    allocate (spn_o(num_bands, num_bands, num_kpts, 3), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating spm_temp in conv_read_spn', stdout)

    allocate (spn_temp(3, (num_bands*(num_bands + 1))/2), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating spm_temp in conv_read_spn', stdout)
    do ik = 1, num_kpts
      read (spn_unit) ((spn_temp(s, m), s=1, 3), m=1, (num_bands*(num_bands + 1))/2)
      counter = 0
      do m = 1, num_bands
        do n = 1, m
          counter = counter + 1
          spn_o(n, m, ik, 1) = spn_temp(1, counter)
          spn_o(n, m, ik, 2) = spn_temp(2, counter)
          spn_o(n, m, ik, 3) = spn_temp(3, counter)
          ! Although each diagonal element of spn_o should be a real number,
          ! actually it has a very small imaginary part.
          ! We skip the conjugation on the diagonal elements so that
          ! the file after formatted <==> unformatted conversions is exactly
          ! the same as the original file, otherwise the diagonal elements
          ! are the conjugations of those of the original file.
          if (m == n) cycle
          spn_o(m, n, ik, 1) = conjg(spn_temp(1, counter))
          spn_o(m, n, ik, 2) = conjg(spn_temp(2, counter))
          spn_o(m, n, ik, 3) = conjg(spn_temp(3, counter))
        end do
      end do
    end do

    close (spn_unit)
    write (stdout, '(1x,a)') "spn: read."

    deallocate (spn_temp, stat=ierr)
    if (ierr /= 0) call io_error('Error in deallocating spm_temp in conv_read_spn', stdout)

    write (stdout, '(1x,a)') 'read done.'

    return

109 call io_error('Error opening '//trim(seedname)//'.spn.fmt in conv_read_spn', stdout)
110 call io_error('Error reading '//trim(seedname)//'.spn.fmt in conv_read_spn', stdout)

  end subroutine conv_read_spn