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