subroutine gyrotropic_get_NOA_k(ws_region, kpt, kweight, eig, del_eig, AA, UU, gyro_NOA_orb, &
num_wann, print_output, fermi_energy_list, wannier_data, &
real_lattice, mp_grid, pw90_gyrotropic, ws_distance, &
wigner_seitz, stdout, error, comm, SS_R, 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_postw90_types, only: pw90_gyrotropic_type, wigner_seitz_type
use w90_constants, only: dp, cmplx_1
use w90_io, only: io_time
use w90_types, only: print_output_type, wannier_data_type, ws_region_type, &
ws_distance_type
use w90_postw90_common, only: pw90common_fourier_R_to_k_new
use w90_spin, only: spin_get_S
use w90_utility, only: utility_rotate
use w90_comms, only: w90_comm_type
implicit none
! arguments
real(kind=dp), allocatable, intent(in) :: fermi_energy_list(:)
type(pw90_gyrotropic_type), intent(in) :: pw90_gyrotropic
type(print_output_type), intent(in) :: print_output
type(ws_region_type), intent(in) :: ws_region
type(wannier_data_type), intent(in) :: wannier_data
type(wigner_seitz_type), intent(inout) :: wigner_seitz
type(ws_distance_type), intent(inout) :: ws_distance
type(w90_error_type), allocatable, intent(out) :: error
type(w90_comm_type), intent(in) :: comm
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: num_wann
integer, intent(in) :: stdout
real(kind=dp), intent(in) :: del_eig(:, :)
real(kind=dp), intent(in) :: eig(:)
real(kind=dp), intent(inout) :: gyro_NOA_orb(:, :, :, :)
real(kind=dp), intent(inout), optional :: gyro_NOA_spn(:, :, :, :)
real(kind=dp), intent(in) :: kpt(3), kweight
real(kind=dp), intent(in) :: real_lattice(3, 3)
complex(kind=dp), allocatable, intent(inout) :: SS_R(:, :, :, :)
complex(kind=dp), intent(in) :: AA(:, :, :)
complex(kind=dp), intent(in) :: UU(:, :)
! local variables
integer :: j, n, l, n1, l1, a, b, c, ab, ifermi
integer :: num_occ, num_unocc, occ_list(num_wann), unocc_list(num_wann)
real(kind=dp) :: wln
real(kind=dp) :: multWe(pw90_gyrotropic%nfreq), multWm(pw90_gyrotropic%nfreq)
complex(kind=dp) :: multW1(pw90_gyrotropic%nfreq)
complex(kind=dp), allocatable :: S_h(:, :, :)
complex(kind=dp), allocatable :: SS(:, :, :)
complex(kind=dp), allocatable :: Bnl_orb(:, :, :, :)
complex(kind=dp), allocatable :: Bnl_spin(:, :, :, :)
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(ws_region, wannier_data, ws_distance, wigner_seitz, &
SS_R(:, :, :, j), kpt, real_lattice, mp_grid, num_wann, &
error, comm, OO=SS(:, :, j))
if (allocated(error)) return
S_h(:, :, j) = utility_rotate(SS(:, :, j), UU, num_wann)
end do
end if
do ifermi = 1, size(fermi_energy_list)
num_occ = 0
num_unocc = 0
do n1 = 1, pw90_gyrotropic%num_bands
n = pw90_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) < pw90_gyrotropic%eigval_max) then
num_unocc = num_unocc + 1
unocc_list(num_unocc) = n
end if
end do
if (num_occ == 0) then
if (print_output%iprint .ge. 2) &
write (stdout, *) "WARNING no occupied bands included in the calculation for kpt=", &
kpt, ", EF[", ifermi, "]=", fermi_energy_list(ifermi), "eV"
cycle
end if
if (num_unocc == 0) then
if (print_output%iprint .ge. 2) &
write (stdout, *) "WARNING no unoccupied bands included in the calculation for kpt=", &
kpt, ", EF[", ifermi, "]=", fermi_energy_list(ifermi), "eV"
cycle
end if
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, pw90_gyrotropic)
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)
end if
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 - pw90_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))
end do ! c
end do ! ab
end do ! l1
end do ! n1
deallocate (Bnl_orb)
if (present(gyro_NOA_spn)) deallocate (Bnl_spin)
end do !ifermi
end subroutine gyrotropic_get_NOA_k
