wham_get_deleig_a Subroutine

  subroutine wham_get_deleig_a(deleig_a, eig, delHH_a, UU, num_wann, pw90_band_deriv_degen, &
                               error, comm)
    !================================================!
    !
    !! Band derivatives dE/dk_a
    !
    !================================================!

    use w90_constants, only: dp !, cmplx_0, cmplx_i
    use w90_utility, only: utility_diagonalize, utility_rotate, utility_rotate_diag
    use w90_postw90_types, only: pw90_band_deriv_degen_type
    use w90_comms, only: w90_comm_type

    ! arguments
    type(pw90_band_deriv_degen_type), intent(in) :: pw90_band_deriv_degen
    integer, intent(in) :: num_wann
    real(kind=dp), intent(in) :: eig(num_wann)
    real(kind=dp), intent(out) :: deleig_a(num_wann)
    complex(kind=dp), intent(in) :: delHH_a(:, :)
    complex(kind=dp), intent(in) :: UU(:, :)
    type(w90_error_type), allocatable, intent(out) :: error
    type(w90_comm_type), intent(in) :: comm

    ! local variables
    integer                       :: i, degen_min, degen_max, dim
    real(kind=dp)                 :: diff
    complex(kind=dp), allocatable :: delHH_bar_a(:, :), U_deg(:, :)

    allocate (delHH_bar_a(num_wann, num_wann))
    allocate (U_deg(num_wann, num_wann))

    if (pw90_band_deriv_degen%use_degen_pert) then

      delHH_bar_a = utility_rotate(delHH_a, UU, num_wann)

      ! Assuming that the energy eigenvalues are stored in eig(:) in
      ! increasing order (diff >= 0)

      i = 0
      do
        i = i + 1
        if (i > num_wann) exit
        if (i + 1 <= num_wann) then
          diff = eig(i + 1) - eig(i)
        else

          ! i-th is the highest band, and it is non-degenerate

          diff = pw90_band_deriv_degen%degen_thr + 1.0_dp
        end if
        if (diff < pw90_band_deriv_degen%degen_thr) then

          ! Bands i and i+1 are degenerate

          degen_min = i
          degen_max = degen_min + 1

          ! See if any higher bands are in the same degenerate group

          do
            if (degen_max + 1 > num_wann) exit
            diff = eig(degen_max + 1) - eig(degen_max)
            if (diff < pw90_band_deriv_degen%degen_thr) then
              degen_max = degen_max + 1
            else
              exit
            end if
          end do

          ! Bands from degen_min to degen_max are degenerate. Diagonalize
          ! the submatrix in Eq.(31) YWVS07 over this degenerate subspace.
          ! The eigenvalues are the band gradients

          dim = degen_max - degen_min + 1
          call utility_diagonalize(delHH_bar_a(degen_min:degen_max, degen_min:degen_max), dim, &
                                   deleig_a(degen_min:degen_max), U_deg(1:dim, 1:dim), error, comm)
          if (allocated(error)) return

          ! Scanned bands up to degen_max

          i = degen_max
        else

          ! Use non-degenerate form [Eq.(27) YWVS07] for current (i-th) band

          deleig_a(i) = real(delHH_bar_a(i, i), dp)
        end if
      end do

    else

      ! Use non-degenerate form for all bands

      deleig_a(:) = real(utility_rotate_diag(delHH_a(:, :), UU, num_wann), dp)

    end if

  end subroutine wham_get_deleig_a