subroutine berry_get_shc_klist(pw90_berry, dis_manifold, fermi_energy_list, kpt_latt, &
pw90_band_deriv_degen, ws_region, pw90_spin_hall, print_output, &
wannier_data, ws_distance, wigner_seitz, AA_R, HH_R, SH_R, SHR_R, &
SR_R, SS_R, SAA_R, SBB_R, u_matrix, v_matrix, eigval, kpt, &
real_lattice, scissors_shift, mp_grid, fermi_n, num_bands, &
num_kpts, num_wann, num_valence_bands, effective_model, &
have_disentangled, seedname, stdout, timer, error, comm, &
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, utility_recip_lattice_base
use w90_comms, only: w90_comm_type
use w90_types, only: print_output_type, wannier_data_type, &
dis_manifold_type, kmesh_info_type, ws_region_type, ws_distance_type, timer_list_type
use w90_postw90_types, only: pw90_berry_mod_type, pw90_spin_hall_type, &
pw90_band_deriv_degen_type, wigner_seitz_type
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
implicit none
! arguments
type(pw90_berry_mod_type), intent(in) :: pw90_berry
type(dis_manifold_type), intent(in) :: dis_manifold
real(kind=dp), allocatable, intent(in) :: fermi_energy_list(:)
real(kind=dp), intent(in) :: kpt_latt(:, :)
type(pw90_band_deriv_degen_type), intent(in) :: pw90_band_deriv_degen
type(print_output_type), intent(in) :: print_output
type(ws_region_type), intent(in) :: ws_region
type(pw90_spin_hall_type), intent(in) :: pw90_spin_hall
type(w90_comm_type), intent(in) :: comm
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(timer_list_type), intent(inout) :: timer
type(w90_error_type), allocatable, intent(out) :: error
integer, intent(in) :: num_wann, num_bands, num_kpts, num_valence_bands, fermi_n
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: stdout
real(kind=dp), intent(in) :: kpt(3)
real(kind=dp), intent(in) :: eigval(:, :)
real(kind=dp), intent(in) :: real_lattice(3, 3)
real(kind=dp), intent(in) :: scissors_shift
complex(kind=dp), intent(in) :: u_matrix(:, :, :), v_matrix(:, :, :)
complex(kind=dp), allocatable, intent(inout) :: AA_R(:, :, :, :) ! <0n|r|Rm>
complex(kind=dp), allocatable, intent(inout) :: HH_R(:, :, :) ! <0n|r|Rm>
complex(kind=dp), allocatable, intent(inout) :: SR_R(:, :, :, :, :) ! <0n|sigma_x,y,z.(r-R)_alpha|Rm>
complex(kind=dp), allocatable, intent(inout) :: SHR_R(:, :, :, :, :) ! <0n|sigma_x,y,z.H.(r-R)_alpha|Rm>
complex(kind=dp), allocatable, intent(inout) :: SH_R(:, :, :, :) ! <0n|sigma_x,y,z.H|Rm>
complex(kind=dp), allocatable, intent(inout) :: SS_R(:, :, :, :) ! <0n|sigma_x,y,z|Rm>
complex(kind=dp), allocatable, intent(inout) :: SAA_R(:, :, :, :, :)
complex(kind=dp), allocatable, intent(inout) :: SBB_R(:, :, :, :, :)
character(len=50), intent(in) :: seedname
logical, intent(in) :: have_disentangled
logical, intent(in) :: effective_model
complex(kind=dp), optional, intent(out) :: shc_k_freq(pw90_berry%kubo_nfreq)
real(kind=dp), optional, intent(out) :: shc_k_fermi(fermi_n)
real(kind=dp), optional, intent(out) :: shc_k_band(num_wann)
! internal vars
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)
logical :: lfreq, lfermi, lband
real(kind=dp) :: recip_lattice(3, 3), volume
integer :: n, m, i, ifreq
! Adaptive smearing
real(kind=dp) :: del_eig(num_wann, 3), joint_level_spacing, eta_smr, Delta_k, vdum(3)
real(kind=dp) :: eig(num_wann)
real(kind=dp) :: occ_fermi(num_wann, fermi_n), occ_freq(num_wann)
real(kind=dp) :: omega, rfac
complex(kind=dp) :: omega_list(pw90_berry%kubo_nfreq)
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.
end if
if (present(shc_k_fermi)) then
shc_k_fermi = 0.0_dp
lfermi = .true.
end if
if (present(shc_k_band)) then
shc_k_band = 0.0_dp
lband = .true.
end if
call wham_get_eig_deleig(dis_manifold, kpt_latt, pw90_band_deriv_degen, ws_region, &
print_output, wannier_data, ws_distance, wigner_seitz, delHH, HH, &
HH_R, u_matrix, UU, v_matrix, del_eig, eig, eigval, kpt, &
real_lattice, scissors_shift, mp_grid, num_bands, num_kpts, num_wann, &
num_valence_bands, effective_model, have_disentangled, seedname, &
stdout, timer, error, comm)
if (allocated(error)) return
call wham_get_D_h(delHH, D_h, UU, eig, num_wann)
! Here I apply a scissor operator to the conduction bands, if required in the input
if (pw90_spin_hall%bandshift) then
eig(pw90_spin_hall%bandshift_firstband:) = eig(pw90_spin_hall%bandshift_firstband:) + pw90_spin_hall%bandshift_energyshift
end if
call pw90common_fourier_R_to_k_vec(ws_region, wannier_data, ws_distance, wigner_seitz, AA_R, &
kpt, real_lattice, mp_grid, num_wann, error, comm, &
OO_true=AA)
if (allocated(error)) return
do i = 1, 3
AA(:, :, i) = utility_rotate(AA(:, :, i), UU, num_wann)
end do
AA = AA + cmplx_i*D_h ! Eq.(25) WYSV06
call berry_get_js_k(ws_region, pw90_spin_hall, wannier_data, ws_distance, wigner_seitz, &
D_h(:, :, pw90_spin_hall%alpha), js_k, SH_R, SHR_R, SR_R, SS_R, SAA_R, &
SBB_R, UU, eig, del_eig(:, pw90_spin_hall%alpha), &
delHH(:, :, pw90_spin_hall%alpha), kpt, real_lattice, mp_grid, num_wann)
! adpt_smr only works with pw90_pw90_berry%kmesh%mesh, so do not use
! adpt_smr in kpath or kslice plots.
if (pw90_berry%kubo_smearing%use_adaptive) then
call utility_recip_lattice_base(real_lattice, recip_lattice, volume)
Delta_k = pw90common_kmesh_spacing(pw90_berry%kmesh%mesh, recip_lattice)
end if
if (lfreq) then
call pw90common_get_occ(fermi_energy_list(1), eig, occ_freq, num_wann)
elseif (lfermi) then
! get occ for different fermi_energy
do i = 1, fermi_n
call pw90common_get_occ(fermi_energy_list(i), eig, occ_fermi(:, i), num_wann)
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) > pw90_berry%kubo_eigval_max .or. eig(n) > pw90_berry%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, pw90_spin_hall%beta)
if (pw90_berry%kubo_smearing%use_adaptive) 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*pw90_berry%kubo_smearing%adaptive_prefactor, &
pw90_berry%kubo_smearing%adaptive_max_width)
else
eta_smr = pw90_berry%kubo_smearing%fixed_width
end if
if (lfreq) then
do ifreq = 1, pw90_berry%kubo_nfreq
cdum = real(pw90_berry%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
end do
if (lfermi) then
do i = 1, fermi_n
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
end do
!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(ws_region, pw90_spin_hall, wannier_data, ws_distance, wigner_seitz, &
D_alpha_h, js_k, SH_R, SHR_R, SR_R, SS_R, SAA_R, SBB_R, UU, eig, &
del_alpha_eig, delHH_alpha, kpt, real_lattice, mp_grid, num_wann)
!================================================!
!
! 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_types, only: print_output_type, wannier_data_type, ws_region_type, &
ws_distance_type
use w90_postw90_types, only: pw90_spin_hall_type, wigner_seitz_type
use w90_postw90_common, only: pw90common_fourier_R_to_k_new, pw90common_fourier_R_to_k_vec
implicit none
! arguments
type(ws_region_type), intent(in) :: ws_region
type(pw90_spin_hall_type), intent(in) :: pw90_spin_hall
type(wannier_data_type), intent(in) :: wannier_data
type(wigner_seitz_type), intent(in) :: wigner_seitz
type(ws_distance_type), intent(inout) :: ws_distance
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: num_wann
real(kind=dp), intent(in) :: kpt(3)
real(kind=dp), intent(in) :: eig(:)
real(kind=dp), intent(in) :: del_alpha_eig(:)
real(kind=dp), intent(in) :: real_lattice(3, 3)
complex(kind=dp), dimension(:, :), intent(in) :: delHH_alpha
complex(kind=dp), intent(in) :: D_alpha_h(:, :)
complex(kind=dp), intent(in) :: UU(:, :)
complex(kind=dp), intent(out) :: js_k(:, :)
complex(kind=dp), allocatable, intent(inout) :: SR_R(:, :, :, :, :) ! <0n|sigma_x,y,z.(r-R)_alpha|Rm>
complex(kind=dp), allocatable, intent(inout) :: SHR_R(:, :, :, :, :) ! <0n|sigma_x,y,z.H.(r-R)_alpha|Rm>
complex(kind=dp), allocatable, intent(inout) :: SH_R(:, :, :, :) ! <0n|sigma_x,y,z.H|Rm>
complex(kind=dp), allocatable, intent(inout) :: SS_R(:, :, :, :) ! <0n|sigma_x,y,z|Rm>
complex(kind=dp), allocatable, intent(inout) :: SAA_R(:, :, :, :, :)
complex(kind=dp), allocatable, intent(inout) :: SBB_R(:, :, :, :, :)
! 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)
!ryoo
complex(kind=dp) :: SAA(num_wann, num_wann, 3, 3)
complex(kind=dp) :: SBB(num_wann, num_wann, 3, 3)
complex(kind=dp) :: VV0(num_wann, num_wann)
complex(kind=dp) :: spinvel0(num_wann, num_wann)
integer :: i
!================================================
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(ws_region, wannier_data, ws_distance, wigner_seitz, &
SS_R(:, :, :, pw90_spin_hall%gamma), kpt, real_lattice, &
mp_grid, num_wann, error, comm, OO=S_w)
if (allocated(error)) return
! QZYZ18 Eq.(30)
S_k = utility_rotate(S_w, UU, num_wann)
if (index(pw90_spin_hall%method, 'qiao') > 0) then !if Qiao
!================================================ K_k ===========
! < u_k | sigma_gamma | \partial_alpha u_k >, QZYZ18 Eq.(26)
! QZYZ18 Eq.(37)
call pw90common_fourier_R_to_k_vec(ws_region, wannier_data, ws_distance, wigner_seitz, &
SR_R(:, :, :, pw90_spin_hall%gamma, :), kpt, &
real_lattice, mp_grid, num_wann, error, comm, &
OO_true=SR_w)
if (allocated(error)) return
! QZYZ18 Eq.(31)
SR_alpha_k = -cmplx_i*utility_rotate(SR_w(:, :, pw90_spin_hall%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(ws_region, wannier_data, ws_distance, wigner_seitz, &
SHR_R(:, :, :, pw90_spin_hall%gamma, :), kpt, &
real_lattice, mp_grid, num_wann, error, comm, &
OO_true=SHR_w)
if (allocated(error)) return
! QZYZ18 Eq.(32)
SHR_alpha_k = -cmplx_i*utility_rotate(SHR_w(:, :, pw90_spin_hall%alpha), UU, num_wann)
! QZYZ18 Eq.(39)
call pw90common_fourier_R_to_k_vec(ws_region, wannier_data, ws_distance, wigner_seitz, &
SH_R, kpt, real_lattice, mp_grid, num_wann, error, &
comm, OO_true=SH_w)
if (allocated(error)) return
! QZYZ18 Eq.(32)
SH_k = utility_rotate(SH_w(:, :, pw90_spin_hall%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)))
else !if Ryoo (PRB RPS19 Eq.(21))
!RPS19 Eqs.(37)-(40)
call pw90common_fourier_R_to_k_new(ws_region, wannier_data, ws_distance, wigner_seitz, &
SAA_R(:, :, :, pw90_spin_hall%gamma, &
pw90_spin_hall%alpha), kpt, real_lattice, mp_grid, &
num_wann, error, comm, &
OO=SAA(:, :, pw90_spin_hall%gamma, pw90_spin_hall%alpha))
if (allocated(error)) return
call pw90common_fourier_R_to_k_new(ws_region, wannier_data, ws_distance, wigner_seitz, &
SBB_R(:, :, :, pw90_spin_hall%gamma, &
pw90_spin_hall%alpha), kpt, real_lattice, mp_grid, &
num_wann, error, comm, &
OO=SBB(:, :, pw90_spin_hall%gamma, pw90_spin_hall%alpha))
if (allocated(error)) return
call pw90common_fourier_R_to_k_new(ws_region, wannier_data, ws_distance, wigner_seitz, &
HH_R, kpt, real_lattice, mp_grid, num_wann, error, &
comm, OO=HH, OO_dx=delHH(:, :, 1), &
OO_dy=delHH(:, :, 2), OO_dz=delHH(:, :, 3))
if (allocated(error)) return
VV0(:, :) = utility_rotate(delHH_alpha(:, :), UU, num_wann)
SAA(:, :, pw90_spin_hall%gamma, pw90_spin_hall%alpha) = &
utility_rotate(SAA(:, :, pw90_spin_hall%gamma, pw90_spin_hall%alpha), UU, num_wann)
SBB(:, :, pw90_spin_hall%gamma, pw90_spin_hall%alpha) = &
utility_rotate(SBB(:, :, pw90_spin_hall%gamma, pw90_spin_hall%alpha), UU, num_wann)
spinVel0(:, :) = matmul(VV0(:, :), S_k(:, :)) + &
matmul(S_k(:, :), VV0(:, :))
do n = 1, num_wann
do m = 1, num_wann !RPS19 Eq.(21) and Eq.(26)
js_k(n, m) = spinVel0(n, m) &
- cmplx_i*(eig(m)*SAA(n, m, pw90_spin_hall%gamma, pw90_spin_hall%alpha) &
- SBB(n, m, pw90_spin_hall%gamma, pw90_spin_hall%alpha))
js_k(n, m) = js_k(n, m) &
+ cmplx_i*(eig(n)*conjg(SAA(m, n, pw90_spin_hall%gamma, pw90_spin_hall%alpha)) &
- conjg(SBB(m, n, pw90_spin_hall%gamma, pw90_spin_hall%alpha)))
end do
end do
js_k = js_k/2.0_dp
end if
end subroutine berry_get_js_k
end subroutine berry_get_shc_klist
