!-*- mode: F90 -*-!
!------------------------------------------------------------!
! Copyright (C) 2026 Wannier Developer Group !
! !
! This library is free software; you can redistribute it !
! and/or modify it under the terms of the GNU Lesser General !
! Public License as published by the Free Software !
! Foundation; either version 2.1 of the License, or (at your !
! option) any later version. !
! !
! This library is distributed in the hope that it will be !
! useful,but WITHOUT ANY WARRANTY; without even the implied !
! warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR !
! PURPOSE. See the GNU Lesser General Public License for !
! more details. !
! !
! You should have received a copy of the GNU Lesser General !
! Public License along with this library; if not, see !
! <https://www.gnu.org/licenses/>. !
! !
! The webpage of the Wannier90 code is !
! <https://www.wannier.org>. !
! !
! The Wannier90 code is hosted on GitHub !
! <https://github.com/wannier-developers/wannier90> !
!------------------------------------------------------------!
! !
! w90_wan_ham: Hamiltonian operations in Wannier basis !
! !
!------------------------------------------------------------!
module w90_wan_ham
!! This module contain operations on the Hamiltonian in the WF basis
use w90_constants, only: dp
use w90_error, only: w90_error_type, set_error_alloc, set_error_dealloc, set_error_fatal, &
set_error_input, set_error_fatal, set_error_file
implicit none
private
public :: wham_get_D_h
public :: wham_get_D_h_P_value
public :: wham_get_eig_deleig
public :: wham_get_eig_deleig_TB_conv
public :: wham_get_eig_UU_HH_AA_sc
public :: wham_get_eig_UU_HH_AA_sc_TB_conv
public :: wham_get_eig_UU_HH_JJlist
public :: wham_get_occ_mat_list
contains
!================================================!
subroutine wham_get_D_h_a(delHH_a, UU, eig, ef, D_h_a, num_wann)
!================================================!
!
!! Compute D^H_a=UU^dag.del_a UU (a=alpha,beta),
!! using Eq.(24) of WYSV06
!
!================================================!
use w90_constants, only: dp, cmplx_0
use w90_utility, only: utility_rotate
use w90_postw90_common, only: pw90common_get_occ
! arguments
complex(kind=dp), intent(in) :: delHH_a(:, :)
complex(kind=dp), intent(in) :: UU(:, :)
complex(kind=dp), intent(out) :: D_h_a(:, :)
real(kind=dp), intent(in) :: eig(:)
real(kind=dp), intent(in) :: ef
integer, intent(in) :: num_wann
! local variables
complex(kind=dp), allocatable :: delHH_a_bar(:, :)
real(kind=dp) :: occ(num_wann)
integer :: n, m
call pw90common_get_occ(ef, eig, occ, num_wann)
allocate (delHH_a_bar(num_wann, num_wann))
delHH_a_bar = utility_rotate(delHH_a, UU, num_wann)
do m = 1, num_wann
do n = 1, num_wann
if (occ(n) > 0.999_dp .and. occ(m) < 0.001_dp) then
D_h_a(n, m) = delHH_a_bar(n, m)/(eig(m) - eig(n))
else
D_h_a(n, m) = cmplx_0
end if
end do
end do
D_h_a = D_h_a - conjg(transpose(D_h_a))
end subroutine wham_get_D_h_a
!================================================!
subroutine wham_get_D_h(delHH, D_h, UU, eig, num_wann)
!================================================!
!
!! Compute D^H_a=UU^dag.del_a UU (a=x,y,z)
!! using Eq.(24) of WYSV06
!
!================================================!
! TO DO: Implement version where energy denominators only connect
! occupied and empty states. In this case probably do not need
! to worry about avoiding small energy denominators
use w90_constants, only: dp, cmplx_0
use w90_utility, only: utility_rotate
! arguments
complex(kind=dp), intent(in) :: delHH(:, :, :)
complex(kind=dp), intent(in) :: UU(:, :)
complex(kind=dp), intent(out) :: D_h(:, :, :)
real(kind=dp), intent(in) :: eig(:)
integer, intent(in) :: num_wann
! local variables
complex(kind=dp), allocatable :: delHH_bar_i(:, :)
integer :: n, m, i
allocate (delHH_bar_i(num_wann, num_wann))
D_h = cmplx_0
do i = 1, 3
delHH_bar_i(:, :) = utility_rotate(delHH(:, :, i), UU, num_wann)
do m = 1, num_wann
do n = 1, num_wann
if (n == m .or. abs(eig(m) - eig(n)) < 1.0e-7_dp) cycle
D_h(n, m, i) = delHH_bar_i(n, m)/(eig(m) - eig(n))
end do
end do
end do
end subroutine wham_get_D_h
!================================================!
subroutine wham_get_D_h_P_value(pw90_berry, delHH, D_h, UU, eig, num_wann)
!================================================!
!
!! Compute D^H_a=UU^dag.del_a UU (a=x,y,z)
!! using Eq.(24) of WYSV06
! and prescription for energy denominator
! from BK81
!
!================================================!
! TO DO: Implement version where energy denominators only connect
! occupied and empty states. In this case probably do not need
! to worry about avoiding small energy denominators
use w90_constants, only: dp, cmplx_0
use w90_postw90_types, only: pw90_berry_mod_type !sc_eta
use w90_utility, only: utility_rotate
! arguments
type(pw90_berry_mod_type), intent(in) :: pw90_berry
integer, intent(in) :: num_wann
real(kind=dp), intent(in) :: eig(:)
complex(kind=dp), intent(in) :: delHH(:, :, :)
complex(kind=dp), intent(in) :: UU(:, :)
complex(kind=dp), intent(out) :: D_h(:, :, :)
! local variables
complex(kind=dp), allocatable :: delHH_bar_i(:, :)
integer :: n, m, i
real(kind=dp) :: deltaE
allocate (delHH_bar_i(num_wann, num_wann))
D_h = cmplx_0
deltaE = 0.d0
do i = 1, 3
delHH_bar_i(:, :) = utility_rotate(delHH(:, :, i), UU, num_wann)
do m = 1, num_wann
do n = 1, num_wann
if (n == m) cycle
deltaE = eig(m) - eig(n)
D_h(n, m, i) = delHH_bar_i(n, m)*(deltaE/(deltaE**(2) + pw90_berry%sc_eta**(2)))
end do
end do
end do
end subroutine wham_get_D_h_P_value
!================================================!
subroutine wham_get_JJp_JJm_list(delHH, UU, eig, JJp_list, JJm_list, num_wann, &
fermi_energy_list, occ)
!================================================!
! !
! Compute JJ^+_a and JJ^-_a (a=Cartesian index) !
! for a list of Fermi energies !
! !
! This routine is a replacement for !
! wham_get_JJp_list and wham_getJJm_list. !
! It computes both lists at once in a more !
! efficient manner. !
! !
! Tsirkin: added the optional occ parameter !
! !
!================================================!
use w90_constants, only: dp, cmplx_0, cmplx_i
use w90_utility, only: utility_rotate_new
! arguments
real(kind=dp), allocatable, intent(in) :: fermi_energy_list(:)
integer, intent(in) :: num_wann
real(kind=dp), intent(in) :: eig(:)
real(kind=dp), intent(in), optional, dimension(:) :: occ
complex(kind=dp), intent(inout) :: delHH(:, :)
complex(kind=dp), intent(in) :: UU(:, :)
complex(kind=dp), intent(out) :: JJm_list(:, :, :)
complex(kind=dp), intent(out) :: JJp_list(:, :, :)
! local variables
integer :: n, m, ife, nfermi_loc
real(kind=dp) :: fe
if (present(occ)) then
nfermi_loc = 1
else
nfermi_loc = 0
if (allocated(fermi_energy_list)) nfermi_loc = size(fermi_energy_list)
end if
call utility_rotate_new(delHH, UU, num_wann)
do ife = 1, nfermi_loc
fe = fermi_energy_list(ife)
do m = 1, num_wann
do n = 1, num_wann
if (present(occ)) then
if (occ(m) < 0.5_dp .and. occ(n) > 0.5_dp) then
JJm_list(n, m, ife) = cmplx_i*delHH(n, m)/(eig(m) - eig(n))
JJp_list(m, n, ife) = cmplx_i*delHH(m, n)/(eig(n) - eig(m))
else
JJm_list(n, m, ife) = cmplx_0
JJp_list(m, n, ife) = cmplx_0
end if
else
if (eig(n) > fe .and. eig(m) < fe) then
JJp_list(n, m, ife) = cmplx_i*delHH(n, m)/(eig(m) - eig(n))
JJm_list(m, n, ife) = cmplx_i*delHH(m, n)/(eig(n) - eig(m))
else
JJp_list(n, m, ife) = cmplx_0
JJm_list(m, n, ife) = cmplx_0
end if
end if
end do
end do
call utility_rotate_new(JJp_list(:, :, ife), UU, num_wann, reverse=.true.)
call utility_rotate_new(JJm_list(:, :, ife), UU, num_wann, reverse=.true.)
end do
end subroutine wham_get_JJp_JJm_list
!================================================!
subroutine wham_get_occ_mat_list(fermi_energy_list, f_list, g_list, UU, num_wann, error, comm, &
eig, occ)
!================================================!
!
!! Occupation matrix f, and g=1-f
!! for a list of Fermi energies
! Tsirkin: !now optionally either eig or occ parameters may be supplied
! (Changed consistently the calls from the Berry module)
!================================================!
use w90_constants, only: dp, cmplx_0, cmplx_1
use w90_postw90_common, only: pw90common_get_occ
use w90_comms, only: w90_comm_type
! arguments
type(w90_error_type), allocatable, intent(out) :: error
type(w90_comm_type), intent(in) :: comm
real(kind=dp), allocatable, intent(in) :: fermi_energy_list(:)
integer, intent(in) :: num_wann
real(kind=dp), intent(in), optional :: eig(:)
real(kind=dp), intent(in), optional :: occ(:)
complex(kind=dp), intent(in) :: UU(:, :)
complex(kind=dp), intent(out) :: f_list(:, :, :)
complex(kind=dp), intent(out) :: g_list(:, :, :)
! local variables
integer :: n, m, i, if, nfermi_loc
real(kind=dp), allocatable :: occ_list(:, :)
if (present(occ)) then
nfermi_loc = 1
else
nfermi_loc = 0
if (allocated(fermi_energy_list)) nfermi_loc = size(fermi_energy_list)
end if
allocate (occ_list(num_wann, nfermi_loc))
if (present(occ) .and. present(eig)) then
call set_error_input(error, 'occ_list and eig cannot be both arguments in get_occ_mat_list', &
comm)
return
elseif (.not. present(occ) .and. .not. present(eig)) then
call set_error_input(error, 'either occ_list or eig must be passed as arguments to get_occ_mat_list', comm)
return
end if
if (present(occ)) then
occ_list(:, 1) = occ(:)
else
do if = 1, nfermi_loc
call pw90common_get_occ(fermi_energy_list(if), eig, occ_list(:, if), num_wann)
end do
end if
f_list = cmplx_0
do if = 1, nfermi_loc
do n = 1, num_wann
do m = 1, num_wann
do i = 1, num_wann
f_list(n, m, if) = f_list(n, m, if) &
+ UU(n, i)*occ_list(i, if)*conjg(UU(m, i))
end do
g_list(n, m, if) = -f_list(n, m, if)
if (m == n) g_list(n, n, if) = g_list(n, n, if) + cmplx_1
end do
end do
end do
end subroutine wham_get_occ_mat_list
!================================================!
subroutine wham_get_deleig_a(deleig_a, eig, delHH_a, UU, num_wann, pw90_band_deriv_degen, &
error, comm)
!================================================!
!
!! Band derivatives dE/dk_a
!
!================================================!
use w90_constants, only: dp !, cmplx_0, cmplx_i
use w90_utility, only: utility_diagonalize, utility_rotate, utility_rotate_diag
use w90_postw90_types, only: pw90_band_deriv_degen_type
use w90_comms, only: w90_comm_type
! arguments
type(pw90_band_deriv_degen_type), intent(in) :: pw90_band_deriv_degen
integer, intent(in) :: num_wann
real(kind=dp), intent(in) :: eig(num_wann)
real(kind=dp), intent(out) :: deleig_a(num_wann)
complex(kind=dp), intent(in) :: delHH_a(:, :)
complex(kind=dp), intent(in) :: UU(:, :)
type(w90_error_type), allocatable, intent(out) :: error
type(w90_comm_type), intent(in) :: comm
! local variables
integer :: i, degen_min, degen_max, dim
real(kind=dp) :: diff
complex(kind=dp), allocatable :: delHH_bar_a(:, :), U_deg(:, :)
allocate (delHH_bar_a(num_wann, num_wann))
allocate (U_deg(num_wann, num_wann))
if (pw90_band_deriv_degen%use_degen_pert) then
delHH_bar_a = utility_rotate(delHH_a, UU, num_wann)
! Assuming that the energy eigenvalues are stored in eig(:) in
! increasing order (diff >= 0)
i = 0
do
i = i + 1
if (i > num_wann) exit
if (i + 1 <= num_wann) then
diff = eig(i + 1) - eig(i)
else
! i-th is the highest band, and it is non-degenerate
diff = pw90_band_deriv_degen%degen_thr + 1.0_dp
end if
if (diff < pw90_band_deriv_degen%degen_thr) then
! Bands i and i+1 are degenerate
degen_min = i
degen_max = degen_min + 1
! See if any higher bands are in the same degenerate group
do
if (degen_max + 1 > num_wann) exit
diff = eig(degen_max + 1) - eig(degen_max)
if (diff < pw90_band_deriv_degen%degen_thr) then
degen_max = degen_max + 1
else
exit
end if
end do
! Bands from degen_min to degen_max are degenerate. Diagonalize
! the submatrix in Eq.(31) YWVS07 over this degenerate subspace.
! The eigenvalues are the band gradients
dim = degen_max - degen_min + 1
call utility_diagonalize(delHH_bar_a(degen_min:degen_max, degen_min:degen_max), dim, &
deleig_a(degen_min:degen_max), U_deg(1:dim, 1:dim), error, comm)
if (allocated(error)) return
! Scanned bands up to degen_max
i = degen_max
else
! Use non-degenerate form [Eq.(27) YWVS07] for current (i-th) band
deleig_a(i) = real(delHH_bar_a(i, i), dp)
end if
end do
else
! Use non-degenerate form for all bands
deleig_a(:) = real(utility_rotate_diag(delHH_a(:, :), UU, num_wann), dp)
end if
end subroutine wham_get_deleig_a
!================================================!
subroutine 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)
!================================================!
!
!! Given a k point, this function returns eigenvalues E and
!! derivatives of the eigenvalues dE/dk_a, using wham_get_deleig_a
!
!================================================!
use w90_constants, only: dp
use w90_postw90_types, only: pw90_band_deriv_degen_type, wigner_seitz_type
use w90_comms, only: w90_comm_type, mpirank
use w90_constants, only: dp, cmplx_0
use w90_get_oper, only: get_HH_R
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
use w90_utility, only: utility_diagonalize
implicit none
! arguments
type(dis_manifold_type), intent(in) :: dis_manifold
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) :: num_wann, num_kpts, num_bands, num_valence_bands
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: stdout
real(kind=dp), intent(in) :: kpt(3)!! the three coordinates of the k point vector (in relative coordinates)
real(kind=dp), intent(out) :: eig(num_wann)!! the calculated eigenvalues at kpt
real(kind=dp), intent(out) :: del_eig(num_wann, 3)
!! the calculated derivatives of the eigenvalues at kpt [first component: band; second component: 1,2,3
!! for the derivatives along the three k directions]
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(out) :: HH(:, :)
!! the Hamiltonian matrix at kpt
complex(kind=dp), intent(out) :: delHH(:, :, :)
!! the delHH matrix (derivative of H) at kpt
complex(kind=dp), intent(out) :: UU(:, :)
!! the rotation matrix that gives the eigenvectors of HH
complex(kind=dp), intent(in) :: u_matrix(:, :, :), v_matrix(:, :, :)
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
! I call it to be sure that it has been called already once,
! and that HH_R contains the actual matrix.
! Further calls should return very fast.
call get_HH_R(dis_manifold, kpt_latt, print_output, wigner_seitz, HH_R, u_matrix, v_matrix, &
eigval, real_lattice, scissors_shift, num_bands, num_kpts, num_wann, &
num_valence_bands, effective_model, have_disentangled, seedname, ws_distance, ws_region, &
stdout, timer, error, comm)
if (allocated(error)) return
call pw90common_fourier_R_to_k(ws_region, wannier_data, ws_distance, wigner_seitz, HH, HH_R, &
kpt, real_lattice, mp_grid, 0, num_wann, error, comm)
if (allocated(error)) return
call utility_diagonalize(HH, num_wann, eig, UU, error, comm)
if (allocated(error)) return
call pw90common_fourier_R_to_k(ws_region, wannier_data, ws_distance, wigner_seitz, &
delHH(:, :, 1), HH_R, kpt, real_lattice, mp_grid, 1, num_wann, &
error, comm)
if (allocated(error)) return
call pw90common_fourier_R_to_k(ws_region, wannier_data, ws_distance, wigner_seitz, &
delHH(:, :, 2), HH_R, kpt, real_lattice, mp_grid, 2, num_wann, &
error, comm)
if (allocated(error)) return
call pw90common_fourier_R_to_k(ws_region, wannier_data, ws_distance, wigner_seitz, &
delHH(:, :, 3), HH_R, kpt, real_lattice, mp_grid, 3, num_wann, &
error, comm)
if (allocated(error)) return
call wham_get_deleig_a(del_eig(:, 1), eig, delHH(:, :, 1), UU, num_wann, &
pw90_band_deriv_degen, error, comm)
if (allocated(error)) return
call wham_get_deleig_a(del_eig(:, 2), eig, delHH(:, :, 2), UU, num_wann, &
pw90_band_deriv_degen, error, comm)
if (allocated(error)) return
call wham_get_deleig_a(del_eig(:, 3), eig, delHH(:, :, 3), UU, num_wann, &
pw90_band_deriv_degen, error, comm)
if (allocated(error)) return
end subroutine wham_get_eig_deleig
!================================================!
subroutine wham_get_eig_deleig_TB_conv(pw90_band_deriv_degen, delHH, UU, eig, del_eig, num_wann, &
error, comm)
!================================================!
! modified version of wham_get_eig_deleig for the TB convention
! avoids recalculating delHH and UU, works with input values
!
!! Given a k point, this function returns eigenvalues E and
!! derivatives of the eigenvalues dE/dk_a, using wham_get_deleig_a
!
!================================================!
use w90_postw90_types, only: pw90_band_deriv_degen_type
use w90_comms, only: w90_comm_type
! arguments
type(pw90_band_deriv_degen_type), intent(in) :: pw90_band_deriv_degen
type(w90_error_type), allocatable, intent(out) :: error
type(w90_comm_type), intent(in) :: comm
integer, intent(in) :: num_wann
real(kind=dp), intent(out) :: del_eig(num_wann, 3)
real(kind=dp), intent(in) :: eig(num_wann)
complex(kind=dp), intent(in) :: delHH(:, :, :)
!! the delHH matrix (derivative of H) at kpt
complex(kind=dp), intent(in) :: UU(:, :)
!! the rotation matrix that gives the eigenvectors of HH
call wham_get_deleig_a(del_eig(:, 1), eig, delHH(:, :, 1), UU, num_wann, &
pw90_band_deriv_degen, error, comm)
if (allocated(error)) return
call wham_get_deleig_a(del_eig(:, 2), eig, delHH(:, :, 2), UU, num_wann, &
pw90_band_deriv_degen, error, comm)
if (allocated(error)) return
call wham_get_deleig_a(del_eig(:, 3), eig, delHH(:, :, 3), UU, num_wann, &
pw90_band_deriv_degen, error, comm)
if (allocated(error)) return
end subroutine wham_get_eig_deleig_TB_conv
!================================================!
subroutine wham_get_eig_UU_HH_JJlist(dis_manifold, fermi_energy_list, kpt_latt, ws_region, &
print_output, wannier_data, ws_distance, wigner_seitz, HH, &
HH_R, JJm_list, JJp_list, 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, occ)
!================================================!
!
!! Wrapper routine used to reduce number of Fourier calls
! Added the optional occ parameter
!
!================================================!
use w90_constants, only: dp
use w90_postw90_common, only: pw90common_fourier_R_to_k_new_second_d, &
pw90common_fourier_R_to_k_new
use w90_get_oper, only: get_HH_R
use w90_utility, only: utility_diagonalize
use w90_types, only: print_output_type, wannier_data_type, dis_manifold_type, &
ws_region_type, ws_distance_type, timer_list_type
use w90_comms, only: w90_comm_type, mpirank
use w90_postw90_types, only: wigner_seitz_type
implicit none
! arguments
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(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_kpts, num_bands, num_valence_bands
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: stdout
real(kind=dp), intent(in) :: kpt(3), real_lattice(3, 3)
real(kind=dp), intent(out) :: eig(:)
real(kind=dp), intent(in) :: eigval(:, :)
real(kind=dp), intent(in) :: scissors_shift
real(kind=dp), intent(in), optional :: occ(:)
complex(kind=dp), intent(out) :: UU(:, :)
complex(kind=dp), intent(out) :: HH(:, :)
complex(kind=dp), intent(out) :: JJp_list(:, :, :, :)
complex(kind=dp), intent(out) :: JJm_list(:, :, :, :)
complex(kind=dp), intent(in) :: u_matrix(:, :, :), v_matrix(:, :, :)
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
integer :: i
complex(kind=dp), allocatable :: delHH(:, :, :)
call get_HH_R(dis_manifold, kpt_latt, print_output, wigner_seitz, HH_R, u_matrix, v_matrix, &
eigval, real_lattice, scissors_shift, num_bands, num_kpts, num_wann, &
num_valence_bands, effective_model, have_disentangled, seedname, ws_distance, ws_region, &
stdout, timer, error, comm)
if (allocated(error)) return
allocate (delHH(num_wann, num_wann, 3))
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
call utility_diagonalize(HH, num_wann, eig, UU, error, comm)
if (allocated(error)) return
do i = 1, 3
if (present(occ)) then
call wham_get_JJp_JJm_list(delHH(:, :, i), UU, eig, JJp_list(:, :, :, i), &
JJm_list(:, :, :, i), num_wann, fermi_energy_list, occ=occ)
else
call wham_get_JJp_JJm_list(delHH(:, :, i), UU, eig, JJp_list(:, :, :, i), &
JJm_list(:, :, :, i), num_wann, fermi_energy_list)
end if
end do
end subroutine wham_get_eig_UU_HH_JJlist
!================================================!
subroutine 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)
!================================================!
!
! modified version of wham_get_eig_UU_HH_AA_sc, calls routines
! satisfying the TB phase convention
!
!================================================!
use w90_constants, only: dp
use w90_get_oper, only: get_HH_R, get_AA_R_effective, get_AA_R
use w90_postw90_common, only: pw90common_fourier_R_to_k_new_second_d_TB_conv
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_utility, only: utility_diagonalize
use w90_postw90_types, only: pw90_berry_mod_type, wigner_seitz_type
use w90_comms, only: w90_comm_type, mpirank
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
real(kind=dp), intent(in) :: kpt_latt(:, :)
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_kpts, num_bands, num_valence_bands
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: stdout
real(kind=dp), intent(out) :: eig(num_wann)
real(kind=dp), intent(in) :: eigval(:, :)
real(kind=dp), intent(in) :: kpt(3), real_lattice(3, 3)
real(kind=dp), intent(in) :: scissors_shift
complex(kind=dp), intent(out) :: UU(:, :)
complex(kind=dp), intent(out) :: HH(:, :)
complex(kind=dp), intent(out) :: HH_da(:, :, :)
complex(kind=dp), intent(out) :: HH_dadb(:, :, :, :)
complex(kind=dp), intent(in) :: u_matrix(:, :, :), v_matrix(:, :, :)
complex(kind=dp), allocatable, intent(inout) :: HH_R(:, :, :) ! <0n|r|Rm>
complex(kind=dp), allocatable, intent(inout) :: AA_R(:, :, :, :) ! <0n|r|Rm>
character(len=50), intent(in) :: seedname
logical, intent(in) :: have_disentangled
logical, intent(in) :: effective_model
call get_HH_R(dis_manifold, kpt_latt, print_output, wigner_seitz, HH_R, u_matrix, v_matrix, &
eigval, real_lattice, scissors_shift, num_bands, num_kpts, num_wann, &
num_valence_bands, effective_model, have_disentangled, seedname, ws_distance, ws_region, &
stdout, timer, error, comm)
if (allocated(error)) return
if (effective_model) then
call get_AA_R_effective(print_output, AA_R, HH_R, wigner_seitz%nrpts, num_wann, seedname, &
stdout, timer, error, comm)
else
call get_AA_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, wannier_data, AA_R, &
v_matrix, eigval, wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, &
num_wann, have_disentangled, seedname, stdout, timer, error, comm)
end if
if (allocated(error)) return
call pw90common_fourier_R_to_k_new_second_d_TB_conv(kpt, HH_R, AA_R, num_wann, ws_region, &
wannier_data, real_lattice, mp_grid, &
ws_distance, wigner_seitz, error, comm, &
OO=HH, OO_da=HH_da(:, :, :), &
OO_dadb=HH_dadb(:, :, :, :))
if (allocated(error)) return
call utility_diagonalize(HH, num_wann, eig, UU, error, comm)
if (allocated(error)) return
end subroutine wham_get_eig_UU_HH_AA_sc_TB_conv
subroutine 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)
!================================================!
!
!! Wrapper routine used to reduce number of Fourier calls
!
!================================================!
use w90_constants, only: dp
use w90_get_oper, only: get_HH_R
use w90_postw90_common, only: pw90common_fourier_R_to_k_new_second_d
use w90_utility, only: utility_diagonalize
use w90_comms, only: w90_comm_type, mpirank
use w90_types, only: print_output_type, wannier_data_type, dis_manifold_type, &
ws_region_type, ws_distance_type, timer_list_type
use w90_postw90_types, only: wigner_seitz_type
implicit none
! arguments
type(dis_manifold_type), intent(in) :: dis_manifold
real(kind=dp), intent(in) :: kpt_latt(:, :)
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_kpts, num_bands, num_valence_bands
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: stdout
real(kind=dp), intent(out) :: eig(num_wann)
real(kind=dp), intent(in) :: eigval(:, :)
real(kind=dp), intent(in) :: kpt(3), real_lattice(3, 3)
real(kind=dp), intent(in) :: scissors_shift
complex(kind=dp), intent(out) :: UU(:, :)
complex(kind=dp), intent(out) :: HH(:, :)
complex(kind=dp), intent(out) :: HH_da(:, :, :)
complex(kind=dp), intent(out) :: HH_dadb(:, :, :, :)
complex(kind=dp), intent(in) :: u_matrix(:, :, :), v_matrix(:, :, :)
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
call get_HH_R(dis_manifold, kpt_latt, print_output, wigner_seitz, HH_R, u_matrix, v_matrix, &
eigval, real_lattice, scissors_shift, num_bands, num_kpts, num_wann, &
num_valence_bands, effective_model, have_disentangled, seedname, ws_distance, ws_region, &
stdout, timer, error, comm)
if (allocated(error)) return
call pw90common_fourier_R_to_k_new_second_d(kpt, HH_R, num_wann, ws_region, wannier_data, &
real_lattice, mp_grid, ws_distance, wigner_seitz, &
error, comm, OO=HH, OO_da=HH_da(:, :, :), &
OO_dadb=HH_dadb(:, :, :, :))
if (allocated(error)) return
call utility_diagonalize(HH, num_wann, eig, UU, error, comm)
if (allocated(error)) return
end subroutine wham_get_eig_UU_HH_AA_sc
end module w90_wan_ham
