subroutine gyrotropic_get_k_list(kpt, kweight, &
gyro_K_spn, gyro_K_orb, gyro_D, gyro_Dw, gyro_C, &
gyro_DOS, gyro_NOA_orb, gyro_NOA_spn, &
eval_K, eval_D, eval_Dw, eval_NOA, eval_spn, eval_C, eval_dos)
!======================================================================!
! !
! Contribution from point k to the GME tensor, Eq.(9) of ZMS16, !
! evaluated in the clean limit of omega.tau >> 1 where it is real. !
! The following two quantities are calculated (sigma = Pauli matrix): !
! !
! gyro_K_spn_k = delta(E_kn-E_f).(d E_{kn}/d k_i).sigma_{kn,j} !
! [units of length] !
! !
! gyro_K_orb_k = delta(E_kn-E_f).(d E_{kn}/d k_i).(2.hbar/e).m^orb_{kn,j} !
! [units of (length^3)*energy] !
! !
! gyro_D_k = delta(E_kn-E_f).(d E_{kn}/d k_i).Omega_{kn,j} !
! [units of length^3] !
! !
! gyro_Dw_k = delta(E_kn-E_f).(d E_{kn}/d k_i).tildeOmega_{kn,j} !
! [units of length^3] !
! !
! gyro_C_k = delta(E_kn-E_f).(d E_{kn}/d k_i).(d E_{kn}/d k_j) !
! [units of energy*length^3] !
! !
! gyro_DOS_k = delta(E_kn-E_f) !
! [units of 1/Energy] !
! !
! gme_NOA_orb_k = ????? !
! !
! gme_NOA_spn_k = ?????? !
! !
!======================================================================!
use w90_constants, only: dp, cmplx_0, cmplx_i
use w90_utility, only: utility_rotate, utility_rotate_diag, utility_w0gauss
use w90_parameters, only: num_wann, fermi_energy_list, &
gyrotropic_smr_index, nfermi, gyrotropic_nfreq, &
gyrotropic_degen_thresh, gyrotropic_smr_max_arg, &
gyrotropic_band_list, gyrotropic_num_bands, &
gyrotropic_smr_fixed_en_width
use w90_postw90_common, only: pw90common_get_occ, &
pw90common_fourier_R_to_k_vec
use w90_wan_ham, only: wham_get_eig_deleig, wham_get_D_h
use w90_get_oper, only: HH_R, SS_R, AA_R
use w90_spin, only: spin_get_S
use w90_io, only: stdout
! Arguments
!
real(kind=dp), intent(in) :: kpt(3), kweight
real(kind=dp), dimension(:, :, :), intent(inout) :: gyro_K_spn, &
gyro_K_orb, &
gyro_D, &
gyro_C
real(kind=dp), dimension(:, :, :, :), intent(inout) :: gyro_Dw, gyro_NOA_spn, gyro_NOA_orb
real(kind=dp), dimension(:), intent(inout) :: gyro_DOS
logical, intent(in) :: eval_K, eval_D, eval_Dw, &
eval_C, eval_NOA, eval_spn, eval_dos
complex(kind=dp), allocatable :: UU(:, :)
complex(kind=dp), allocatable :: HH(:, :)
complex(kind=dp), allocatable :: delHH(:, :, :)
complex(kind=dp), allocatable :: SS(:, :, :)
complex(kind=dp), allocatable :: AA(:, :, :)
complex(kind=dp), allocatable :: D_h(:, :, :)
real(kind=dp), allocatable :: curv_w_nk(:, :, :)
integer :: i, j, n, n1, m1, m, ifermi
real(kind=dp) :: delta, occ(num_wann), &
eig(num_wann), del_eig(num_wann, 3), &
S(num_wann, 3), eta_smr, arg, &
orb_nk(3), curv_nk(3), &
imf_k(3, 3, 1), img_k(3, 3, 1), imh_k(3, 3, 1)
logical :: got_spin, got_orb_n
allocate (UU(num_wann, num_wann))
allocate (HH(num_wann, num_wann))
allocate (delHH(num_wann, num_wann, 3))
allocate (D_h(num_wann, num_wann, 3))
if (eval_spn) allocate (SS(num_wann, num_wann, 3))
call wham_get_eig_deleig(kpt, eig, del_eig, HH, delHH, UU)
if (eval_Dw .or. eval_NOA) then
allocate (AA(num_wann, num_wann, 3))
call wham_get_D_h(delHH, UU, eig, D_h)
call pw90common_fourier_R_to_k_vec(kpt, AA_R, OO_true=AA)
do i = 1, 3
AA(:, :, i) = utility_rotate(AA(:, :, i), UU, num_wann)
enddo
AA = AA + cmplx_i*D_h ! Eq.(25) WYSV06
endif
if (eval_Dw) allocate (curv_w_nk(num_wann, gyrotropic_nfreq, 3))
eta_smr = gyrotropic_smr_fixed_en_width
got_spin = .false.
do n1 = 1, gyrotropic_num_bands
n = gyrotropic_band_list(n1)
!
! ***ADJUSTABLE PARAMETER***
! avoid degeneracies
!---------------------------------------------------
if (n > 1) then
if (eig(n) - eig(n - 1) <= gyrotropic_degen_thresh) cycle
endif
if (n < num_wann) then
if (eig(n + 1) - eig(n) <= gyrotropic_degen_thresh) cycle
endif
!---------------------------------------------------
got_orb_n = .false.
do ifermi = 1, nfermi
arg = (eig(n) - fermi_energy_list(ifermi))/eta_smr
!
! To save time: far from the Fermi surface, negligible contribution
!
!-------------------------
if (abs(arg) > gyrotropic_smr_max_arg) cycle
!-------------------------
!
! Spin is computed for all bands simultaneously
!
if (eval_spn .and. .not. got_spin) then
call spin_get_S(kpt, S)
got_spin = .true. ! Do it for only one value of ifermi and n
endif
! Orbital quantities are computed for each band separately
if (.not. got_orb_n) then
if (eval_K) then
! Fake occupations: band n occupied, others empty
occ = 0.0_dp
occ(n) = 1.0_dp
call berry_get_imfgh_klist(kpt, imf_k, img_k, imh_k, occ)
do i = 1, 3
orb_nk(i) = sum(imh_k(:, i, 1)) - sum(img_k(:, i, 1))
curv_nk(i) = sum(imf_k(:, i, 1))
enddo
else if (eval_D) then
occ = 0.0_dp
occ(n) = 1.0_dp
call berry_get_imf_klist(kpt, imf_k, occ)
do i = 1, 3
curv_nk(i) = sum(imf_k(:, i, 1))
enddo
got_orb_n = .true. ! Do it for only one value of ifermi
endif
if (eval_Dw) call gyrotropic_get_curv_w_k(eig, AA, curv_w_nk)
got_orb_n = .true. ! Do it for only one value of ifermi
endif
!
delta = utility_w0gauss(arg, gyrotropic_smr_index)/eta_smr*kweight ! Broadened delta(E_nk-E_f)
!
! Loop over Cartesian tensor components
!
do j = 1, 3
if (eval_K .and. eval_spn) gyro_K_spn(:, j, ifermi) = &
gyro_K_spn(:, j, ifermi) + del_eig(n, :)*S(n, j)*delta
if (eval_K) gyro_K_orb(:, j, ifermi) = &
gyro_K_orb(:, j, ifermi) + del_eig(n, :)*orb_nk(j)*delta
if (eval_D) gyro_D(:, j, ifermi) = &
gyro_D(:, j, ifermi) + del_eig(n, :)*curv_nk(j)*delta
if (eval_Dw) then
do i = 1, 3
gyro_Dw(i, j, ifermi, :) = &
gyro_Dw(i, j, ifermi, :) + del_eig(n, i)*delta*curv_w_nk(n, :, j)
enddo
endif
if (eval_C) gyro_C(:, j, ifermi) = &
gyro_C(:, j, ifermi) + del_eig(n, :)*del_eig(n, j)*delta
enddo !j
if (eval_dos) gyro_DOS(ifermi) = gyro_DOS(ifermi) + delta
enddo !ifermi
enddo !n
if (eval_NOA) then
if (eval_spn) then
call gyrotropic_get_NOA_k(kpt, kweight, eig, del_eig, AA, UU, gyro_NOA_orb, gyro_NOA_spn)
else
call gyrotropic_get_NOA_k(kpt, kweight, eig, del_eig, AA, UU, gyro_NOA_orb)
endif
endif
end subroutine gyrotropic_get_k_list