gyrotropic_get_curv_w_k Subroutine

  subroutine gyrotropic_get_curv_w_k(eig, AA, curv_w_k, pw90_gyrotropic)
    !================================================!
    !
    ! calculation of the band-resolved
    ! frequency-dependent   berry curvature
    !
    ! tildeOmega(w)=
    !     -eps_{bcd}sum_m ( w_mn^2/(wmn^2-w^2)) *Im[A_{nm,c}A_{mn,d}
    !
    !================================================!

    use w90_postw90_types, only: pw90_gyrotropic_type
    use w90_constants, only: dp

    implicit none

    ! arguments
    type(pw90_gyrotropic_type), intent(in) :: pw90_gyrotropic
    real(kind=dp), intent(in) :: eig(:)
    real(kind=dp), intent(out) :: curv_w_k(:, :, :)  ! (num_wann,n_freq,3)
    complex(kind=dp), intent(in) :: AA(:, :, :)

    ! local variables
    real(kind=dp), allocatable    :: multWre(:)
    integer          :: i, n, m, n1, m1
    real(kind=dp)    :: wmn

    allocate (multWre(pw90_gyrotropic%nfreq))
    curv_w_k(:, :, :) = 0_dp

    do n1 = 1, pw90_gyrotropic%num_bands
      n = pw90_gyrotropic%band_list(n1)
      do m1 = 1, pw90_gyrotropic%num_bands
        m = pw90_gyrotropic%band_list(m1)
        if (n == m) cycle
        wmn = eig(m) - eig(n)
        multWre(:) = real(wmn**2/(wmn**2 - pw90_gyrotropic%freq_list(:)**2))
        do i = 1, 3
          curv_w_k(n, :, i) = curv_w_k(n, :, i) - &
                              2_dp*aimag(AA(n, m, alpha_A(i))*AA(m, n, beta_A(i)))*multWre
        end do
      end do !m
    end do !n

  end subroutine gyrotropic_get_curv_w_k