wham_get_deleig_a Subroutine

  subroutine wham_get_deleig_a(deleig_a, eig, delHH_a, UU)
    !==========================!
    !                          !
    !! 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_parameters, only: num_wann, use_degen_pert, degen_thr

    ! Arguments
    !
    real(kind=dp), intent(out) :: deleig_a(num_wann)
    real(kind=dp), intent(in)  :: eig(num_wann)
    complex(kind=dp), dimension(:, :), intent(in)  :: delHH_a
    complex(kind=dp), dimension(:, :), intent(in)  :: UU

    ! Misc/Dummy
    !
    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 (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 = degen_thr + 1.0_dp
        end if
        if (diff < 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 < 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))
          !
          ! 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