subroutine berry_get_kdotp(kdotp, dis_manifold, kpt_latt, print_output, pw90_berry, &
pw90_band_deriv_degen, wannier_data, ws_distance, wigner_seitz, &
ws_region, HH_R, u_matrix, v_matrix, eigval, real_lattice, &
scissors_shift, mp_grid, num_bands, num_kpts, num_wann, &
num_valence_bands, effective_model, have_disentangled, seedname, &
stdout, timer, error, comm)
!================================================!
! Extracts k.p expansion coefficients using quasi-degenerate
! (Lowdin) perturbation theory, adapted to the Wannier formalism,
! see Appendix in IAdJS19 for details
!================================================!
use w90_constants, only: dp, cmplx_0, cmplx_i
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
use w90_utility, only: utility_rotate
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_mod_type, &
pw90_spin_hall_type, pw90_band_deriv_degen_type, pw90_oper_read_type, wigner_seitz_type, &
kpoint_dist_type
use w90_comms, only: w90_comm_type
implicit none
! Arguments
type(pw90_berry_mod_type), intent(in) :: pw90_berry
type(dis_manifold_type), intent(in) :: dis_manifold
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(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_comm_type), intent(in) :: comm
type(w90_error_type), allocatable, intent(out) :: error
complex(kind=dp), intent(out), dimension(:, :, :, :, :) :: kdotp
!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), intent(in) :: u_matrix(:, :, :), v_matrix(:, :, :)
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), intent(in) :: kpt_latt(:, :)
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: num_wann, num_kpts, num_bands, num_valence_bands
integer, intent(in) :: stdout
character(len=50), intent(in) :: seedname
logical, intent(in) :: have_disentangled
logical, intent(in) :: effective_model
complex(kind=dp), allocatable :: UU(:, :)
complex(kind=dp), allocatable :: HH_da(:, :, :), HH_da_bar(:, :, :)
complex(kind=dp), allocatable :: HH_dadb(:, :, :, :), HH_dadb_bar(:, :, :, :)
complex(kind=dp), allocatable :: HH(:, :), HH_bar(:, :)
real(kind=dp), allocatable :: eig(:)
real(kind=dp), allocatable :: eig_da(:, :)
complex(kind=dp), allocatable :: D_h(:, :, :)
! local variables
!real(kind=dp) :: DeltaE_n, DeltaE_m
integer :: kdotp_num_bands
integer :: i, a, b, n, m, r !, c, bc,if, ifreq, istart, iend
logical :: break_loop
allocate (UU(num_wann, num_wann))
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 (HH_bar(num_wann, num_wann))
allocate (eig(num_wann))
allocate (eig_da(num_wann, 3))
allocate (D_h(num_wann, num_wann, 3))
! Gather W-gauge matrix objects !
! get Hamiltonian and its first and second derivatives
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, &
pw90_berry%kdotp_kpoint, 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 eigenvalues and their k-derivatives
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, &
pw90_berry%kdotp_kpoint, 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 D_h (Eq. (24) WYSV06)
call wham_get_D_h_P_value(pw90_berry, HH_da, D_h, UU, eig, num_wann)
! rotate quantities from W to H gauge
HH_bar(:, :) = utility_rotate(HH(:, :), UU, num_wann)
do a = 1, 3
! first derivative of Hamiltonian dH_da
HH_da_bar(:, :, a) = utility_rotate(HH_da(:, :, a), UU, num_wann)
do b = 1, 3
! 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
kdotp_num_bands = size(pw90_berry%kdotp_bands)
! loop on initial and final bands in k.p set (subset A in IAdJS19)
do n = 1, kdotp_num_bands
do m = 1, kdotp_num_bands
! zeroth order term
if (n == m) kdotp(n, m, 1, 1, 1) = eig(pw90_berry%kdotp_bands(n))
! first order term
do a = 1, 3
kdotp(n, m, 2, a, 1) = HH_da_bar(pw90_berry%kdotp_bands(n), pw90_berry%kdotp_bands(m), a)
end do
! second order term
do a = 1, 3
do b = 1, 3
! add contribution independent of other states
kdotp(n, m, 3, a, b) = 0.5*(HH_dadb_bar(pw90_berry%kdotp_bands(n), &
pw90_berry%kdotp_bands(m), a, b))
! add contribution dependent on other states (subset B in IAdJS19)
do r = 1, num_wann
! cycle for bands in the k.p set (subset A)
break_loop = .false.
do i = 1, kdotp_num_bands
if (r == pw90_berry%kdotp_bands(i)) break_loop = .true.
end do
if (break_loop) cycle
kdotp(n, m, 3, a, b) = kdotp(n, m, 3, a, b) + &
0.5*HH_da_bar(pw90_berry%kdotp_bands(n), r, a) &
*HH_da_bar(r, pw90_berry%kdotp_bands(m), b) &
*((eig(pw90_berry%kdotp_bands(n)) - eig(r))**(-1) &
+ (eig(pw90_berry%kdotp_bands(m)) - eig(r))**(-1))
end do
end do
end do
end do ! bands
end do ! bands
end subroutine berry_get_kdotp
