computes <0n|H|Rm>, in eV (pwscf uses Ry, but pw2wannier90 converts to eV)
subroutine get_HH_R
!======================================================
!
!! computes <0n|H|Rm>, in eV
!! (pwscf uses Ry, but pw2wannier90 converts to eV)
!
!======================================================
use w90_constants, only: dp, cmplx_0
use w90_io, only: io_error, stdout, io_stopwatch, &
io_file_unit, seedname
use w90_parameters, only: num_wann, ndimwin, num_kpts, num_bands, &
eigval, u_matrix, have_disentangled, &
timing_level, scissors_shift, &
num_valence_bands, effective_model, &
real_lattice
use w90_postw90_common, only: nrpts, rpt_origin, v_matrix, ndegen, irvec, crvec
use w90_comms, only: on_root, comms_bcast
integer :: i, j, n, m, ii, ik, winmin_q, file_unit, &
ir, io, idum, ivdum(3), ivdum_old(3)
integer, allocatable :: num_states(:)
real(kind=dp) :: rdum_real, rdum_imag
complex(kind=dp), allocatable :: HH_q(:, :, :)
logical :: new_ir
!ivo
complex(kind=dp), allocatable :: sciss_q(:, :, :)
complex(kind=dp), allocatable :: sciss_R(:, :, :)
real(kind=dp) :: sciss_shift
if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_HH_R', 1)
if (.not. allocated(HH_R)) then
allocate (HH_R(num_wann, num_wann, nrpts))
else
if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_HH_R', 2)
return
end if
! Real-space Hamiltonian H(R) is read from file
!
if (effective_model) then
HH_R = cmplx_0
if (on_root) then
write (stdout, '(/a)') ' Reading real-space Hamiltonian from file ' &
//trim(seedname)//'_HH_R.dat'
file_unit = io_file_unit()
open (file_unit, file=trim(seedname)//'_HH_R.dat', form='formatted', &
status='old', err=101)
read (file_unit, *) ! header
read (file_unit, *) idum ! num_wann
read (file_unit, *) idum ! nrpts
ir = 1
new_ir = .true.
ivdum_old(:) = 0
n = 1
do
read (file_unit, '(5I5,2F12.6)', iostat=io) ivdum(1:3), j, i, &
rdum_real, rdum_imag
if (io < 0) exit ! reached end of file
if (i < 1 .or. i > num_wann .or. j < 1 .or. j > num_wann) then
write (stdout, *) 'num_wann=', num_wann, ' i=', i, ' j=', j
call io_error &
('Error in get_HH_R: orbital indices out of bounds')
endif
if (n > 1) then
if (ivdum(1) /= ivdum_old(1) .or. ivdum(2) /= ivdum_old(2) .or. &
ivdum(3) /= ivdum_old(3)) then
ir = ir + 1
new_ir = .true.
else
new_ir = .false.
endif
endif
ivdum_old = ivdum
! Note that the same (j,i,ir) may occur more than once in
! the file seedname_HH_R.dat, hence the addition instead
! of a simple equality. (This has to do with the way the
! Berlijn effective Hamiltonian algorithm is
! implemented.)
HH_R(j, i, ir) = HH_R(j, i, ir) + cmplx(rdum_real, rdum_imag, kind=dp)
if (new_ir) then
irvec(:, ir) = ivdum(:)
if (ivdum(1) == 0 .and. ivdum(2) == 0 .and. ivdum(3) == 0) rpt_origin = ir
endif
n = n + 1
enddo
close (file_unit)
if (ir /= nrpts) then
write (stdout, *) 'ir=', ir, ' nrpts=', nrpts
call io_error('Error in get_HH_R: inconsistent nrpts values')
endif
do ir = 1, nrpts
crvec(:, ir) = matmul(transpose(real_lattice), irvec(:, ir))
end do
ndegen(:) = 1 ! This is assumed when reading HH_R from file
!
! TODO: Implement scissors in this case? Need to choose a
! uniform k-mesh (the scissors correction is applied in
! k-space) and then proceed as below, Fourier transforming
! back to real space and adding to HH_R, Hopefully the
! result converges (rapidly) with the k-mesh density, but
! one should check
!
if (abs(scissors_shift) > 1.0e-7_dp) &
call io_error( &
'Error in get_HH_R: scissors shift not implemented for ' &
//'effective_model=T')
endif
call comms_bcast(HH_R(1, 1, 1), num_wann*num_wann*nrpts)
call comms_bcast(ndegen(1), nrpts)
call comms_bcast(irvec(1, 1), 3*nrpts)
call comms_bcast(crvec(1, 1), 3*nrpts)
if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_HH_R', 2)
return
endif
! Everything below is only executed if effective_model==False (default)
! Real-space Hamiltonian H(R) is calculated by Fourier
! transforming H(q) defined on the ab-initio reciprocal mesh
!
allocate (HH_q(num_wann, num_wann, num_kpts))
allocate (num_states(num_kpts))
HH_q = cmplx_0
do ik = 1, num_kpts
if (have_disentangled) then
num_states(ik) = ndimwin(ik)
else
num_states(ik) = num_wann
endif
call get_win_min(ik, winmin_q)
do m = 1, num_wann
do n = 1, m
do i = 1, num_states(ik)
ii = winmin_q + i - 1
HH_q(n, m, ik) = HH_q(n, m, ik) &
+ conjg(v_matrix(i, n, ik))*eigval(ii, ik) &
*v_matrix(i, m, ik)
enddo
HH_q(m, n, ik) = conjg(HH_q(n, m, ik))
enddo
enddo
enddo
call fourier_q_to_R(HH_q, HH_R)
! Scissors correction for an insulator: shift conduction bands upwards by
! scissors_shift eV
!
if (num_valence_bands > 0 .and. abs(scissors_shift) > 1.0e-7_dp) then
allocate (sciss_R(num_wann, num_wann, nrpts))
allocate (sciss_q(num_wann, num_wann, num_kpts))
sciss_q = cmplx_0
do ik = 1, num_kpts
do j = 1, num_wann
do i = 1, j
do m = 1, num_valence_bands
sciss_q(i, j, ik) = sciss_q(i, j, ik) - &
conjg(u_matrix(m, i, ik))*u_matrix(m, j, ik)
enddo
sciss_q(j, i, ik) = conjg(sciss_q(i, j, ik))
enddo
enddo
enddo
call fourier_q_to_R(sciss_q, sciss_R)
do n = 1, num_wann
sciss_R(n, n, rpt_origin) = sciss_R(n, n, rpt_origin) + 1.0_dp
end do
sciss_R = sciss_R*scissors_shift
HH_R = HH_R + sciss_R
endif
if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_HH_R', 2)
return
101 call io_error('Error in get_HH_R: problem opening file '// &
trim(seedname)//'_HH_R.dat')
end subroutine get_HH_R