gyrotropic_get_NOA_Bnl_orb Subroutine

  subroutine gyrotropic_get_NOA_Bnl_orb(eig, del_eig, AA, num_occ, occ_list, num_unocc, &
                                        unocc_list, Bnl, pw90_gyrotropic)
    !================================================!
    !
    ! Calculating the matrix
    ! B_{nl,ac}(num_occ,num_unocc,3,3)=
    !      -sum_m(  (E_m-E_n)A_nma*Amlc +(E_l-E_m)A_nmc*A_mla -
    !      -i( del_a (E_n+E_l) A_nlc
    !   in units eV*Ang^2
    !================================================!

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

    implicit none

    ! arguments
    type(pw90_gyrotropic_type), intent(in) :: pw90_gyrotropic
    integer, intent(in) :: num_occ, num_unocc
    integer, intent(in) :: occ_list(:), unocc_list(:)
    real(kind=dp), intent(in) :: eig(:)               !  n
    real(kind=dp), intent(in) :: del_eig(:, :)        !  n
    complex(kind=dp), intent(in) :: AA(:, :, :)       !  n,l,a
    complex(kind=dp), intent(out) :: Bnl(:, :, :, :)  !  n,l,a,c

    ! local variables
    integer n, m, l, a, c, n1, m1, l1
    Bnl(:, :, :, :) = cmplx_0

    do a = 1, 3
      do c = 1, 3
        do n1 = 1, num_occ
          n = occ_list(n1)
          do l1 = 1, num_unocc
            l = unocc_list(l1)
            Bnl(n1, l1, a, c) = -cmplx_i*(del_eig(n, a) + del_eig(l, a))*AA(n, l, c)
            do m1 = 1, pw90_gyrotropic%num_bands
              m = pw90_gyrotropic%band_list(m1)
              Bnl(n1, l1, a, c) = Bnl(n1, l1, a, c) + &
                                  (eig(n) - eig(m))*AA(n, m, a)*AA(m, l, c) - &
                                  (eig(l) - eig(m))*AA(n, m, c)*AA(m, l, a)
            end do ! m1
          end do !l1
        end do !n1
      end do !c
    end do !a

  end subroutine gyrotropic_get_NOA_Bnl_orb