!
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=dp), | intent(in) | :: | kpt(3) | |||
| real(kind=dp), | intent(out), | optional | :: | shc_k_fermi(nfermi) | ||
| complex(kind=dp), | intent(out), | optional | :: | shc_k_freq(kubo_nfreq) | ||
| real(kind=dp), | intent(out), | optional | :: | shc_k_band(num_wann) |
subroutine berry_get_shc_klist(kpt, shc_k_fermi, shc_k_freq, shc_k_band)
!====================================================================!
! !
! Contribution from a k-point to the spin Hall conductivity on a list
! of Fermi energies or a list of frequencies or a list of energy bands
! sigma_{alpha,beta}^{gamma}(k), alpha, beta, gamma = 1, 2, 3
! (x, y, z, respectively)
! i.e. the Berry curvature-like term of QZYZ18 Eq.(3) & (4).
! The unit is angstrom^2, similar to that of Berry curvature of AHC.
!
! Note the berry_get_js_k() has not been multiplied by hbar/2 (as
! required by spin operator) and not been divided by hbar (as required
! by the velocity operator). The second velocity operator has not been
! divided by hbar as well. But these two hbar required by velocity
! operators are canceled by the preceding hbar^2 of QZYZ18 Eq.(3).
!
! shc_k_fermi: return a list for different Fermi energies
! shc_k_freq: return a list for different frequencies
! shc_k_band: return a list for each energy band
!
! Junfeng Qiao (18/8/2018) !
!====================================================================!
use w90_constants, only: dp, cmplx_0, cmplx_i
use w90_utility, only: utility_rotate
use w90_parameters, only: num_wann, kubo_eigval_max, kubo_nfreq, &
kubo_freq_list, kubo_adpt_smr, kubo_smr_fixed_en_width, &
kubo_adpt_smr_max, kubo_adpt_smr_fac, berry_kmesh, &
fermi_energy_list, nfermi, shc_alpha, shc_beta, shc_gamma, &
shc_bandshift, shc_bandshift_firstband, shc_bandshift_energyshift
use w90_postw90_common, only: pw90common_get_occ, &
pw90common_fourier_R_to_k_vec, pw90common_kmesh_spacing
use w90_wan_ham, only: wham_get_D_h, wham_get_eig_deleig
use w90_get_oper, only: AA_R
!use w90_comms, only: my_node_id
!!!
! args
real(kind=dp), intent(in) :: kpt(3)
real(kind=dp), optional, intent(out) :: shc_k_fermi(nfermi)
complex(kind=dp), optional, intent(out) :: shc_k_freq(kubo_nfreq)
real(kind=dp), optional, intent(out) :: shc_k_band(num_wann)
! internal vars
logical :: lfreq, lfermi, lband
complex(kind=dp), allocatable :: HH(:, :)
complex(kind=dp), allocatable :: delHH(:, :, :)
complex(kind=dp), allocatable :: UU(:, :)
complex(kind=dp), allocatable :: D_h(:, :, :)
complex(kind=dp), allocatable :: AA(:, :, :)
complex(kind=dp) :: js_k(num_wann, num_wann)
! Adaptive smearing
!
real(kind=dp) :: del_eig(num_wann, 3), joint_level_spacing, &
eta_smr, Delta_k, vdum(3)
integer :: n, m, i, ifreq
real(kind=dp) :: eig(num_wann)
real(kind=dp) :: occ_fermi(num_wann, nfermi), occ_freq(num_wann)
complex(kind=dp) :: omega_list(kubo_nfreq)
real(kind=dp) :: omega, rfac
complex(kind=dp) :: prod, cdum, cfac
allocate (HH(num_wann, num_wann))
allocate (delHH(num_wann, num_wann, 3))
allocate (UU(num_wann, num_wann))
allocate (D_h(num_wann, num_wann, 3))
allocate (AA(num_wann, num_wann, 3))
lfreq = .false.
lfermi = .false.
lband = .false.
if (present(shc_k_freq)) then
shc_k_freq = 0.0_dp
lfreq = .true.
endif
if (present(shc_k_fermi)) then
shc_k_fermi = 0.0_dp
lfermi = .true.
endif
if (present(shc_k_band)) then
shc_k_band = 0.0_dp
lband = .true.
endif
call wham_get_eig_deleig(kpt, eig, del_eig, HH, delHH, UU)
call wham_get_D_h(delHH, UU, eig, D_h)
! Here I apply a scissor operator to the conduction bands, if required in the input
if (shc_bandshift) then
eig(shc_bandshift_firstband:) = eig(shc_bandshift_firstband:) + shc_bandshift_energyshift
end if
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
call berry_get_js_k(kpt, eig, del_eig(:, shc_alpha), &
D_h(:, :, shc_alpha), UU, js_k)
! adpt_smr only works with berry_kmesh, so do not use
! adpt_smr in kpath or kslice plots.
if (kubo_adpt_smr) then
Delta_k = pw90common_kmesh_spacing(berry_kmesh)
endif
if (lfreq) then
call pw90common_get_occ(eig, occ_freq, fermi_energy_list(1))
elseif (lfermi) then
! get occ for different fermi_energy
do i = 1, nfermi
call pw90common_get_occ(eig, occ_fermi(:, i), fermi_energy_list(i))
end do
end if
do n = 1, num_wann
! get Omega_{n,alpha beta}^{gamma}
if (lfreq) then
omega_list = cmplx_0
else if (lfermi .or. lband) then
omega = 0.0_dp
end if
do m = 1, num_wann
if (m == n) cycle
if (eig(m) > kubo_eigval_max .or. eig(n) > kubo_eigval_max) cycle
rfac = eig(m) - eig(n)
!this will calculate AHC
!prod = -rfac*cmplx_i*AA(n, m, shc_alpha) * rfac*cmplx_i*AA(m, n, shc_beta)
prod = js_k(n, m)*cmplx_i*rfac*AA(m, n, shc_beta)
if (kubo_adpt_smr) then
! Eq.(35) YWVS07
vdum(:) = del_eig(m, :) - del_eig(n, :)
joint_level_spacing = sqrt(dot_product(vdum(:), vdum(:)))*Delta_k
eta_smr = min(joint_level_spacing*kubo_adpt_smr_fac, &
kubo_adpt_smr_max)
else
eta_smr = kubo_smr_fixed_en_width
endif
if (lfreq) then
do ifreq = 1, kubo_nfreq
cdum = real(kubo_freq_list(ifreq), dp) + cmplx_i*eta_smr
cfac = -2.0_dp/(rfac**2 - cdum**2)
omega_list(ifreq) = omega_list(ifreq) + cfac*aimag(prod)
end do
else if (lfermi .or. lband) then
rfac = -2.0_dp/(rfac**2 + eta_smr**2)
omega = omega + rfac*aimag(prod)
end if
enddo
if (lfermi) then
do i = 1, nfermi
shc_k_fermi(i) = shc_k_fermi(i) + occ_fermi(n, i)*omega
end do
else if (lfreq) then
shc_k_freq = shc_k_freq + occ_freq(n)*omega_list
else if (lband) then
shc_k_band(n) = omega
end if
enddo
!if (lfermi) then
! write (*, '(3(f9.6,1x),f16.8,1x,1E16.8)') &
! kpt(1), kpt(2), kpt(3), fermi_energy_list(1), shc_k_fermi(1)
!end if
return
contains
!===========================================================!
! PRIVATE PROCEDURES !
!===========================================================!
subroutine berry_get_js_k(kpt, eig, del_alpha_eig, D_alpha_h, UU, js_k)
!====================================================================!
! !
! Contribution from point k to the
! <psi_k | 1/2*(sigma_gamma*v_alpha + v_alpha*sigma_gamma) | psi_k>
!
! QZYZ18 Eq.(23) without hbar/2 (required by spin operator) and
! not divided by hbar (required by velocity operator)
!
! Junfeng Qiao (8/7/2018)
! !
!====================================================================!
use w90_constants, only: dp, cmplx_0, cmplx_i
use w90_utility, only: utility_rotate
use w90_parameters, only: num_wann, shc_alpha, shc_gamma
use w90_postw90_common, only: pw90common_fourier_R_to_k_new, &
pw90common_fourier_R_to_k_vec
use w90_get_oper, only: SS_R, SR_R, SHR_R, SH_R
! args
real(kind=dp), intent(in) :: kpt(3)
real(kind=dp), dimension(:), intent(in) :: eig
real(kind=dp), dimension(:), intent(in) :: del_alpha_eig
complex(kind=dp), dimension(:, :), intent(in) :: D_alpha_h
complex(kind=dp), dimension(:, :), intent(in) :: UU
complex(kind=dp), dimension(:, :), intent(out) :: js_k
! internal vars
complex(kind=dp) :: B_k(num_wann, num_wann)
complex(kind=dp) :: K_k(num_wann, num_wann)
complex(kind=dp) :: L_k(num_wann, num_wann)
complex(kind=dp) :: S_w(num_wann, num_wann)
complex(kind=dp) :: S_k(num_wann, num_wann)
complex(kind=dp) :: SR_w(num_wann, num_wann, 3)
complex(kind=dp) :: SR_alpha_k(num_wann, num_wann)
complex(kind=dp) :: SHR_w(num_wann, num_wann, 3)
complex(kind=dp) :: SHR_alpha_k(num_wann, num_wann)
complex(kind=dp) :: SH_w(num_wann, num_wann, 3)
complex(kind=dp) :: SH_k(num_wann, num_wann)
complex(kind=dp) :: eig_mat(num_wann, num_wann)
complex(kind=dp) :: del_eig_mat(num_wann, num_wann)
!===========
js_k = cmplx_0
!=========== S_k ===========
! < u_k | sigma_gamma | u_k >, QZYZ18 Eq.(25)
! QZYZ18 Eq.(36)
call pw90common_fourier_R_to_k_new(kpt, SS_R(:, :, :, shc_gamma), OO=S_w)
! QZYZ18 Eq.(30)
S_k = utility_rotate(S_w, UU, num_wann)
!=========== K_k ===========
! < u_k | sigma_gamma | \partial_alpha u_k >, QZYZ18 Eq.(26)
! QZYZ18 Eq.(37)
call pw90common_fourier_R_to_k_vec(kpt, SR_R(:, :, :, shc_gamma, :), OO_true=SR_w)
! QZYZ18 Eq.(31)
SR_alpha_k = -cmplx_i*utility_rotate(SR_w(:, :, shc_alpha), UU, num_wann)
K_k = SR_alpha_k + matmul(S_k, D_alpha_h)
!=========== L_k ===========
! < u_k | sigma_gamma.H | \partial_alpha u_k >, QZYZ18 Eq.(27)
! QZYZ18 Eq.(38)
call pw90common_fourier_R_to_k_vec(kpt, SHR_R(:, :, :, shc_gamma, :), OO_true=SHR_w)
! QZYZ18 Eq.(32)
SHR_alpha_k = -cmplx_i*utility_rotate(SHR_w(:, :, shc_alpha), UU, num_wann)
! QZYZ18 Eq.(39)
call pw90common_fourier_R_to_k_vec(kpt, SH_R, OO_true=SH_w)
! QZYZ18 Eq.(32)
SH_k = utility_rotate(SH_w(:, :, shc_gamma), UU, num_wann)
L_k = SHR_alpha_k + matmul(SH_k, D_alpha_h)
!=========== B_k ===========
! < \psi_nk | sigma_gamma v_alpha | \psi_mk >, QZYZ18 Eq.(24)
B_k = cmplx_0
do i = 1, num_wann
eig_mat(i, :) = eig(:)
del_eig_mat(i, :) = del_alpha_eig(:)
end do
! note * is not matmul
B_k = del_eig_mat*S_k + eig_mat*K_k - L_k
!=========== js_k ===========
! QZYZ18 Eq.(23)
! note the S in SR_R,SHR_R,SH_R of get_SHC_R is sigma,
! to get spin current, we need to multiply it by hbar/2,
! also we need to divide it by hbar to recover the velocity
! operator, these are done outside of this subroutine
js_k = 1.0_dp/2.0_dp*(B_k + conjg(transpose(B_k)))
end subroutine berry_get_js_k
end subroutine berry_get_shc_klist