subroutine get_SAA_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, SS_R, SAA_R, v_matrix, &
scissors_shift, wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, num_wann, &
have_disentangled, seedname, stdout, timer, error, comm)
!================================================!
!
! SAA_ab(R) = <0|s_a.(r-R)_b|R> is the Fourier transform of
! SAA_ab(k) = <u|s_a|del_b u> (a,b=x,y,z)
!
!================================================!
use w90_postw90_types, only: pw90_berry_mod_type, wigner_seitz_type
use w90_types, only: dis_manifold_type, kmesh_info_type, ws_distance_type, ws_region_type, &
print_output_type, timer_list_type
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
type(wigner_seitz_type), intent(in) :: wigner_seitz
type(ws_distance_type), intent(in) :: ws_distance
type(ws_region_type), intent(in) :: ws_region
type(print_output_type), intent(in) :: print_output
type(timer_list_type), intent(inout) :: timer
type(w90_comm_type), intent(in) :: comm
type(w90_error_type), allocatable, intent(out) :: error
integer, intent(in) :: num_bands, num_kpts, num_wann, stdout
!real(kind=dp), intent(in) :: eigval(:, :)
real(kind=dp), intent(in) :: scissors_shift
real(kind=dp), intent(in) :: kpt_latt(:, :)
complex(kind=dp), intent(in) :: v_matrix(:, :, :)
complex(kind=dp), allocatable, intent(in) :: SS_R(:, :, :, :)
complex(kind=dp), allocatable, intent(inout) :: SAA_R(:, :, :, :, :) !<0n|sigma_x,y,z.(r-R)_alpha|Rm>
complex(kind=dp), allocatable :: SAA_R_temp(:, :, :)
complex(kind=dp), allocatable :: SAA_R_b(:, :, :, :, :)
logical, intent(in) :: have_disentangled
character(len=50), intent(in) :: seedname
! local variables
integer :: i, j, ii, jj, ir, m, n, a, b, nn2, ik, nb_tmp, nkp_tmp, &
nntot_tmp, sIu_in, qb2, winmin_q, winmin_qb2, w
integer :: ipol
integer, allocatable :: num_states(:)
integer, allocatable :: counts(:), displs(:)
complex(kind=dp), allocatable :: SAA_q_b(:, :, :, :, :, :)
complex(kind=dp), allocatable :: SAA_q(:, :, :, :, :)
complex(kind=dp), allocatable :: SAA_q_loc(:, :, :)
complex(kind=dp), allocatable :: Ho_q_qb2(:, :, :)
complex(kind=dp), allocatable :: H_q_qb2(:, :)
real(kind=dp), allocatable :: r0(:, :, :)
complex(kind=dp), allocatable :: phase1(:, :), phase2(:)
character(len=60) :: header
logical :: on_root = .false.
if (mpirank(comm) == 0) on_root = .true.
if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
call io_stopwatch_start('get_oper: get_SAA_R', timer)
if (.not. allocated(SAA_R)) then
allocate (SAA_R(num_wann, num_wann, wigner_seitz%nrpts_pw90, 3, 3))
else
if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
call io_stopwatch_stop('get_oper: get_SAA_R', timer)
return
end if
if (on_root) then
allocate (SAA_q_b(num_wann, num_wann, num_kpts, kmesh_info%nntot, 3, 3))
SAA_R = cmplx_0
else
allocate (SAA_q_b(1, 1, 1, kmesh_info%nntot, 3, 3))
end if
if (on_root) then
if (abs(scissors_shift) > 1.0e-7_dp) then
call set_error_fatal(error, 'Error: scissors correction not yet implemented for SAA_R', comm)
return
end if
allocate (Ho_q_qb2(num_bands, num_bands, 3))
allocate (H_q_qb2(num_wann, num_wann))
allocate (num_states(num_kpts))
allocate (phase1(num_wann, num_wann))
if (pw90_berry%transl_inv_full) then
allocate (r0(num_wann, num_wann, 3))
do j = 1, num_wann
do i = 1, num_wann
r0(i, j, :) = (wigner_seitz%wannier_centres_from_AA_R(:, i) + &
wigner_seitz%wannier_centres_from_AA_R(:, j))/2.0_dp
end do
end do
end if
do ik = 1, num_kpts
if (have_disentangled) then
num_states(ik) = dis_manifold%ndimwin(ik)
else
num_states(ik) = num_wann
end if
end do
open (newunit=sIu_in, file=trim(seedname)//".sIu", form='unformatted', &
status='old', action='read', err=113)
write (stdout, '(/a)', advance='no') &
' Reading sIu overlaps from '//trim(seedname)//'.sIu in get_SAA_R: '
read (sIu_in, err=114, end=114) header
write (stdout, '(a)') trim(header)
read (sIu_in, err=114, end=114) nb_tmp, nkp_tmp, nntot_tmp
if (nb_tmp .ne. num_bands) then
call set_error_fatal(error, trim(seedname)//'.sIu has not the right number of bands', comm)
return
end if
if (nkp_tmp .ne. num_kpts) then
call set_error_fatal(error, trim(seedname)//'.sIu has not the right number of k-points', comm)
return
end if
if (nntot_tmp .ne. kmesh_info%nntot) then
call set_error_fatal(error, trim(seedname)//'.sIu has not the right number of nearest neighbours', comm)
return
end if
SAA_q_b = cmplx_0
do ik = 1, num_kpts
call get_win_min(num_bands, dis_manifold, ik, winmin_q, have_disentangled)
do nn2 = 1, kmesh_info%nntot
if (pw90_berry%transl_inv_full) then
phase1 = (r0(:, :, 1)*kmesh_info%bk(1, nn2, ik) + &
r0(:, :, 2)*kmesh_info%bk(2, nn2, ik) + &
r0(:, :, 3)*kmesh_info%bk(3, nn2, ik))
phase1 = exp(cmplx_i*phase1)
else
phase1 = cmplx_1
end if
qb2 = kmesh_info%nnlist(ik, nn2)
call get_win_min(num_bands, dis_manifold, qb2, winmin_qb2, have_disentangled)
do ipol = 1, 3
!
! Read from .sIu file the matrices <u_q|s_a|u_{q+b2}>
! between the original ab initio eigenstates
!
read (sIu_in, err=114, end=114) &
((Ho_q_qb2(n, m, ipol), n=1, num_bands), m=1, num_bands)
! pw2wannier90 is coded a bit strangely, so here we take the transpose
Ho_q_qb2(:, :, ipol) = transpose(Ho_q_qb2(:, :, ipol))
end do
H_q_qb2(:, :) = cmplx_0
do ipol = 1, 3
do m = 1, num_wann
do n = 1, num_wann
do i = 1, num_states(ik)
ii = winmin_q + i - 1
do j = 1, num_states(qb2)
jj = winmin_qb2 + j - 1
H_q_qb2(n, m) = H_q_qb2(n, m) &
+ conjg(v_matrix(i, n, ik)) &
*Ho_q_qb2(ii, jj, ipol) &
*v_matrix(j, m, qb2)
end do
end do
end do
end do
do b = 1, 3
nn2o = kmesh_info%nninv(nn2, ik)
SAA_q_b(:, :, ik, nn2o, ipol, b) = SAA_q_b(:, :, ik, nn2o, ipol, b) + &
cmplx_i*phase1(:, :)*kmesh_info%wb(nn2)*kmesh_info%bk(b, nn2, ik)*H_q_qb2(:, :)
end do
! enddo !nn1
end do !ipol
end do !nn2
end do !ik
close (sIu_in)
end if !on_root
if (pw90_berry%transl_inv_full) then
if (.not. allocated(SS_R)) then
call set_error_fatal(error, 'transl_inv_full=T for SAA_R needs SS_R', comm)
end if
allocate (counts(0:mpisize(comm) - 1))
allocate (displs(0:mpisize(comm) - 1))
w = num_wann*num_wann
call comms_array_split(num_kpts, counts, displs, comm)
allocate (SAA_q_loc(num_wann, num_wann, counts(mpirank(comm))))
allocate (SAA_R_temp(num_wann, num_wann, wigner_seitz%nrpts))
if (on_root) then
allocate (SAA_R_b(num_wann, num_wann, wigner_seitz%nrpts_pw90, 3, 3))
allocate (phase2(wigner_seitz%nrpts_pw90))
end if
do b = 1, 3
do ipol = 1, 3
do nn2 = 1, kmesh_info%nntot
call comms_scatterv(SAA_q_loc, w*counts(mpirank(comm)), SAA_q_b(:, :, :, nn2, ipol, b), w*counts, w*displs, error, comm)
call fourier_loc_q_to_R(num_kpts, counts, displs, mpirank(comm), wigner_seitz%nrpts, wigner_seitz%irvec, &
kpt_latt, SAA_q_loc, SAA_R_temp)
call comms_reduce(SAA_R_temp(1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts, 'SUM', error, comm)
if (on_root) then
! Apply degeneracy factor and reorder according to the wigner-seitz vectors
call operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, SAA_R_temp, SAA_R_b(:, :, :, ipol, b))
phase2 = -0.5_dp*(wigner_seitz%crvec_pw90(1, :)*kmesh_info%bk(1, nn2, 1) + &
wigner_seitz%crvec_pw90(2, :)*kmesh_info%bk(2, nn2, 1) + &
wigner_seitz%crvec_pw90(3, :)*kmesh_info%bk(3, nn2, 1))
phase2 = exp(cmplx_i*phase2)
SAA_R(:, :, :, ipol, b) = SAA_R(:, :, :, ipol, b) + &
SAA_R_b(:, :, :, ipol, b)*spread(spread(phase2, 1, num_wann), 1, num_wann)
end if
end do
end do
end do
deallocate (SAA_q_loc)
deallocate (SAA_R_temp)
if (on_root) then
deallocate (phase2)
deallocate (SAA_q_b)
deallocate (SAA_R_b)
do b = 1, 3
do ipol = 1, 3
do ir = 1, wigner_seitz%nrpts_pw90
SAA_R(:, :, ir, ipol, b) = SAA_R(:, :, ir, ipol, b) + &
(r0(:, :, b) - 0.5_dp*wigner_seitz%crvec_pw90(b, ir))*SS_R(:, :, ir, ipol)
end do
end do
end do
end if
else
allocate (SAA_q(num_wann, num_wann, num_kpts, 3, 3))
if (on_root) then
SAA_q = sum(SAA_q_b, 4)
deallocate (SAA_q_b)
end if
!
allocate (counts(0:mpisize(comm) - 1))
allocate (displs(0:mpisize(comm) - 1))
w = num_wann*num_wann
call comms_array_split(num_kpts, counts, displs, comm)
allocate (SAA_q_loc(num_wann, num_wann, counts(mpirank(comm))))
allocate (SAA_R_temp(num_wann, num_wann, wigner_seitz%nrpts))
do b = 1, 3
do ipol = 1, 3
call comms_scatterv(SAA_q_loc, w*counts(mpirank(comm)), SAA_q(:, :, :, ipol, b), w*counts, w*displs, error, comm)
call fourier_loc_q_to_R(num_kpts, counts, displs, mpirank(comm), wigner_seitz%nrpts, wigner_seitz%irvec, &
kpt_latt, SAA_q_loc, SAA_R_temp)
call comms_reduce(SAA_R_temp(1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts, 'SUM', error, comm)
if (on_root) then
! Apply degeneracy factor and reorder according to the wigner-seitz vectors
call operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, SAA_R_temp, SAA_R(:, :, :, ipol, b))
end if
end do
end do
deallocate (SAA_q_loc)
deallocate (SAA_q)
deallocate (SAA_R_temp)
end if
call comms_bcast(SAA_R(1, 1, 1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts_pw90*3*3, error, comm)
if (allocated(error)) return
if (print_output%timing_level > 1 .and. print_output%iprint > 0) call io_stopwatch_stop('get_oper: get_SAA_R', timer)
return
deallocate (SAA_R_temp)
113 call set_error_file(error, 'Error: Problem opening input file '//trim(seedname)//'.sIu', comm)
return
114 call set_error_file(error, 'Error: Problem reading input file '//trim(seedname)//'.sIu', comm)
return
end subroutine get_SAA_R
