CC_ab(R) = <0|r_a.H.(r-R)_b|R> is the Fourier transform of
CC_ab(k) =
subroutine get_CC_R
!=============================================================
!
!! CC_ab(R) = <0|r_a.H.(r-R)_b|R> is the Fourier transform of
!! CC_ab(k) = <del_a u|H|del_b u> (a,b=x,y,z)
!
!=============================================================
use w90_constants, only: dp, cmplx_0
use w90_parameters, only: num_kpts, nntot, nnlist, num_wann, &
num_bands, ndimwin, wb, bk, &
have_disentangled, timing_level, &
scissors_shift, uHu_formatted
use w90_postw90_common, only: nrpts, v_matrix
use w90_io, only: stdout, io_error, io_stopwatch, io_file_unit, &
seedname
use w90_comms, only: on_root, comms_bcast
integer :: i, j, ii, jj, m, n, a, b, nn1, nn2, ik, nb_tmp, nkp_tmp, &
nntot_tmp, uHu_in, qb1, qb2, winmin_qb1, winmin_qb2
integer, allocatable :: num_states(:)
complex(kind=dp), allocatable :: CC_q(:, :, :, :, :)
complex(kind=dp), allocatable :: Ho_qb1_q_qb2(:, :)
complex(kind=dp), allocatable :: H_qb1_q_qb2(:, :)
real(kind=dp) :: c_real, c_img
character(len=60) :: header
if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_CC_R', 1)
if (.not. allocated(CC_R)) then
allocate (CC_R(num_wann, num_wann, nrpts, 3, 3))
else
if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_CC_R', 2)
return
end if
if (on_root) then
if (abs(scissors_shift) > 1.0e-7_dp) &
call io_error('Error: scissors correction not yet implemented for CC_R')
allocate (Ho_qb1_q_qb2(num_bands, num_bands))
allocate (H_qb1_q_qb2(num_wann, num_wann))
allocate (CC_q(num_wann, num_wann, num_kpts, 3, 3))
allocate (num_states(num_kpts))
do ik = 1, num_kpts
if (have_disentangled) then
num_states(ik) = ndimwin(ik)
else
num_states(ik) = num_wann
endif
enddo
uHu_in = io_file_unit()
if (uHu_formatted) then
open (unit=uHu_in, file=trim(seedname)//".uHu", form='formatted', &
status='old', action='read', err=105)
write (stdout, '(/a)', advance='no') &
' Reading uHu overlaps from '//trim(seedname)//'.uHu in get_CC_R: '
read (uHu_in, *, err=106, end=106) header
write (stdout, '(a)') trim(header)
read (uHu_in, *, err=106, end=106) nb_tmp, nkp_tmp, nntot_tmp
else
open (unit=uHu_in, file=trim(seedname)//".uHu", form='unformatted', &
status='old', action='read', err=105)
write (stdout, '(/a)', advance='no') &
' Reading uHu overlaps from '//trim(seedname)//'.uHu in get_CC_R: '
read (uHu_in, err=106, end=106) header
write (stdout, '(a)') trim(header)
read (uHu_in, err=106, end=106) nb_tmp, nkp_tmp, nntot_tmp
endif
if (nb_tmp .ne. num_bands) &
call io_error &
(trim(seedname)//'.uHu has not the right number of bands')
if (nkp_tmp .ne. num_kpts) &
call io_error &
(trim(seedname)//'.uHu has not the right number of k-points')
if (nntot_tmp .ne. nntot) &
call io_error &
(trim(seedname)//'.uHu has not the right number of nearest neighbours')
CC_q = cmplx_0
do ik = 1, num_kpts
do nn2 = 1, nntot
qb2 = nnlist(ik, nn2)
call get_win_min(qb2, winmin_qb2)
do nn1 = 1, nntot
qb1 = nnlist(ik, nn1)
call get_win_min(qb1, winmin_qb1)
!
! Read from .uHu file the matrices <u_{q+b1}|H_q|u_{q+b2}>
! between the original ab initio eigenstates
!
if (uHu_formatted) then
do m = 1, num_bands
do n = 1, num_bands
read (uHu_in, *, err=106, end=106) c_real, c_img
Ho_qb1_q_qb2(n, m) = cmplx(c_real, c_img, dp)
end do
end do
else
read (uHu_in, err=106, end=106) &
((Ho_qb1_q_qb2(n, m), n=1, num_bands), m=1, num_bands)
endif
! pw2wannier90 is coded a bit strangely, so here we take the transpose
Ho_qb1_q_qb2 = transpose(Ho_qb1_q_qb2)
! old code here
!do m=1,num_bands
! do n=1,num_bands
! read(uHu_in,err=106,end=106) Ho_qb1_q_qb2(m,n)
! end do
!end do
!
! Transform to projected subspace, Wannier gauge
!
call get_gauge_overlap_matrix( &
qb1, num_states(qb1), &
qb2, num_states(qb2), &
Ho_qb1_q_qb2, H_qb1_q_qb2)
do b = 1, 3
do a = 1, b
CC_q(:, :, ik, a, b) = CC_q(:, :, ik, a, b) + wb(nn1)*bk(a, nn1, ik) &
*wb(nn2)*bk(b, nn2, ik)*H_qb1_q_qb2(:, :)
enddo
enddo
enddo !nn1
enddo !nn2
do b = 1, 3
do a = 1, b
CC_q(:, :, ik, b, a) = conjg(transpose(CC_q(:, :, ik, a, b)))
enddo
enddo
enddo !ik
close (uHu_in)
do b = 1, 3
do a = 1, 3
call fourier_q_to_R(CC_q(:, :, :, a, b), CC_R(:, :, :, a, b))
enddo
enddo
endif !on_root
call comms_bcast(CC_R(1, 1, 1, 1, 1), num_wann*num_wann*nrpts*3*3)
if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_CC_R', 2)
return
105 call io_error &
('Error: Problem opening input file '//trim(seedname)//'.uHu')
106 call io_error &
('Error: Problem reading input file '//trim(seedname)//'.uHu')
end subroutine get_CC_R