subroutine berry_get_sc_klist(pw90_berry, dis_manifold, fermi_energy_list, kmesh_info, kpt_latt, &
ws_region, print_output, pw90_band_deriv_degen, wannier_data, &
ws_distance, wigner_seitz, AA_R, HH_R, u_matrix, v_matrix, eigval, &
kpt, real_lattice, sc_k_list, scissors_shift, mp_grid, num_bands, &
num_kpts, num_wann, num_valence_bands, effective_model, &
have_disentangled, seedname, stdout, timer, error, comm)
!================================================!
!
! Contribution from point k to the nonlinear shift current
! [integrand of Eq.8 IATS18]
! Notation correspondence with IATS18:
! AA_da_bar <--> \mathbbm{b}
! AA_bar <--> \mathbbm{a}
! HH_da_bar <--> \mathbbm{v}
! HH_dadb_bar <--> \mathbbm{w}
! D_h(n,m) <--> \mathbbm{v}_{nm} * Re[1/(E_{m}-E_{n}+i*sc_eta)]
! D_h_no_eta(n,m) <--> \mathbbm{v}_{nm} / (E_{m}-E_{n})
! sum_AD <--> summatory of Eq. 32 IATS18
! sum_HD <--> summatory of Eq. 30 IATS18
! eig_da(n)-eig_da(m) <--> \mathbbm{Delta}_{nm}
!
!================================================!
use w90_constants, only: dp, cmplx_0, cmplx_i
use w90_utility, only: utility_re_tr, utility_im_tr, utility_w0gauss, utility_w0gauss_vec
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_band_deriv_degen_type, wigner_seitz_type
use w90_postw90_common, only: pw90common_fourier_R_to_k_vec_dadb, &
pw90common_fourier_R_to_k_new_second_d, pw90common_get_occ, &
pw90common_kmesh_spacing, pw90common_fourier_R_to_k_vec_dadb_TB_conv
use w90_wan_ham, only: wham_get_D_h, &
wham_get_eig_UU_HH_AA_sc, wham_get_eig_deleig, wham_get_D_h_P_value, &
wham_get_eig_deleig_TB_conv, wham_get_eig_UU_HH_AA_sc_TB_conv
use w90_comms, only: w90_comm_type
use w90_utility, only: utility_rotate, utility_zdotu, utility_recip_lattice_base
implicit none
! arguments
type(pw90_berry_mod_type), intent(in) :: pw90_berry
type(dis_manifold_type), intent(in) :: dis_manifold
type(kmesh_info_type), intent(in) :: kmesh_info
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) :: num_wann, num_bands, num_kpts, num_valence_bands
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: stdout
real(kind=dp), intent(in) :: kpt(3)
real(kind=dp), intent(out) :: sc_k_list(:, :, :)
real(kind=dp), intent(in) :: eigval(:, :)
real(kind=dp), intent(in) :: real_lattice(3, 3)
real(kind=dp), intent(in) :: scissors_shift
real(kind=dp), allocatable, intent(in) :: fermi_energy_list(:)
real(kind=dp), intent(in) :: kpt_latt(:, :)
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>
character(len=50), intent(in) :: seedname
logical, intent(in) :: have_disentangled
logical, intent(in) :: effective_model
! local variables
complex(kind=dp), allocatable :: UU(:, :)
complex(kind=dp), allocatable :: AA(:, :, :), AA_bar(:, :, :)
complex(kind=dp), allocatable :: AA_da(:, :, :, :), AA_da_bar(:, :, :, :)
complex(kind=dp), allocatable :: HH_da(:, :, :), HH_da_bar(:, :, :)
complex(kind=dp), allocatable :: HH_dadb(:, :, :, :), HH_dadb_bar(:, :, :, :)
complex(kind=dp), allocatable :: HH(:, :)
complex(kind=dp), allocatable :: D_h(:, :, :), D_h_no_eta(:, :, :)
real(kind=dp), allocatable :: eig(:)
real(kind=dp), allocatable :: eig_da(:, :)
real(kind=dp), allocatable :: occ(:)
real(kind=dp) :: recip_lattice(3, 3), volume
complex(kind=dp) :: sum_AD(3, 3), sum_HD(3, 3), r_mn(3), gen_r_nm(3)
integer :: a, b, c, bc, n, m, istart, iend
integer :: p ! i, if, r, ifreq
real(kind=dp) :: I_nm(3, 6)
real(kind=dp) :: omega(pw90_berry%kubo_nfreq), delta(pw90_berry%kubo_nfreq), joint_level_spacing
real(kind=dp) :: eta_smr, Delta_k, vdum(3), occ_fac, wstep, wmin, wmax
allocate (UU(num_wann, num_wann))
allocate (AA(num_wann, num_wann, 3))
allocate (AA_bar(num_wann, num_wann, 3))
allocate (AA_da(num_wann, num_wann, 3, 3))
allocate (AA_da_bar(num_wann, num_wann, 3, 3))
allocate (HH_da(num_wann, num_wann, 3))
allocate (HH_da_bar(num_wann, num_wann, 3))
allocate (HH_dadb(num_wann, num_wann, 3, 3))
allocate (HH_dadb_bar(num_wann, num_wann, 3, 3))
allocate (HH(num_wann, num_wann))
allocate (D_h(num_wann, num_wann, 3))
allocate (D_h_no_eta(num_wann, num_wann, 3))
allocate (eig(num_wann))
allocate (occ(num_wann))
allocate (eig_da(num_wann, 3))
! Initialize shift current array at point k
sc_k_list = 0.d0
! Gather W-gauge matrix objects !
! choose the convention for the FT sums
if (pw90_berry%sc_phase_conv .eq. 1) then ! use Wannier centres in the FT exponentials (so called TB convention)
! get Hamiltonian and its first and second derivatives
! Note that below we calculate the UU matrix--> we have to use the same UU from here on for
! maintaining the gauge-covariance of the whole matrix element
call wham_get_eig_UU_HH_AA_sc_TB_conv(pw90_berry, dis_manifold, kmesh_info, kpt_latt, &
ws_region, print_output, wannier_data, ws_distance, &
wigner_seitz, AA_R, HH, HH_da, HH_dadb, HH_R, &
u_matrix, UU, v_matrix, 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
! get position operator and its derivative
! note that AA_da(:,:,a,b) \propto \sum_R exp(iRk)*iR_{b}*<0|r_{a}|R>
call pw90common_fourier_R_to_k_vec_dadb_TB_conv(ws_region, wannier_data, ws_distance, &
wigner_seitz, AA_R, kpt, real_lattice, &
mp_grid, num_wann, error, comm, OO_da=AA, &
OO_dadb=AA_da)
if (allocated(error)) return
! get eigenvalues and their k-derivatives
call wham_get_eig_deleig_TB_conv(pw90_band_deriv_degen, HH_da, UU, eig, eig_da, num_wann, &
error, comm)
if (allocated(error)) return
elseif (pw90_berry%sc_phase_conv .eq. 2) then ! do not use Wannier centres in the FT exponentials (usual W90 convention)
! same as above
call wham_get_eig_UU_HH_AA_sc(dis_manifold, kpt_latt, ws_region, print_output, wannier_data, &
ws_distance, wigner_seitz, HH, HH_da, HH_dadb, HH_R, u_matrix, UU, &
v_matrix, 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 pw90common_fourier_R_to_k_vec_dadb(ws_region, wannier_data, ws_distance, wigner_seitz, &
AA_R, kpt, real_lattice, mp_grid, num_wann, error, &
comm, OO_da=AA, OO_dadb=AA_da)
if (allocated(error)) return
call wham_get_eig_deleig(dis_manifold, kpt_latt, pw90_band_deriv_degen, ws_region, print_output, wannier_data, &
ws_distance, wigner_seitz, HH_da, HH, HH_R, u_matrix, UU, v_matrix, &
eig_da, 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
end if
! get electronic occupations
call pw90common_get_occ(fermi_energy_list(1), eig, occ, num_wann)
! get D_h (Eq. (24) WYSV06)
call wham_get_D_h_P_value(pw90_berry, HH_da, D_h, UU, eig, num_wann)
call wham_get_D_h(HH_da, D_h_no_eta, UU, eig, num_wann)
! calculate k-spacing in case of adaptive smearing
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
! rotate quantities from W to H gauge (we follow wham_get_D_h for delHH_bar_i)
do a = 1, 3
! Berry connection A
AA_bar(:, :, a) = utility_rotate(AA(:, :, a), UU, num_wann)
! first derivative of Hamiltonian dH_da
HH_da_bar(:, :, a) = utility_rotate(HH_da(:, :, a), UU, num_wann)
do b = 1, 3
! derivative of Berry connection dA_da
AA_da_bar(:, :, a, b) = utility_rotate(AA_da(:, :, a, b), UU, num_wann)
! second derivative of Hamiltonian d^{2}H_dadb
HH_dadb_bar(:, :, a, b) = utility_rotate(HH_dadb(:, :, a, b), UU, num_wann)
end do
end do
! setup for frequency-related quantities
omega = real(pw90_berry%kubo_freq_list(:), dp)
wmin = omega(1)
wmax = omega(pw90_berry%kubo_nfreq)
wstep = omega(2) - omega(1)
! loop on initial and final bands
do n = 1, num_wann
do m = 1, num_wann
! cycle diagonal matrix elements and bands above the maximum
if (n == m) cycle
if (eig(m) > pw90_berry%kubo_eigval_max .or. eig(n) > pw90_berry%kubo_eigval_max) cycle
! setup T=0 occupation factors
occ_fac = (occ(n) - occ(m))
if (abs(occ_fac) < 1e-10) cycle
! set delta function smearing
if (pw90_berry%kubo_smearing%use_adaptive) then
vdum(:) = eig_da(m, :) - eig_da(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
! restrict to energy window spanning [-sc_w_thr*eta_smr,+sc_w_thr*eta_smr]
! outside this range, the two delta functions are virtually zero
if (((eig(n) - eig(m) + pw90_berry%sc_w_thr*eta_smr < wmin) .or. &
(eig(n) - eig(m) - pw90_berry%sc_w_thr*eta_smr > wmax)) .and. &
((eig(m) - eig(n) + pw90_berry%sc_w_thr*eta_smr < wmin) .or. &
(eig(m) - eig(n) - pw90_berry%sc_w_thr*eta_smr > wmax))) cycle
! first compute the two sums over intermediate states between AA_bar and HH_da_bar with D_h
! appearing in Eqs. (30) and (32) of IATS18
sum_AD = cmplx_0
sum_HD = cmplx_0
do a = 1, 3
do c = 1, 3
! Note that we substract diagonal elements in AA_bar and
! HH_da_bar to match the convention in IATS18
! (diagonals in D_h are automatically zero, so we do not substract them)
sum_AD(c, a) = (utility_zdotu(AA_bar(n, :, c), D_h(:, m, a)) - AA_bar(n, n, c)*D_h(n, m, a)) &
- (utility_zdotu(D_h(n, :, a), AA_bar(:, m, c)) - D_h(n, m, a)*AA_bar(m, m, c))
sum_HD(c, a) = (utility_zdotu(HH_da_bar(n, :, c), D_h(:, m, a)) - HH_da_bar(n, n, c)*D_h(n, m, a)) &
- (utility_zdotu(D_h(n, :, a), HH_da_bar(:, m, c)) - D_h(n, m, a)*HH_da_bar(m, m, c))
end do
end do
! dipole matrix element
r_mn(:) = AA_bar(m, n, :) + cmplx_i*D_h_no_eta(m, n, :)
! loop over direction of generalized derivative
do a = 1, 3
! store generalized derivative as an array on the additional spatial index,
! its composed of 8 terms in total, see Eq (34) combined with (30) and
! (32) of IATS18
gen_r_nm(:) = (AA_da_bar(n, m, :, a) &
+ ((AA_bar(n, n, :) - AA_bar(m, m, :))*D_h_no_eta(n, m, a) + &
(AA_bar(n, n, a) - AA_bar(m, m, a))*D_h_no_eta(n, m, :)) &
- cmplx_i*AA_bar(n, m, :)*(AA_bar(n, n, a) - AA_bar(m, m, a)) &
+ sum_AD(:, a) &
+ cmplx_i*(HH_dadb_bar(n, m, :, a) &
+ sum_HD(:, a) &
+ (D_h_no_eta(n, m, :)*(eig_da(n, a) - eig_da(m, a)) + &
D_h_no_eta(n, m, a)*(eig_da(n, :) - eig_da(m, :)))) &
/(eig(m) - eig(n)))
! Correction term due to finite sc_eta
! See Eq. (19) of Phys. Rev. B 103, 247101 (2021)
if (pw90_berry%sc_use_eta_corr) then
do p = 1, num_wann
if (p == n .or. p == m) cycle
gen_r_nm(:) = gen_r_nm(:) &
- pw90_berry%sc_eta**2/((eig(p) - eig(m))**2 &
+ pw90_berry%sc_eta**2)/(eig(n) - eig(m)) &
*(AA_bar(n, p, :)*HH_da_bar(p, m, a) &
- (HH_da_bar(n, p, :) + cmplx_i*(eig(n) &
- eig(p))*AA_bar(n, p, :))*AA_bar(p, m, a)) &
+ pw90_berry%sc_eta**2/((eig(n) - eig(p))**2 &
+ pw90_berry%sc_eta**2)/(eig(n) - eig(m)) &
*(HH_da_bar(n, p, a)*AA_bar(p, m, :) &
- AA_bar(n, p, a)*(HH_da_bar(p, m, :) + cmplx_i*(eig(p) - eig(m))*AA_bar(p, m, :)))
end do
end if
! loop over the remaining two indexes of the matrix product.
! Note that shift current is symmetric under b <--> c exchange,
! so we avoid computing all combinations using alpha_S and beta_S
do bc = 1, 6
b = alpha_S(bc)
c = beta_S(bc)
I_nm(a, bc) = aimag(r_mn(b)*gen_r_nm(c) + r_mn(c)*gen_r_nm(b))
end do ! bc
end do ! a
! compute delta(E_nm-w)
! choose energy window spanning [-sc_w_thr*eta_smr,+sc_w_thr*eta_smr]
istart = max(int((eig(n) - eig(m) - pw90_berry%sc_w_thr*eta_smr - wmin)/wstep + 1), 1)
iend = min(int((eig(n) - eig(m) + pw90_berry%sc_w_thr*eta_smr - wmin)/wstep + 1), pw90_berry%kubo_nfreq)
! multiply matrix elements with delta function for the relevant frequencies
if (istart <= iend) then
delta = 0.0
delta(istart:iend) = &
utility_w0gauss_vec((eig(m) - eig(n) + omega(istart:iend))/eta_smr, &
pw90_berry%kubo_smearing%type_index, error, comm)/eta_smr
if (allocated(error)) return
call DGER(18, iend - istart + 1, occ_fac, I_nm, 1, delta(istart:iend), 1, sc_k_list(:, :, istart:iend), 18)
end if
! same for delta(E_mn-w)
istart = max(int((eig(m) - eig(n) - pw90_berry%sc_w_thr*eta_smr - wmin)/wstep + 1), 1)
iend = min(int((eig(m) - eig(n) + pw90_berry%sc_w_thr*eta_smr - wmin)/wstep + 1), pw90_berry%kubo_nfreq)
if (istart <= iend) then
delta = 0.0
delta(istart:iend) = &
utility_w0gauss_vec((eig(n) - eig(m) + omega(istart:iend))/eta_smr, &
pw90_berry%kubo_smearing%type_index, error, comm)/eta_smr
if (allocated(error)) return
call DGER(18, iend - istart + 1, occ_fac, I_nm, 1, delta(istart:iend), 1, sc_k_list(:, :, istart:iend), 18)
end if
end do ! bands
end do ! bands
end subroutine berry_get_sc_klist
