subroutine gyrotropic_get_k_list(ws_region, num_valence_bands, have_disentangled, 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, num_wann, print_output, &
fermi_energy_list, wannier_data, eigval, real_lattice, mp_grid, &
num_bands, num_kpts, u_matrix, v_matrix, dis_manifold, kpt_latt, &
pw90_gyrotropic, scissors_shift, effective_model, pw90_band_deriv_degen, &
ws_distance, wigner_seitz, stdout, seedname, timer, error, &
comm, HH_R, AA_R, BB_R, CC_R, SS_R)
!======================================================================!
! !
! 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_comms, only: w90_comm_type, mpirank
use w90_constants, only: dp, cmplx_0, cmplx_i
use w90_postw90_types, only: pw90_gyrotropic_type, pw90_band_deriv_degen_type, wigner_seitz_type
use w90_types, only: dis_manifold_type, print_output_type, &
wannier_data_type, ws_region_type, ws_distance_type, timer_list_type
use w90_postw90_common, only: pw90common_fourier_R_to_k_new_second_d, &
pw90common_fourier_R_to_k_vec
use w90_spin, only: spin_get_S
use w90_utility, only: utility_diagonalize, utility_rotate, utility_rotate_diag, &
utility_w0gauss, utility_recip_lattice_base
use w90_wan_ham, only: wham_get_eig_deleig, wham_get_D_h
implicit none
! arguments
type(dis_manifold_type), intent(in) :: dis_manifold
real(kind=dp), allocatable, intent(in) :: fermi_energy_list(:)
type(pw90_gyrotropic_type), intent(in) :: pw90_gyrotropic
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(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) :: mp_grid(3)
integer, intent(in) :: num_bands, num_kpts, num_wann, num_valence_bands
integer, intent(in) :: stdout
real(kind=dp), allocatable, intent(inout) :: gyro_DOS(:)
real(kind=dp), allocatable, intent(inout) :: gyro_Dw(:, :, :, :)
real(kind=dp), allocatable, intent(inout) :: gyro_NOA_spn(:, :, :, :)
real(kind=dp), allocatable, intent(inout) :: gyro_NOA_orb(:, :, :, :)
real(kind=dp), allocatable, intent(inout) :: gyro_K_spn(:, :, :)
real(kind=dp), allocatable, intent(inout) :: gyro_K_orb(:, :, :)
real(kind=dp), allocatable, intent(inout) :: gyro_D(:, :, :)
real(kind=dp), allocatable, intent(inout) :: gyro_C(:, :, :)
real(kind=dp), intent(in) :: eigval(:, :)
real(kind=dp), intent(in) :: kpt(3), kweight
real(kind=dp), intent(in) :: real_lattice(3, 3)
real(kind=dp), intent(in) :: scissors_shift
complex(kind=dp), allocatable, intent(inout) :: AA_R(:, :, :, :)
complex(kind=dp), allocatable, intent(inout) :: BB_R(:, :, :, :)
complex(kind=dp), allocatable, intent(inout) :: CC_R(:, :, :, :, :)
complex(kind=dp), allocatable, intent(inout) :: HH_R(:, :, :)
complex(kind=dp), allocatable, intent(inout) :: SS_R(:, :, :, :)
complex(kind=dp), intent(in) :: u_matrix(:, :, :), v_matrix(:, :, :)
character(len=50), intent(in) :: seedname
logical, intent(in) :: eval_K, eval_D, eval_Dw, eval_C, eval_NOA, eval_spn, eval_dos
logical, intent(in) :: have_disentangled
logical, intent(in) :: effective_model
! local variables
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, ifermi, fermi_n
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
if (pw90_gyrotropic%smearing%use_adaptive) then
call set_error_input(error, 'Adaptive smearing not allowed in Gyrotropic', comm)
return
end if
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(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
if (eval_Dw .or. eval_NOA) then
allocate (AA(num_wann, num_wann, 3))
call wham_get_D_h(delHH, D_h, UU, eig, num_wann)
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
end if
if (eval_Dw) allocate (curv_w_nk(num_wann, pw90_gyrotropic%nfreq, 3))
eta_smr = pw90_gyrotropic%smearing%fixed_width
got_spin = .false.
do n1 = 1, pw90_gyrotropic%num_bands
n = pw90_gyrotropic%band_list(n1)
!
! ***ADJUSTABLE PARAMETER***
! avoid degeneracies
!---------------------------------------------------
if (n > 1) then
if (eig(n) - eig(n - 1) <= pw90_gyrotropic%degen_thresh) cycle
end if
if (n < num_wann) then
if (eig(n + 1) - eig(n) <= pw90_gyrotropic%degen_thresh) cycle
end if
!---------------------------------------------------
fermi_n = size(fermi_energy_list)
got_orb_n = .false.
do ifermi = 1, fermi_n
arg = (eig(n) - fermi_energy_list(ifermi))/eta_smr
!
! To save time: far from the Fermi surface, negligible contribution
!
!-------------------------
if (abs(arg) > pw90_gyrotropic%smearing%max_arg) cycle
!-------------------------
!
! Spin is computed for all bands simultaneously
!
if (eval_spn .and. .not. got_spin) then
call spin_get_S(kpt, S, num_wann, ws_region, wannier_data, real_lattice, mp_grid, &
ws_distance, HH_R, SS_R, wigner_seitz, error, comm)
if (allocated(error)) return
got_spin = .true. ! Do it for only one value of ifermi and n
end if
! 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(dis_manifold, fermi_energy_list, kpt_latt, ws_region, &
print_output, wannier_data, ws_distance, wigner_seitz, &
AA_R, BB_R, CC_R, HH_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, &
imf_k, img_k, imh_k, occ)
if (allocated(error)) return
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))
end do
else if (eval_D) then
occ = 0.0_dp
occ(n) = 1.0_dp
call berry_get_imf_klist(dis_manifold, fermi_energy_list, kpt_latt, ws_region, &
print_output, wannier_data, ws_distance, wigner_seitz, AA_R, &
BB_R, CC_R, HH_R, u_matrix, v_matrix, eigval, kpt, &
real_lattice, imf_k, scissors_shift, mp_grid, num_bands, &
num_kpts, num_wann, num_valence_bands, effective_model, &
have_disentangled, seedname, stdout, timer, error, comm, occ)
if (allocated(error)) return
do i = 1, 3
curv_nk(i) = sum(imf_k(:, i, 1))
end do
got_orb_n = .true. ! Do it for only one value of ifermi
end if
if (eval_Dw) call gyrotropic_get_curv_w_k(eig, AA, curv_w_nk, pw90_gyrotropic)
got_orb_n = .true. ! Do it for only one value of ifermi
end if
!
delta = utility_w0gauss(arg, pw90_gyrotropic%smearing%type_index, error, comm) &
/eta_smr*kweight ! Broadened delta(E_nk-E_f)
if (allocated(error)) return
!
! 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)
end do
end if
if (eval_C) gyro_C(:, j, ifermi) = &
gyro_C(:, j, ifermi) + del_eig(n, :)*del_eig(n, j)*delta
end do !j
if (eval_dos) gyro_DOS(ifermi) = gyro_DOS(ifermi) + delta
end do !ifermi
end do !n
if (eval_NOA) then
if (eval_spn) then
call 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)
if (allocated(error)) return
else
call 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)
if (allocated(error)) return
end if
end if
end subroutine gyrotropic_get_k_list
