subroutine gyrotropic_get_NOA_k(kpt, kweight, eig, del_eig, AA, UU, gyro_NOA_orb, gyro_NOA_spn)
!====================================================================!
! !
! Contribution from point k to the real (antisymmetric) part !
! of the natural complex interband optical conductivity !
! !
! Re gyro_NOA_orb = SUM_{n,l}^{oe} hbar^{-2}/(w_nl^2-w^2) *
! Re ( A_lnb Bnlac -Alna Bnlbc)
! -SUM_{n,l}^{oe} hbar^{-2}(3*w_ln^2-w^2)/(w_nl^2-w^2)^2 *
! Im ( A_lnb Bnlac -Alna Bnlac)nm_a A_mn_b )
! [units of Ang^3/eV] !
! [units of Ang^3] !
! Re gyro_NOA_spn_{ab,c} = SUM_{n,l}^{oe} hbar^{-2}/(w_nl^2-w^2) *
! Re ( A_lnb Bnlac -Alna Bnlbc)
! [units of Ang/eV^2] !
!
! here a,b defined as epsilon_{abd}=1 (and NOA_dc tensor is saved) !
!====================================================================!
use w90_constants, only: dp, cmplx_0, cmplx_i, pi, cmplx_1
use w90_utility, only: utility_rotate
use w90_parameters, only: num_wann, gyrotropic_nfreq, gyrotropic_freq_list, &
fermi_energy_list, nfermi, gyrotropic_eigval_max, &
gyrotropic_num_bands, gyrotropic_band_list, iprint
use w90_comms, only: on_root
use w90_io, only: stdout, io_time, io_error
use w90_postw90_common, only: pw90common_fourier_R_to_k_new
use w90_get_oper, only: SS_R
use w90_spin, only: spin_get_S
! Arguments
!
real(kind=dp), intent(in) :: kpt(3), kweight
real(kind=dp), dimension(:, :, :, :), optional, intent(inout) :: gyro_NOA_spn
real(kind=dp), dimension(:, :, :, :), intent(inout) :: gyro_NOA_orb
complex(kind=dp), dimension(:, :, :), intent(in) :: AA
complex(kind=dp), dimension(:, :), intent(in) :: UU
real(kind=dp), dimension(:), intent(in) :: eig
real(kind=dp), dimension(:, :), intent(in) :: del_eig
complex(kind=dp), allocatable :: Bnl_orb(:, :, :, :)
complex(kind=dp), allocatable :: Bnl_spin(:, :, :, :)
complex(kind=dp), allocatable :: SS(:, :, :)
complex(kind=dp), allocatable :: S_h(:, :, :)
complex(kind=dp) :: multW1(gyrotropic_nfreq)
real(kind=dp) :: multWe(gyrotropic_nfreq), multWm(gyrotropic_nfreq)
integer :: num_occ, num_unocc, occ_list(num_wann), unocc_list(num_wann)
real(kind=dp) :: wmn
integer :: i, j, n, l, n1, l1, a, b, c, ab, ifermi
real(kind=dp) :: wln
if (present(gyro_NOA_spn)) then
allocate (SS(num_wann, num_wann, 3))
allocate (S_h(num_wann, num_wann, 3))
do j = 1, 3 ! spin direction
call pw90common_fourier_R_to_k_new(kpt, SS_R(:, :, :, j), OO=SS(:, :, j))
S_h(:, :, j) = utility_rotate(SS(:, :, j), UU, num_wann)
enddo
endif
do ifermi = 1, nfermi
num_occ = 0
num_unocc = 0
do n1 = 1, gyrotropic_num_bands
n = gyrotropic_band_list(n1)
if (eig(n) < fermi_energy_list(ifermi)) then
num_occ = num_occ + 1
occ_list(num_occ) = n
elseif (eig(n) < gyrotropic_eigval_max) then
num_unocc = num_unocc + 1
unocc_list(num_unocc) = n
endif
enddo
if (num_occ == 0) then
if (iprint .ge. 2) &
write (stdout, *) "WARNING no occupied bands included in the calculation for kpt=", &
kpt, ", EF[", ifermi, "]=", fermi_energy_list(ifermi), "eV"
cycle
endif
if (num_unocc == 0) then
if (iprint .ge. 2) &
write (stdout, *) "WARNING no unoccupied bands included in the calculation for kpt=", &
kpt, ", EF[", ifermi, "]=", fermi_energy_list(ifermi), "eV"
cycle
endif
allocate (Bnl_orb(num_occ, num_unocc, 3, 3))
call gyrotropic_get_NOA_Bnl_orb(eig, del_eig, AA, num_occ, occ_list, num_unocc, unocc_list, Bnl_orb)
if (present(gyro_NOA_spn)) then
allocate (Bnl_spin(num_occ, num_unocc, 3, 3))
call gyrotropic_get_NOA_Bnl_spin(S_h, num_occ, occ_list, num_unocc, unocc_list, Bnl_spin)
endif
do n1 = 1, num_occ
n = occ_list(n1)
do l1 = 1, num_unocc
l = unocc_list(l1)
wln = eig(l) - eig(n)
multW1(:) = cmplx_1/(wln*wln - gyrotropic_freq_list(:)**2)
multWm(:) = real(multW1)*kweight
multWe(:) = real(-multW1(:)*(2*wln**2*multW1(:) + cmplx_1))*kweight
do ab = 1, 3
a = alpha_A(ab)
b = beta_A(ab)
do c = 1, 3
gyro_NOA_orb(ab, c, ifermi, :) = &
gyro_NOA_orb(ab, c, ifermi, :) + &
multWm(:)*real(AA(l, n, b)*Bnl_orb(n1, l1, a, c) - AA(l, n, a)*Bnl_orb(n1, l1, b, c)) + &
multWe(:)*(del_eig(n, c) + del_eig(l, c))*aimag(AA(n, l, a)*AA(l, n, b))
if (present(gyro_NOA_spn)) &
gyro_NOA_spn(ab, c, ifermi, :) = &
gyro_NOA_spn(ab, c, ifermi, :) + &
multWm(:)*real(AA(l, n, b)*Bnl_spin(n1, l1, a, c) - AA(l, n, a)*Bnl_spin(n1, l1, b, c))
enddo ! c
enddo ! ab
enddo ! l1
enddo ! n1
deallocate (Bnl_orb)
if (present(gyro_NOA_spn)) deallocate (Bnl_spin)
enddo !ifermi
end subroutine gyrotropic_get_NOA_k