gyrotropic_get_curv_w_k Subroutine

  subroutine gyrotropic_get_curv_w_k(eig, AA, curv_w_k)
    !======================================================================!
    !                                                                      !
    ! 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_constants, only: dp
    use w90_parameters, only: num_wann, gyrotropic_nfreq, gyrotropic_freq_list, &
      gyrotropic_band_list, gyrotropic_num_bands
    ! Arguments
    !
    real(kind=dp), intent(in)                    :: eig(:)
    complex(kind=dp), intent(in)                 :: AA(:, :, :)
    real(kind=dp), dimension(:, :, :), intent(out) :: curv_w_k  ! (num_wann,n_freq,3)

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

    allocate (multWre(gyrotropic_nfreq))
    curv_w_k(:, :, :) = 0_dp

    do n1 = 1, gyrotropic_num_bands
      n = gyrotropic_band_list(n1)
      do m1 = 1, gyrotropic_num_bands
        m = gyrotropic_band_list(m1)
        if (n == m) cycle
        wmn = eig(m) - eig(n)
        multWre(:) = real(wmn**2/(wmn**2 - 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
        enddo
      enddo !m
    enddo !n

  end subroutine gyrotropic_get_curv_w_k