!-*- 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_get_oper: matrix elements of various operators !
! !
!------------------------------------------------------------!
module w90_get_oper
!================================================
!! Finds the Wannier matrix elements of various operators,
!! starting from k-space matrices generated by an interface
!! (e.g., pw2wannier90) to an ab initio package
!! (e.g., quantum-espresso)
!================================================
use w90_comms, only: comms_bcast, comms_reduce, comms_array_split, comms_scatterv, &
w90_comm_type, mpirank, mpisize
use w90_constants, only: dp, cmplx_0, cmplx_i, cmplx_1, twopi
use w90_io, only: io_stopwatch_start, io_stopwatch_stop
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
public
private :: fourier_q_to_R
private :: get_win_min
integer :: nno, nn1o, nn2o
contains
!================================================!
! PUBLIC PROCEDURES
!================================================!
!================================================
subroutine 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)
!================================================
!
!! computes <0n|H|Rm>, in eV
!! (pwscf uses Ry, but pw2wannier90 converts to eV)
!
!================================================
use w90_postw90_types, only: wigner_seitz_type
use w90_types, only: dis_manifold_type, print_output_type, timer_list_type, &
ws_distance_type, ws_region_type
implicit none
! arguments
type(dis_manifold_type), intent(in) :: dis_manifold
type(print_output_type), intent(in) :: print_output
type(w90_comm_type), intent(in) :: comm
type(wigner_seitz_type), intent(inout) :: wigner_seitz
type(ws_distance_type), intent(in) :: ws_distance
type(ws_region_type), intent(in) :: ws_region
type(timer_list_type), intent(inout) :: timer
type(w90_error_type), allocatable, intent(out) :: error
integer, intent(in) :: num_bands, num_kpts, num_wann, num_valence_bands, stdout
real(kind=dp), intent(in) :: eigval(:, :), real_lattice(3, 3)
real(kind=dp), intent(in) :: scissors_shift
real(kind=dp), intent(in) :: kpt_latt(:, :)
complex(kind=dp), intent(in) :: u_matrix(:, :, :), v_matrix(:, :, :)
complex(kind=dp), allocatable, intent(inout) :: HH_R(:, :, :) ! <0n|r|Rm>
complex(kind=dp), allocatable :: HH_R_temp(:, :, :)
character(len=50), intent(in) :: seedname
logical, intent(in) :: have_disentangled
logical, intent(in) :: effective_model
! local variables
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
complex(kind=dp), allocatable :: sciss_q(:, :, :)
complex(kind=dp), allocatable :: sciss_R(:, :, :)
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_HH_R', timer)
if (wigner_seitz%nrpts < 1) then
call set_error_fatal(error, 'Error: wigner_setiz%nrpts incorrect at get_HH_R() ', comm)
return
end if
allocate (HH_R_temp(num_wann, num_wann, wigner_seitz%nrpts))
if (.not. allocated(HH_R)) then
allocate (HH_R(num_wann, num_wann, wigner_seitz%nrpts_pw90))
else
if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
call io_stopwatch_stop('get_oper: get_HH_R', timer)
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'
open (newunit=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 set_error_fatal(error, 'Error in get_HH_R: orbital indices out of bounds', comm)
return
end if
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.
end if
end if
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
wigner_seitz%irvec(:, ir) = ivdum(:)
if (ivdum(1) == 0 .and. ivdum(2) == 0 .and. ivdum(3) == 0) wigner_seitz%rpt_origin = ir
end if
n = n + 1
end do
close (file_unit)
if (ir /= wigner_seitz%nrpts) then
write (stdout, *) 'ir=', ir, ' nrpts=', wigner_seitz%nrpts
call set_error_fatal(error, 'Error in get_HH_R: inconsistent nrpts values', comm)
return
end if
do ir = 1, wigner_seitz%nrpts
wigner_seitz%crvec(:, ir) = matmul(transpose(real_lattice), wigner_seitz%irvec(:, ir))
end do
wigner_seitz%ndegen(:) = 1 ! This is assumed when reading HH_R from file
!
wigner_seitz%nrpts_pw90 = wigner_seitz%nrpts
wigner_seitz%irvec_pw90 = wigner_seitz%irvec
wigner_seitz%crvec_pw90 = wigner_seitz%crvec
!
! 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_temp, Hopefully the
! result converges (rapidly) with the k-mesh density, but
! one should check
!
if (abs(scissors_shift) > 1.0e-7_dp) then
call set_error_input(error, 'Error in get_HH_R: scissors shift not implemented for ' &
//'effective_model=T', comm)
return
end if
end if
call comms_bcast(HH_R(1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts, error, comm)
if (allocated(error)) return
call comms_bcast(wigner_seitz%ndegen(1), wigner_seitz%nrpts, error, comm)
if (allocated(error)) return
call comms_bcast(wigner_seitz%irvec(1, 1), 3*wigner_seitz%nrpts, error, comm)
if (allocated(error)) return
call comms_bcast(wigner_seitz%crvec(1, 1), 3*wigner_seitz%nrpts, error, comm)
if (allocated(error)) return
if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
call io_stopwatch_stop('get_oper: get_HH_R', timer)
return
end if
! 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) = dis_manifold%ndimwin(ik)
else
num_states(ik) = num_wann
end if
call get_win_min(num_bands, dis_manifold, ik, winmin_q, have_disentangled)
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)
end do
HH_q(m, n, ik) = conjg(HH_q(n, m, ik))
end do
end do
end do
call fourier_q_to_R(num_kpts, wigner_seitz%nrpts, wigner_seitz%irvec, kpt_latt, HH_q, HH_R_temp)
! 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, wigner_seitz%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)
end do
sciss_q(j, i, ik) = conjg(sciss_q(i, j, ik))
end do
end do
end do
call fourier_q_to_R(num_kpts, wigner_seitz%nrpts, wigner_seitz%irvec, kpt_latt, sciss_q, &
sciss_R)
do n = 1, num_wann
sciss_R(n, n, wigner_seitz%rpt_origin) = sciss_R(n, n, wigner_seitz%rpt_origin) + 1.0_dp
end do
sciss_R = sciss_R*scissors_shift
HH_R_temp = HH_R_temp + sciss_R
end if
! Apply degeneracy factor and reorder according to the wigner-seitz vectors
call operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, HH_R_temp, HH_R)
if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
call io_stopwatch_stop('get_oper: get_HH_R', timer)
return
101 call set_error_file(error, 'Error in get_HH_R: problem opening file '// &
trim(seedname)//'_HH_R.dat', comm)
return
end subroutine get_HH_R
!================================================
subroutine get_AA_R_effective(print_output, AA_R, HH_R, nrpts, num_wann, seedname, stdout, timer, error, comm)
!================================================
!
!! AA_a(R) = <0|r_a|R> is the Fourier transform
!! of the Berrry connection AA_a(k) = i<u|del_a u>
!! (a=x,y,z)
!
!================================================
use w90_types, only: print_output_type, timer_list_type, &
ws_distance_type, ws_region_type
implicit none
! arguments
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_wann, nrpts, stdout
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
! local variables
integer :: n, m, i, j, file_unit, ir, io, ivdum(3), ivdum_old(3)
real(kind=dp) :: rdum1_real, rdum1_imag, rdum2_real, rdum2_imag, &
rdum3_real, rdum3_imag
logical :: nn_found
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_AA_R', timer)
if (.not. allocated(AA_R)) then
allocate (AA_R(num_wann, num_wann, nrpts, 3))
else
if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
call io_stopwatch_stop('get_oper: get_AA_R', timer)
return
end if
! Real-space position matrix elements read from file
!
if (.not. allocated(HH_R)) then
call set_error_fatal(error, 'Error in get_AA_R: Must read file'//trim(seedname)//'_HH_R.dat first', comm)
return
end if
AA_R = cmplx_0
if (on_root) then
write (stdout, '(/a)') ' Reading position matrix elements from file ' &
//trim(seedname)//'_AA_R.dat'
open (newunit=file_unit, file=trim(seedname)//'_AA_R.dat', form='formatted', &
status='old', err=103)
read (file_unit, *) ! header
ir = 1
ivdum_old(:) = 0
n = 1
do
read (file_unit, '(5I5,6F12.6)', iostat=io) &
ivdum(1:3), j, i, rdum1_real, rdum1_imag, &
rdum2_real, rdum2_imag, rdum3_real, rdum3_imag
if (io < 0) exit
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 set_error_fatal(error, 'Error in get_AA_R: orbital indices out of bounds', comm)
return
end if
if (n > 1) then
if (ivdum(1) /= ivdum_old(1) .or. ivdum(2) /= ivdum_old(2) .or. &
ivdum(3) /= ivdum_old(3)) ir = ir + 1
end if
ivdum_old = ivdum
AA_R(j, i, ir, 1) = AA_R(j, i, ir, 1) + cmplx(rdum1_real, rdum1_imag, kind=dp)
AA_R(j, i, ir, 2) = AA_R(j, i, ir, 2) + cmplx(rdum2_real, rdum2_imag, kind=dp)
AA_R(j, i, ir, 3) = AA_R(j, i, ir, 3) + cmplx(rdum3_real, rdum3_imag, kind=dp)
n = n + 1
end do
close (file_unit)
! AA_R may not contain the same number of R-vectors as HH_R
! (e.g., if a diagonal representation of the position matrix
! elements is used, but it cannot be larger
if (ir > nrpts) then
write (stdout, *) 'ir=', ir, ' nrpts=', nrpts
call set_error_fatal(error, 'Error in get_AA_R: inconsistent nrpts values', comm)
return
end if
end if
call comms_bcast(AA_R(1, 1, 1, 1), num_wann*num_wann*nrpts*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_AA_R', timer)
return ! careful, don't fall into the below!
101 call set_error_file(error, 'Error: Problem opening input file '//trim(seedname)//'.mmn', comm)
return
102 call set_error_file(error, 'Error: Problem reading input file '//trim(seedname)//'.mmn', comm)
return
103 call set_error_file(error, 'Error in get_AA_R: problem opening file '//trim(seedname)//'_AA_R.dat', comm)
return
end subroutine get_AA_R_effective
!================================================
subroutine get_AA_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, wann_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)
!================================================
!
!! AA_a(R) = <0|r_a|R> is the Fourier transform
!! of the Berrry connection AA_a(k) = i<u|del_a u>
!! (a=x,y,z)
!
!================================================
use w90_postw90_types, only: pw90_berry_mod_type, pw90_oper_read_type, pw90_spin_hall_type, &
wigner_seitz_type
use w90_types, only: dis_manifold_type, kmesh_info_type, print_output_type, timer_list_type, &
ws_distance_type, ws_region_type, wannier_data_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(inout) :: 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(wannier_data_type), intent(in) :: wann_data
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) :: kpt_latt(:, :)
complex(kind=dp), intent(in) :: v_matrix(:, :, :)
complex(kind=dp), allocatable, intent(inout) :: AA_R(:, :, :, :) ! <0n|r|Rm>
complex(kind=dp), allocatable :: AA_R_temp(:, :, :)
complex(kind=dp), allocatable :: AA_R_b(:, :, :, :)
logical, intent(in) :: have_disentangled
character(len=50), intent(in) :: seedname
! local variables
complex(kind=dp), allocatable :: AA_q_b(:, :, :, :, :)
complex(kind=dp), allocatable :: AA_q(:, :, :, :)
complex(kind=dp), allocatable :: AA_q_b_diag(:, :, :)
complex(kind=dp), allocatable :: AA_q_loc(:, :, :)
complex(kind=dp), allocatable :: S_o(:, :)
complex(kind=dp), allocatable :: S(:, :)
integer :: n, m, i, j, &
ik, ik2, ik_prev, nn, inn, nnl, nnm, nnn, &
idir, ncount, nn_count, mmn_in, &
nb_tmp, nkp_tmp, nntot_tmp, file_unit, &
ir, io, w
integer, allocatable :: num_states(:)
integer, allocatable :: counts(:), displs(:)
real(kind=dp) :: m_real, m_imag, rdum1_real, rdum1_imag, &
rdum2_real, rdum2_imag, rdum3_real, rdum3_imag
real(kind=dp), allocatable :: r0(:, :, :)
complex(kind=dp), allocatable :: phase1(:, :), phase2(:)
logical :: nn_found
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_AA_R', timer)
if (.not. allocated(wigner_seitz%wannier_centres_from_AA_R)) then
allocate (wigner_seitz%wannier_centres_from_AA_R(3, num_wann))
end if
if (.not. allocated(AA_R)) then
allocate (AA_R(num_wann, num_wann, wigner_seitz%nrpts_pw90, 3))
else
if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
call io_stopwatch_stop('get_oper: get_AA_R', timer)
return
end if
! Real-space position matrix elements calculated by Fourier
! transforming overlap matrices defined on the ab-initio
! reciprocal mesh
!
! Do everything on root, broadcast AA_R at the end (smaller than S_o)
!
if (on_root) then
allocate (AA_q_b(num_wann, num_wann, num_kpts, kmesh_info%nntot, 3))
AA_R = cmplx_0
else
allocate (AA_q_b(1, 1, 1, kmesh_info%nntot, 3))
end if
!
if (on_root) then
allocate (S_o(num_bands, num_bands))
allocate (S(num_wann, num_wann))
allocate (AA_q_b_diag(num_wann, kmesh_info%nntot, 3))
allocate (num_states(num_kpts))
wigner_seitz%wannier_centres_from_AA_R(:, :) = 0.d0
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=mmn_in, file=trim(seedname)//'.mmn', &
form='formatted', status='old', action='read', err=101)
write (stdout, '(/a)', advance='no') &
' Reading overlaps from '//trim(seedname)//'.mmn in get_AA_R : '
! Read the comment line (header)
read (mmn_in, '(a)', err=102, end=102) header
write (stdout, '(a)') trim(header)
! Read the number of bands, k-points and nearest neighbours
read (mmn_in, *, err=102, end=102) nb_tmp, nkp_tmp, nntot_tmp
! Checks
if (nb_tmp .ne. num_bands) then
call set_error_fatal(error, trim(seedname)//'.mmn has wrong number of bands', comm)
return
end if
if (nkp_tmp .ne. num_kpts) then
call set_error_fatal(error, trim(seedname)//'.mmn has wrong number of k-points', comm)
return
end if
if (nntot_tmp .ne. kmesh_info%nntot) then
call set_error_fatal(error, trim(seedname)//'.mmn has wrong number of nearest neighbours', comm)
return
end if
AA_q_b = cmplx_0
ik_prev = 0
! Composite loop over k-points ik (outer loop) and neighbors ik2 (inner)
do ncount = 1, num_kpts*kmesh_info%nntot
!
!Read from .mmn file the original overlap matrix
! S_o=<u_ik|u_ik2> between ab initio eigenstates
!
read (mmn_in, *, err=102, end=102) ik, ik2, nnl, nnm, nnn
do n = 1, num_bands
do m = 1, num_bands
read (mmn_in, *, err=102, end=102) m_real, m_imag
S_o(m, n) = cmplx(m_real, m_imag, kind=dp)
end do
end do
!debug
!OK
!if(ik.ne.ik_prev .and.ik_prev.ne.0) then
! if(nn_count.ne.nntot)&
! write(stdout,*) 'something wrong in get_AA_R!'
!endif
!enddebug
if (ik .ne. ik_prev) nn_count = 0
nn = 0
nn_found = .false.
do inn = 1, kmesh_info%nntot
if ((ik2 .eq. kmesh_info%nnlist(ik, inn)) .and. &
(nnl .eq. kmesh_info%nncell(1, ik, inn)) .and. &
(nnm .eq. kmesh_info%nncell(2, ik, inn)) .and. &
(nnn .eq. kmesh_info%nncell(3, ik, inn))) then
if (.not. nn_found) then
nn_found = .true.
nn = inn
else
call set_error_fatal(error, 'Error reading '//trim(seedname)//'.mmn.&
& More than one matching nearest neighbour found', comm)
return
end if
end if
end do
if (nn .eq. 0) then
write (stdout, '(/a,i8,2i5,i4,2x,3i3)') ' Error reading '//trim(seedname)//'.mmn:', &
ncount, ik, ik2, nn, nnl, nnm, nnn
call set_error_fatal(error, 'Neighbour not found', comm)
return
end if
nn_count = nn_count + 1 !Check: can also be place after nn=inn (?)
! Wannier-gauge overlap matrix S in the projected subspace
!
call get_gauge_overlap_matrix(num_bands, num_wann, eigval, v_matrix, dis_manifold, &
ik, num_states(ik), kmesh_info%nnlist(ik, nn), &
num_states(kmesh_info%nnlist(ik, nn)), S_o, &
have_disentangled, S)
! save the wannier centers (diagonals of AA_R_temp) to wannier_centres_from_AA_R
! used in pw90common_fourier_R_to_k_new_second_d_TB_conv
do i = 1, num_wann
wigner_seitz%wannier_centres_from_AA_R(:, i) = &
wigner_seitz%wannier_centres_from_AA_R(:, i) &
- kmesh_info%wb(nn)*kmesh_info%bk(:, nn, ik)*aimag(log(S(i, i)))/num_kpts
end do
! Berry connection matrix
!
if (pw90_berry%transl_inv .and. ik .ne. ik_prev) AA_q_b_diag(:, :, :) = cmplx_0
nno = nn
if (pw90_berry%transl_inv_full) nno = kmesh_info%nninv(nn, ik) ! reorder AA_q_b nn indices required for transl_inv_full method
do idir = 1, 3
AA_q_b(:, :, ik, nno, idir) = AA_q_b(:, :, ik, nno, idir) &
+ cmplx_i*kmesh_info%wb(nn)*kmesh_info%bk(idir, nn, ik)*S(:, :)
if (pw90_berry%transl_inv) then
!
! Rewrite band-diagonal elements a la Eq.(31) of MV97
!
do i = 1, num_wann
AA_q_b_diag(i, nno, idir) = AA_q_b_diag(i, nno, idir) &
- kmesh_info%wb(nn)*kmesh_info%bk(idir, nn, ik) &
*aimag(log(S(i, i)))
end do
end if
end do
do idir = 1, 3
if (pw90_berry%transl_inv) then
do n = 1, num_wann
AA_q_b(n, n, ik, nno, idir) = AA_q_b_diag(n, nno, idir)
end do
end if
end do
ik_prev = ik
end do !ncount
close (mmn_in)
if (sum((wigner_seitz%wannier_centres_from_AA_R - wann_data%centres)**2) > 1.0e-08) then
if (pw90_berry%guiding_centres) then
write (stdout, '(/a)', advance='no') &
' Computed and read Wannier centres different. This can happen for guiding_centres=T'
wigner_seitz%wannier_centres_from_AA_R = wann_data%centres
else
call set_error_fatal(error, 'Computed and read Wannier centres different.', comm)
end if
end if
if (pw90_berry%transl_inv_full) then
allocate (r0(num_wann, num_wann, 3))
allocate (phase1(num_wann, num_wann))
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
do nn = 1, kmesh_info%nntot
do ik = 1, num_kpts
phase1 = (r0(:, :, 1)*kmesh_info%bk(1, nn, ik) + &
r0(:, :, 2)*kmesh_info%bk(2, nn, ik) + &
r0(:, :, 3)*kmesh_info%bk(3, nn, ik))
phase1 = exp(cmplx_i*phase1)
nno = kmesh_info%nninv(nn, ik)
do idir = 1, 3
AA_q_b(:, :, ik, nno, idir) = AA_q_b(:, :, ik, nno, idir)*phase1(:, :)
end do
end do ! ik
end do ! nn
deallocate (phase1)
end if
end if !on_root
if (pw90_berry%transl_inv_full) then
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 (AA_q_loc(num_wann, num_wann, counts(mpirank(comm))))
allocate (AA_R_temp(num_wann, num_wann, wigner_seitz%nrpts))
if (on_root) then
allocate (AA_R_b(num_wann, num_wann, wigner_seitz%nrpts_pw90, 3))
allocate (phase2(wigner_seitz%nrpts_pw90))
end if
do idir = 1, 3
do nn = 1, kmesh_info%nntot
call comms_scatterv(AA_q_loc, w*counts(mpirank(comm)), AA_q_b(:, :, :, nn, idir), 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, AA_q_loc, AA_R_temp)
call comms_reduce(AA_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, AA_R_temp, AA_R_b(:, :, :, idir))
phase2 = -0.5_dp*(wigner_seitz%crvec_pw90(1, :)*kmesh_info%bk(1, nn, 1) + &
wigner_seitz%crvec_pw90(2, :)*kmesh_info%bk(2, nn, 1) + &
wigner_seitz%crvec_pw90(3, :)*kmesh_info%bk(3, nn, 1))
phase2 = exp(cmplx_i*phase2)
AA_R(:, :, :, idir) = AA_R(:, :, :, idir) + AA_R_b(:, :, :, idir)*spread(spread(phase2, 1, num_wann), 1, num_wann)
end if
end do
end do
deallocate (AA_q_loc)
deallocate (AA_R_temp)
if (on_root) then
deallocate (phase2)
deallocate (AA_q_b)
deallocate (AA_R_b)
do ir = 1, wigner_seitz%nrpts_pw90
if ((wigner_seitz%irvec_pw90(1, ir) .eq. 0) .and. &
(wigner_seitz%irvec_pw90(2, ir) .eq. 0) .and. &
(wigner_seitz%irvec_pw90(3, ir) .eq. 0)) then
do i = 1, num_wann
AA_R(i, i, ir, :) = wigner_seitz%wannier_centres_from_AA_R(:, i)
end do
exit
end if
end do
end if
else
allocate (AA_q(num_wann, num_wann, num_kpts, 3))
if (on_root) then
AA_q = sum(AA_q_b, 4)
deallocate (AA_q_b)
! Since Eq.(44) WYSV06 does not preserve the Hermiticity of the
! Berry potential matrix, take Hermitean part (whether this
! makes a difference or not for e.g. the AHC, depends on which
! expression is used to evaluate the Berry curvature.
! See comments in berry_wanint.F90)
!
do idir = 1, 3
do ik = 1, num_kpts
AA_q(:, :, ik, idir) = &
0.5_dp*(AA_q(:, :, ik, idir) &
+ conjg(transpose(AA_q(:, :, ik, idir))))
end do
end do
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 (AA_q_loc(num_wann, num_wann, counts(mpirank(comm))))
allocate (AA_R_temp(num_wann, num_wann, wigner_seitz%nrpts))
do idir = 1, 3
call comms_scatterv(AA_q_loc, w*counts(mpirank(comm)), AA_q(:, :, :, idir), 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, AA_q_loc, AA_R_temp)
call comms_reduce(AA_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, AA_R_temp, AA_R(:, :, :, idir))
end if
end do
deallocate (AA_q_loc)
deallocate (AA_q)
deallocate (AA_R_temp)
end if
call comms_bcast(AA_R(1, 1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts_pw90*3, error, comm)
call comms_bcast(wigner_seitz%wannier_centres_from_AA_R(1, 1), num_wann*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_AA_R', timer)
return
101 call set_error_file(error, 'Error: Problem opening input file '//trim(seedname)//'.mmn', comm)
return
102 call set_error_file(error, 'Error: Problem reading input file '//trim(seedname)//'.mmn', comm)
return
103 call set_error_file(error, 'Error in get_AA_R: problem opening file '//trim(seedname)//'_AA_R.dat', comm)
return
end subroutine get_AA_R
!================================================
subroutine get_BB_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, HH_R, BB_R, v_matrix, eigval, &
scissors_shift, wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, num_wann, &
have_disentangled, seedname, stdout, timer, error, comm)
!================================================
!
!! BB_a(R)=<0n|H(r-R)|Rm> is the Fourier transform of
!! BB_a(k) = i<u|H|del_a u> (a=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(inout) :: HH_R(:, :, :)
complex(kind=dp), allocatable, intent(inout) :: BB_R(:, :, :, :) ! <0|H(r-R)|R>
complex(kind=dp), allocatable :: BB_R_temp(:, :, :)
complex(kind=dp), allocatable :: BB_R_b(:, :, :, :)
logical, intent(in) :: have_disentangled
character(len=50), intent(in) :: seedname
! local variables
integer :: idir, n, m, nn, i, j, ir, &
ik, ik2, inn, nnl, nnm, nnn, &
winmin_q, winmin_qb, ncount, &
nb_tmp, nkp_tmp, nntot_tmp, mmn_in, w
complex(kind=dp), allocatable :: S_o(:, :)
complex(kind=dp), allocatable :: BB_q(:, :, :, :)
complex(kind=dp), allocatable :: BB_q_loc(:, :, :)
complex(kind=dp), allocatable :: BB_q_b(:, :, :, :, :)
complex(kind=dp), allocatable :: H_q_qb(:, :)
real(kind=dp), allocatable :: r0(:, :, :)
complex(kind=dp), allocatable :: phase1(:, :), phase2(:)
integer, allocatable :: num_states(:)
integer, allocatable :: counts(:), displs(:)
real(kind=dp) :: m_real, m_imag
logical :: nn_found
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_BB_R', timer)
if (.not. allocated(BB_R)) then
allocate (BB_R(num_wann, num_wann, wigner_seitz%nrpts_pw90, 3))
else
if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
call io_stopwatch_stop('get_oper: get_BB_R', timer)
return
end if
if (on_root) then
allocate (BB_q_b(num_wann, num_wann, num_kpts, kmesh_info%nntot, 3))
BB_R = cmplx_0
else
allocate (BB_q_b(1, 1, 1, kmesh_info%nntot, 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 BB_R', comm)
return
end if
allocate (S_o(num_bands, num_bands))
allocate (H_q_qb(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=mmn_in, file=trim(seedname)//'.mmn', &
form='formatted', status='old', action='read', err=103)
write (stdout, '(/a)', advance='no') &
' Reading overlaps from '//trim(seedname)//'.mmn in get_BB_R : '
! Read the comment line (header)
read (mmn_in, '(a)', err=104, end=104) header
write (stdout, '(a)') trim(header)
! Read the number of bands, k-points and nearest neighbours
read (mmn_in, *, err=104, end=104) nb_tmp, nkp_tmp, nntot_tmp
! Checks
if (nb_tmp .ne. num_bands) then
call set_error_fatal(error, trim(seedname)//'.mmn has wrong number of bands', comm)
return
end if
if (nkp_tmp .ne. num_kpts) then
call set_error_fatal(error, trim(seedname)//'.mmn has wrong number of k-points', comm)
return
end if
if (nntot_tmp .ne. kmesh_info%nntot) then
call set_error_fatal(error, trim(seedname)//'.mmn has wrong number of nearest neighbours', comm)
return
end if
BB_q_b = cmplx_0
do ncount = 1, num_kpts*kmesh_info%nntot
!
!Read from .mmn file the original overlap matrix
! S_o=<u_ik|u_ik2> between ab initio eigenstates
!
read (mmn_in, *, err=104, end=104) ik, ik2, nnl, nnm, nnn
do n = 1, num_bands
do m = 1, num_bands
read (mmn_in, *, err=104, end=104) m_real, m_imag
S_o(m, n) = cmplx(m_real, m_imag, kind=dp)
end do
end do
nn = 0
nn_found = .false.
do inn = 1, kmesh_info%nntot
if ((ik2 .eq. kmesh_info%nnlist(ik, inn)) .and. &
(nnl .eq. kmesh_info%nncell(1, ik, inn)) .and. &
(nnm .eq. kmesh_info%nncell(2, ik, inn)) .and. &
(nnn .eq. kmesh_info%nncell(3, ik, inn))) then
if (.not. nn_found) then
nn_found = .true.
nn = inn
else
call set_error_fatal(error, 'Error reading '//trim(seedname)//'.mmn.&
& More than one matching nearest neighbour found', comm)
return
end if
end if
end do
if (nn .eq. 0) then
write (stdout, '(/a,i8,2i5,i4,2x,3i3)') ' Error reading '//trim(seedname)//'.mmn:', &
ncount, ik, ik2, nn, nnl, nnm, nnn
call set_error_fatal(error, 'Neighbour not found', comm)
return
end if
call get_win_min(num_bands, dis_manifold, ik, winmin_q, have_disentangled)
call get_win_min(num_bands, dis_manifold, kmesh_info%nnlist(ik, nn), winmin_qb, &
have_disentangled)
call get_gauge_overlap_matrix(num_bands, num_wann, eigval, v_matrix, dis_manifold, &
ik, num_states(ik), kmesh_info%nnlist(ik, nn), &
num_states(kmesh_info%nnlist(ik, nn)), S_o, &
have_disentangled, H=H_q_qb)
if (pw90_berry%transl_inv_full) then
phase1 = (r0(:, :, 1)*kmesh_info%bk(1, nn, ik) + &
r0(:, :, 2)*kmesh_info%bk(2, nn, ik) + &
r0(:, :, 3)*kmesh_info%bk(3, nn, ik))
phase1 = exp(cmplx_i*phase1)
else
phase1 = cmplx_1
end if
do idir = 1, 3
nno = kmesh_info%nninv(nn, ik)
BB_q_b(:, :, ik, nno, idir) = BB_q_b(:, :, ik, nno, idir) &
+ cmplx_i*phase1(:, :)*kmesh_info%wb(nn)*kmesh_info%bk(idir, nn, ik) &
*H_q_qb(:, :)
end do
end do !ncount
close (mmn_in)
deallocate (phase1)
end if !on_root
if (pw90_berry%transl_inv_full) then
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 (BB_q_loc(num_wann, num_wann, counts(mpirank(comm))))
allocate (BB_R_temp(num_wann, num_wann, wigner_seitz%nrpts))
if (on_root) then
allocate (BB_R_b(num_wann, num_wann, wigner_seitz%nrpts_pw90, 3))
allocate (phase2(wigner_seitz%nrpts_pw90))
end if
do idir = 1, 3
do nn = 1, kmesh_info%nntot
call comms_scatterv(BB_q_loc, w*counts(mpirank(comm)), BB_q_b(:, :, :, nn, idir), 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, BB_q_loc, BB_R_temp)
call comms_reduce(BB_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, BB_R_temp, BB_R_b(:, :, :, idir))
phase2 = -0.5_dp*(wigner_seitz%crvec_pw90(1, :)*kmesh_info%bk(1, nn, 1) + &
wigner_seitz%crvec_pw90(2, :)*kmesh_info%bk(2, nn, 1) + &
wigner_seitz%crvec_pw90(3, :)*kmesh_info%bk(3, nn, 1))
phase2 = exp(cmplx_i*phase2)
BB_R(:, :, :, idir) = BB_R(:, :, :, idir) + BB_R_b(:, :, :, idir)*spread(spread(phase2, 1, num_wann), 1, num_wann)
end if
end do
end do
deallocate (BB_q_loc)
deallocate (BB_R_temp)
if (on_root) then
deallocate (phase2)
deallocate (BB_q_b)
deallocate (BB_R_b)
do idir = 1, 3
do ir = 1, wigner_seitz%nrpts_pw90
BB_R(:, :, ir, idir) = BB_R(:, :, ir, idir) + &
(r0(:, :, idir) - 0.5_dp*wigner_seitz%crvec_pw90(idir, ir))*HH_R(:, :, ir)
end do
end do
end if
else
allocate (BB_q(num_wann, num_wann, num_kpts, 3))
if (on_root) then
BB_q = sum(BB_q_b, 4)
deallocate (BB_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 (BB_q_loc(num_wann, num_wann, counts(mpirank(comm))))
allocate (BB_R_temp(num_wann, num_wann, wigner_seitz%nrpts))
do idir = 1, 3
call comms_scatterv(BB_q_loc, w*counts(mpirank(comm)), BB_q(:, :, :, idir), 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, BB_q_loc, BB_R_temp)
call comms_reduce(BB_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, BB_R_temp, BB_R(:, :, :, idir))
end if
end do
deallocate (BB_q_loc)
deallocate (BB_q)
deallocate (BB_R_temp)
end if
call comms_bcast(BB_R(1, 1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts_pw90*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_BB_R', timer)
return
103 call set_error_file(error, 'Error: Problem opening input file '//trim(seedname)//'.mmn', comm)
return
104 call set_error_file(error, 'Error: Problem reading input file '//trim(seedname)//'.mmn', comm)
return
end subroutine get_BB_R
!================================================
subroutine get_CC_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, pw90_oper_read, HH_R, &
BB_R, CC_R, v_matrix, eigval, scissors_shift, wigner_seitz, ws_distance, ws_region, &
num_bands, num_kpts, num_wann, have_disentangled, seedname, stdout, timer, error, comm)
!================================================
!
!! 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_postw90_types, only: pw90_berry_mod_type, pw90_oper_read_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
use w90_utility, only: utility_compar
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(pw90_oper_read_type), intent(in) :: pw90_oper_read
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(inout) :: HH_R(:, :, :)
complex(kind=dp), allocatable, intent(inout) :: BB_R(:, :, :, :)
complex(kind=dp), allocatable, intent(inout) :: CC_R(:, :, :, :, :) ! <0|r_alpha.H(r-R)_beta|R>
complex(kind=dp), allocatable :: CC_R_temp(:, :, :)
complex(kind=dp), allocatable :: CC_R_b(:, :, :, :, :)
logical, intent(in) :: have_disentangled
character(len=50), intent(in) :: seedname
! local variables
integer :: m, n, a, b, nn1, nn2, ik, nb_tmp, nkp_tmp, i, j, ir, ir2, &
nntot_tmp, uHu_in, qb1, qb2, winmin_qb1, winmin_qb2, &
ifpos, ifneg, w
integer, allocatable :: num_states(:)
integer, allocatable :: counts(:), displs(:)
complex(kind=dp), allocatable :: CC_q(:, :, :, :, :)
complex(kind=dp), allocatable :: CC_q_b(:, :, :, :, :, :, :)
complex(kind=dp), allocatable :: CC_q_loc(:, :, :)
complex(kind=dp), allocatable :: Ho_qb1_q_qb2(:, :)
complex(kind=dp), allocatable :: H_qb1_q_qb2(:, :)
real(kind=dp), allocatable :: r0(:, :, :)
complex(kind=dp), allocatable :: phase1(:, :), phase2(:)
real(kind=dp) :: c_real, c_img
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_CC_R', timer)
if (.not. allocated(CC_R)) then
allocate (CC_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_CC_R', timer)
return
end if
if (on_root) then
allocate (CC_q_b(num_wann, num_wann, num_kpts, kmesh_info%nntot, kmesh_info%nntot, 3, 3))
CC_R = cmplx_0
else
allocate (CC_q_b(1, 1, 1, kmesh_info%nntot, 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 CC_R', comm)
return
end if
allocate (Ho_qb1_q_qb2(num_bands, num_bands))
allocate (H_qb1_q_qb2(num_wann, num_wann))
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
allocate (num_states(num_kpts))
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
if (pw90_oper_read%uHu_formatted) then
open (newunit=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 (newunit=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
end if
if (nb_tmp .ne. num_bands) then
call set_error_fatal(error, trim(seedname)//'.uHu 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)//'.uHu 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)//'.uHu has not the right number of nearest neighbours', comm)
return
end if
CC_q_b = cmplx_0
do ik = 1, num_kpts
do nn2 = 1, kmesh_info%nntot
qb2 = kmesh_info%nnlist(ik, nn2)
call get_win_min(num_bands, dis_manifold, qb2, winmin_qb2, have_disentangled)
do nn1 = 1, kmesh_info%nntot
qb1 = kmesh_info%nnlist(ik, nn1)
call get_win_min(num_bands, dis_manifold, qb1, winmin_qb1, have_disentangled)
!
! Read from .uHu file the matrices <u_{q+b1}|H_q|u_{q+b2}>
! between the original ab initio eigenstates
!
if (pw90_oper_read%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)
end if
! 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(num_bands, num_wann, eigval, v_matrix, dis_manifold, &
qb1, num_states(qb1), qb2, num_states(qb2), &
Ho_qb1_q_qb2, have_disentangled, H_qb1_q_qb2)
if (pw90_berry%transl_inv_full) then
phase1 = -(r0(:, :, 1)*kmesh_info%bk(1, nn1, ik) + &
r0(:, :, 2)*kmesh_info%bk(2, nn1, ik) + &
r0(:, :, 3)*kmesh_info%bk(3, nn1, ik)) &
+ (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
do b = 1, 3
do a = 1, b
nn1o = kmesh_info%nninv(nn1, ik)
nn2o = kmesh_info%nninv(nn2, ik)
CC_q_b(:, :, ik, nn1o, nn2o, a, b) = CC_q_b(:, :, ik, nn1o, nn2o, a, b) &
+ phase1(:, :)*kmesh_info%wb(nn1)*kmesh_info%bk(a, nn1, ik) &
*kmesh_info%wb(nn2)*kmesh_info%bk(b, nn2, ik)*H_qb1_q_qb2(:, :)
end do
end do
end do !nn1
end do !nn2
end do !ik
close (uHu_in)
deallocate (phase1)
end if !on_root
if (pw90_berry%transl_inv_full) then
if (.not. allocated(HH_R)) then
call set_error_fatal(error, 'transl_inv_full=T for CC_R needs HH_R', comm)
end if
if (.not. allocated(BB_R)) then
call set_error_fatal(error, 'transl_inv_full=T for CC_R needs BB_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 (CC_q_loc(num_wann, num_wann, counts(mpirank(comm))))
allocate (CC_R_temp(num_wann, num_wann, wigner_seitz%nrpts))
if (on_root) then
allocate (CC_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 a = 1, 3
do nn2 = 1, kmesh_info%nntot
do nn1 = 1, kmesh_info%nntot
call comms_scatterv(CC_q_loc, w*counts(mpirank(comm)), CC_q_b(:, :, :, nn1, nn2, a, 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, CC_q_loc, CC_R_temp)
call comms_reduce(CC_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, &
CC_R_temp, CC_R_b(:, :, :, a, b))
phase2 = -0.5_dp*(wigner_seitz%crvec_pw90(1, :)*kmesh_info%bk(1, nn1, 1) + &
wigner_seitz%crvec_pw90(2, :)*kmesh_info%bk(2, nn1, 1) + &
wigner_seitz%crvec_pw90(3, :)*kmesh_info%bk(3, nn1, 1)) &
- 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)
CC_R(:, :, :, a, b) = CC_R(:, :, :, a, b) + CC_R_b(:, :, :, a, b)* &
spread(spread(phase2, 1, num_wann), 1, num_wann)
end if
end do
end do
end do
end do
deallocate (CC_q_loc)
deallocate (CC_R_temp)
if (on_root) then
deallocate (phase2)
deallocate (CC_q_b)
deallocate (CC_R_b)
do b = 1, 3
do a = 1, 3
do ir = 1, wigner_seitz%nrpts_pw90
CC_R(:, :, ir, a, b) = CC_R(:, :, ir, a, b) + &
(r0(:, :, a) + 0.5_dp*wigner_seitz%crvec_pw90(a, ir))* &
BB_R(:, :, ir, b)
do ir2 = 1, wigner_seitz%nrpts_pw90
call utility_compar(wigner_seitz%crvec_pw90(1, ir), &
wigner_seitz%crvec_pw90(1, ir2), ifpos, ifneg)
if (ifneg .eq. 1) then
CC_R(:, :, ir, a, b) = CC_R(:, :, ir, a, b) + &
conjg(transpose(BB_R(:, :, ir2, a)))* &
(r0(:, :, b) - 0.5_dp*wigner_seitz%crvec_pw90(b, ir))
exit
end if
end do
CC_R(:, :, ir, a, b) = CC_R(:, :, ir, a, b) + &
(r0(:, :, a) + 0.5_dp*wigner_seitz%crvec_pw90(a, ir))* &
wigner_seitz%crvec_pw90(b, ir)*HH_R(:, :, ir)
end do
end do
end do
end if
else
allocate (CC_q(num_wann, num_wann, num_kpts, 3, 3))
if (on_root) then
CC_q = sum(sum(CC_q_b, 5), 4)
deallocate (CC_q_b)
!
do b = 1, 3
do a = 1, b
do ik = 1, num_kpts
CC_q(:, :, ik, b, a) = conjg(transpose(CC_q(:, :, ik, a, b)))
end do
end do
end do
!
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 (CC_q_loc(num_wann, num_wann, counts(mpirank(comm))))
allocate (CC_R_temp(num_wann, num_wann, wigner_seitz%nrpts))
do b = 1, 3
do a = 1, 3
call comms_scatterv(CC_q_loc, w*counts(mpirank(comm)), CC_q(:, :, :, a, 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, CC_q_loc, CC_R_temp)
call comms_reduce(CC_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, CC_R_temp, CC_R(:, :, :, a, b))
end if
end do
end do
deallocate (CC_q_loc)
deallocate (CC_q)
deallocate (CC_R_temp)
end if
call comms_bcast(CC_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_CC_R', timer)
return
105 call set_error_file(error, 'Error: Problem opening input file '//trim(seedname)//'.uHu', comm)
return
106 call set_error_file(error, 'Error: Problem reading input file '//trim(seedname)//'.uHu', comm)
return
end subroutine get_CC_R
!================================================
subroutine get_FF_R(num_bands, num_kpts, num_wann, wigner_seitz, ws_distance, ws_region, v_matrix, &
FF_R, dis_manifold, kmesh_info, kpt_latt, print_output, have_disentangled, stdout, &
seedname, timer, error, comm)
!================================================
!
!! FF_ab(R) = <0|r_a.(r-R)_b|R> is the Fourier transform of
!! FF_ab(k) = <del_a u|del_b u> (a=alpha,b=beta)
!
!================================================
use w90_postw90_types, only: 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(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) :: kpt_latt(:, :)
complex(kind=dp), intent(in) :: v_matrix(:, :, :)
complex(kind=dp), allocatable, intent(inout) :: FF_R(:, :, :, :, :) ! <0|r_alpha.(r-R)_beta|R>
complex(kind=dp), allocatable :: FF_R_temp(:, :, :, :, :)
character(len=50), intent(in) :: seedname
logical, intent(in) :: have_disentangled
! local variables
integer :: i, j, ii, jj, m, n, a, b, nn1, nn2, ik, nb_tmp, nkp_tmp, nntot_tmp, &
uIu_in, qb1, qb2, winmin_qb1, winmin_qb2
integer, allocatable :: num_states(:)
complex(kind=dp), allocatable :: FF_q(:, :, :, :, :)
complex(kind=dp), allocatable :: Lo_qb1_q_qb2(:, :)
complex(kind=dp), allocatable :: L_qb1_q_qb2(:, :)
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_FF_R', timer)
allocate (FF_R_temp(num_wann, num_wann, wigner_seitz%nrpts, 3, 3))
if (.not. allocated(FF_R)) then
allocate (FF_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_FF_R', timer)
return
end if
if (on_root) then
allocate (Lo_qb1_q_qb2(num_bands, num_bands))
allocate (L_qb1_q_qb2(num_wann, num_wann))
allocate (FF_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) = dis_manifold%ndimwin(ik)
else
num_states(ik) = num_wann
end if
end do
open (newunit=uIu_in, file=TRIM(seedname)//".uIu", form='unformatted', &
status='old', action='read', err=107)
write (stdout, '(/a)', advance='no') &
' Reading uIu overlaps from '//trim(seedname)//'.uIu in get_FF_R: '
read (uIu_in, err=108, end=108) header
write (stdout, '(a)') trim(header)
read (uIu_in, err=108, end=108) nb_tmp, nkp_tmp, nntot_tmp
if (nb_tmp .ne. num_bands) then
call set_error_fatal(error, trim(seedname)//'.uIu 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)//'.uIu 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)//'.uIu has not the right number of nearest neighbours', comm)
return
end if
FF_q = cmplx_0
do ik = 1, num_kpts
do nn2 = 1, kmesh_info%nntot
qb2 = kmesh_info%nnlist(ik, nn2)
call get_win_min(num_bands, dis_manifold, qb2, winmin_qb2, have_disentangled)
do nn1 = 1, kmesh_info%nntot
qb1 = kmesh_info%nnlist(ik, nn1)
call get_win_min(num_bands, dis_manifold, qb1, winmin_qb1, have_disentangled)
!
! Read from .uIu file the matrices <u_{q+b1}|u_{q+b2}>
! between the original ab initio eigenstates
!
! do m=1,num_bands
! do n=1,num_bands
! read(uIu_in,err=108,end=108) Lo_qb1_q_qb2(m,n)
! end do
! end do
read (uIu_in, err=108, end=108) &
((Lo_qb1_q_qb2(n, m), n=1, num_bands), m=1, num_bands)
!
! **************************************************************
! 2013-08-09: Do we need to take a transpose here?! SEE get_CC_R
Lo_qb1_q_qb2 = transpose(Lo_qb1_q_qb2) ! added 2013-08-09 (?)
! **************************************************************
!
! Transform to projected subspace, Wannier gauge
!
L_qb1_q_qb2(:, :) = cmplx_0
do m = 1, num_wann
do n = 1, num_wann
do i = 1, num_states(qb1)
ii = winmin_qb1 + i - 1
do j = 1, num_states(qb2)
jj = winmin_qb2 + j - 1
L_qb1_q_qb2(n, m) = L_qb1_q_qb2(n, m) &
+ conjg(v_matrix(i, n, qb1)) &
*Lo_qb1_q_qb2(ii, jj) &
*v_matrix(j, m, qb2)
end do
end do
end do
end do
do b = 1, 3
do a = 1, b
FF_q(:, :, ik, a, b) = FF_q(:, :, ik, a, b) &
+ kmesh_info%wb(nn1)*kmesh_info%bk(a, nn1, ik) &
*kmesh_info%wb(nn2)*kmesh_info%bk(b, nn2, ik)*L_qb1_q_qb2(:, :)
end do
end do
end do !nn1
end do !nn2
do b = 1, 3
do a = 1, b
FF_q(:, :, ik, b, a) = conjg(transpose(FF_q(:, :, ik, a, b)))
end do
end do
end do !ik
close (uIu_in)
do b = 1, 3
do a = 1, 3
call fourier_q_to_R(num_kpts, wigner_seitz%nrpts, wigner_seitz%irvec, kpt_latt, &
FF_q(:, :, :, a, b), FF_R_temp(:, :, :, a, b))
end do
end do
do b = 1, 3
do a = 1, 3
call operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, &
FF_R_temp(:, :, :, a, b), FF_R(:, :, :, a, b))
end do
end do
end if !on_root
call comms_bcast(FF_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_FF_R', timer)
return
deallocate (FF_R_temp)
107 call set_error_file(error, 'Error: Problem opening input file '//trim(seedname)//'.uIu', comm)
return
108 call set_error_file(error, 'Error: Problem reading input file '//trim(seedname)//'.uIu', comm)
return
end subroutine get_FF_R
!================================================
subroutine get_SS_R(dis_manifold, kpt_latt, print_output, pw90_oper_read, SS_R, v_matrix, &
eigval, wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, num_wann, &
have_disentangled, seedname, stdout, timer, error, comm)
!================================================
!
!! Wannier representation of the Pauli matrices: <0n|sigma_a|Rm>
!! (a=x,y,z)
!
!================================================
use w90_postw90_types, only: pw90_oper_read_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(dis_manifold_type), intent(in) :: dis_manifold
type(pw90_oper_read_type), intent(in) :: pw90_oper_read
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) :: stdout, num_bands, num_kpts, num_wann
real(kind=dp), intent(in) :: eigval(:, :)
real(kind=dp), intent(in) :: kpt_latt(:, :)
complex(kind=dp), intent(in) :: v_matrix(:, :, :)
complex(kind=dp), allocatable, intent(inout) :: SS_R(:, :, :, :) ! <0n|sigma_x,y,z|Rm>
complex(kind=dp), allocatable :: SS_R_temp(:, :, :, :)
character(len=50), intent(in) :: seedname
logical, intent(in) :: have_disentangled
! local variables
complex(kind=dp), allocatable :: spn_o(:, :, :, :), SS_q(:, :, :, :), spn_temp(:, :)
real(kind=dp) :: s_real, s_img
integer, allocatable :: num_states(:)
integer :: m, n, spn_in, ik, is, &
nb_tmp, nkp_tmp, ierr, s, counter
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_SS_R', timer)
allocate (SS_R_temp(num_wann, num_wann, wigner_seitz%nrpts, 3))
if (.not. allocated(SS_R)) then
allocate (SS_R(num_wann, num_wann, wigner_seitz%nrpts_pw90, 3))
else
return ! been here before
end if
if (on_root) then
allocate (spn_o(num_bands, num_bands, num_kpts, 3))
allocate (SS_q(num_wann, num_wann, num_kpts, 3))
allocate (num_states(num_kpts))
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
! Read from .spn file the original spin matrices <psi_nk|sigma_i|psi_mk>
! (sigma_i = Pauli matrix) between ab initio eigenstates
!
if (pw90_oper_read%spn_formatted) then
open (newunit=spn_in, file=trim(seedname)//'.spn', form='formatted', &
status='old', err=109)
write (stdout, '(/a)', advance='no') &
' Reading spin matrices from '//trim(seedname)//'.spn in get_SS_R : '
read (spn_in, *, err=110, end=110) header
write (stdout, '(a)') trim(header)
read (spn_in, *, err=110, end=110) nb_tmp, nkp_tmp
else
open (newunit=spn_in, file=trim(seedname)//'.spn', form='unformatted', &
status='old', err=109)
write (stdout, '(/a)', advance='no') &
' Reading spin matrices from '//trim(seedname)//'.spn in get_SS_R : '
read (spn_in, err=110, end=110) header
write (stdout, '(a)') trim(header)
read (spn_in, err=110, end=110) nb_tmp, nkp_tmp
end if
if (nb_tmp .ne. num_bands) then
call set_error_fatal(error, trim(seedname)//'.spn has wrong number of bands', comm)
return
end if
if (nkp_tmp .ne. num_kpts) then
call set_error_fatal(error, trim(seedname)//'.spn has wrong number of k-points', comm)
return
end if
if (pw90_oper_read%spn_formatted) then
do ik = 1, num_kpts
do m = 1, num_bands
do n = 1, m
read (spn_in, *, err=110, end=110) s_real, s_img
spn_o(n, m, ik, 1) = cmplx(s_real, s_img, dp)
read (spn_in, *, err=110, end=110) s_real, s_img
spn_o(n, m, ik, 2) = cmplx(s_real, s_img, dp)
read (spn_in, *, err=110, end=110) s_real, s_img
spn_o(n, m, ik, 3) = cmplx(s_real, s_img, dp)
! Read upper-triangular part, now build the rest
spn_o(m, n, ik, 1) = conjg(spn_o(n, m, ik, 1))
spn_o(m, n, ik, 2) = conjg(spn_o(n, m, ik, 2))
spn_o(m, n, ik, 3) = conjg(spn_o(n, m, ik, 3))
end do
end do
end do
else
allocate (spn_temp(3, (num_bands*(num_bands + 1))/2), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating spm_temp in get_SS_R', comm)
return
end if
do ik = 1, num_kpts
read (spn_in) ((spn_temp(s, m), s=1, 3), m=1, (num_bands*(num_bands + 1))/2)
counter = 0
do m = 1, num_bands
do n = 1, m
counter = counter + 1
spn_o(n, m, ik, 1) = spn_temp(1, counter)
spn_o(m, n, ik, 1) = conjg(spn_temp(1, counter))
spn_o(n, m, ik, 2) = spn_temp(2, counter)
spn_o(m, n, ik, 2) = conjg(spn_temp(2, counter))
spn_o(n, m, ik, 3) = spn_temp(3, counter)
spn_o(m, n, ik, 3) = conjg(spn_temp(3, counter))
end do
end do
end do
deallocate (spn_temp, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating spm_temp in get_SS_R', comm)
return
end if
end if
close (spn_in)
! Transform to projected subspace, Wannier gauge
!
SS_q(:, :, :, :) = cmplx_0
do ik = 1, num_kpts
do is = 1, 3
call get_gauge_overlap_matrix(num_bands, num_wann, eigval, v_matrix, dis_manifold, &
ik, num_states(ik), ik, num_states(ik), &
spn_o(:, :, ik, is), have_disentangled, SS_q(:, :, ik, is))
end do !is
end do !ik
call fourier_q_to_R(num_kpts, wigner_seitz%nrpts, wigner_seitz%irvec, kpt_latt, SS_q(:, :, :, 1), SS_R_temp(:, :, :, 1))
call fourier_q_to_R(num_kpts, wigner_seitz%nrpts, wigner_seitz%irvec, kpt_latt, SS_q(:, :, :, 2), SS_R_temp(:, :, :, 2))
call fourier_q_to_R(num_kpts, wigner_seitz%nrpts, wigner_seitz%irvec, kpt_latt, SS_q(:, :, :, 3), SS_R_temp(:, :, :, 3))
call operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, SS_R_temp(:, :, :, 1), SS_R(:, :, :, 1))
call operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, SS_R_temp(:, :, :, 2), SS_R(:, :, :, 2))
call operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, SS_R_temp(:, :, :, 3), SS_R(:, :, :, 3))
end if !on_root
call comms_bcast(SS_R(1, 1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts_pw90*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_SS_R', timer)
return
deallocate (SS_R_temp)
109 call set_error_file(error, 'Error: Problem opening input file '//trim(seedname)//'.spn', comm)
return
110 call set_error_file(error, 'Error: Problem reading input file '//trim(seedname)//'.spn', comm)
return
end subroutine get_SS_R
!================================================
subroutine get_SHC_R(dis_manifold, kmesh_info, kpt_latt, print_output, pw90_oper_read, &
pw90_spin_hall, SH_R, SHR_R, SR_R, v_matrix, eigval, scissors_shift, &
wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, num_wann, &
num_valence_bands, have_disentangled, seedname, stdout, timer, error, comm)
!================================================
!
!! Compute several matrices for spin Hall conductivity
!! SR_R = <0n|sigma_{x,y,z}.(r-R)_alpha|Rm>
!! SHR_R = <0n|sigma_{x,y,z}.H.(r-R)_alpha|Rm>
!! SH_R = <0n|sigma_{x,y,z}.H|Rm>
!
!================================================
use w90_postw90_types, only: pw90_oper_read_type, pw90_spin_hall_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(dis_manifold_type), intent(in) :: dis_manifold
type(kmesh_info_type), intent(in) :: kmesh_info
type(pw90_oper_read_type), intent(in) :: pw90_oper_read
type(print_output_type), intent(in) :: print_output
type(pw90_spin_hall_type), intent(in) :: pw90_spin_hall
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(timer_list_type), intent(inout) :: timer
type(w90_comm_type), intent(in) :: comm
type(w90_error_type), allocatable, intent(out) :: error
integer, intent(in) :: stdout, num_bands, num_kpts, num_wann, num_valence_bands
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(inout) :: SR_R(:, :, :, :, :) ! <0n|sigma_x,y,z.(r-R)_alpha|Rm>
complex(kind=dp), allocatable, intent(inout) :: SHR_R(:, :, :, :, :) ! <0n|sigma_x,y,z.H.(r-R)_alpha|Rm>
complex(kind=dp), allocatable, intent(inout) :: SH_R(:, :, :, :) ! <0n|sigma_x,y,z.H|Rm>
complex(kind=dp), allocatable :: SR_R_temp(:, :, :, :, :)
complex(kind=dp), allocatable :: SHR_R_temp(:, :, :, :, :)
complex(kind=dp), allocatable :: SH_R_temp(:, :, :, :)
character(len=50), intent(in) :: seedname
logical, intent(in) :: have_disentangled
! local variables
complex(kind=dp), allocatable :: SR_q(:, :, :, :, :)
complex(kind=dp), allocatable :: SHR_q(:, :, :, :, :)
complex(kind=dp), allocatable :: SH_q(:, :, :, :)
complex(kind=dp), allocatable :: S_o(:, :)
complex(kind=dp), allocatable :: spn_o(:, :, :, :), spn_temp(:, :)
complex(kind=dp), allocatable :: H_o(:, :, :)
complex(kind=dp), allocatable :: SH_o(:, :, :, :)
complex(kind=dp) :: SM_o(num_bands, num_bands, 3)
complex(kind=dp) :: SHM_o(num_bands, num_bands, 3)
complex(kind=dp) :: SS_q(num_wann, num_wann, 3)
complex(kind=dp) :: SM_q(num_wann, num_wann, 3)
complex(kind=dp) :: SHM_q(num_wann, num_wann, 3)
real(kind=dp) :: s_real, s_img
integer :: spn_in, counter, ierr, s, is
integer :: n, m, &
ik, ik2, ik_prev, nn, inn, nnl, nnm, nnn, &
idir, ncount, nn_count, mmn_in, &
nb_tmp, nkp_tmp, nntot_tmp
integer, allocatable :: num_states(:)
real(kind=dp) :: m_real, m_imag
logical :: nn_found
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_SHC_R', timer)
allocate (SR_R_temp(num_wann, num_wann, wigner_seitz%nrpts, 3, 3))
allocate (SHR_R_temp(num_wann, num_wann, wigner_seitz%nrpts, 3, 3))
allocate (SH_R_temp(num_wann, num_wann, wigner_seitz%nrpts, 3))
if (.not. allocated(SR_R)) then
allocate (SR_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_SHC_R', timer)
return
end if
if (.not. allocated(SHR_R)) then
allocate (SHR_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_SHC_R', timer)
return
end if
if (.not. allocated(SH_R)) then
allocate (SH_R(num_wann, num_wann, wigner_seitz%nrpts_pw90, 3))
else
if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
call io_stopwatch_stop('get_oper: get_SHC_R', timer)
return
end if
! start copying from get_SS_R, Junfeng Qiao
! read spn file
if (on_root) then
allocate (spn_o(num_bands, num_bands, num_kpts, 3))
allocate (num_states(num_kpts))
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
! Read from .spn file the original spin matrices <psi_nk|sigma_i|psi_mk>
! (sigma_i = Pauli matrix) between ab initio eigenstates
!
if (pw90_oper_read%spn_formatted) then
open (newunit=spn_in, file=trim(seedname)//'.spn', form='formatted', &
status='old', err=109)
write (stdout, '(/a)', advance='no') &
' Reading spin matrices from '//trim(seedname)//'.spn in get_SHC_R : '
read (spn_in, *, err=110, end=110) header
write (stdout, '(a)') trim(header)
read (spn_in, *, err=110, end=110) nb_tmp, nkp_tmp
else
open (newunit=spn_in, file=trim(seedname)//'.spn', form='unformatted', &
status='old', err=109)
write (stdout, '(/a)', advance='no') &
' Reading spin matrices from '//trim(seedname)//'.spn in get_SHC_R : '
read (spn_in, err=110, end=110) header
write (stdout, '(a)') trim(header)
read (spn_in, err=110, end=110) nb_tmp, nkp_tmp
end if
if (nb_tmp .ne. num_bands) then
call set_error_fatal(error, trim(seedname)//'.spn has wrong number of bands', comm)
return
end if
if (nkp_tmp .ne. num_kpts) then
call set_error_fatal(error, trim(seedname)//'.spn has wrong number of k-points', comm)
return
end if
if (pw90_oper_read%spn_formatted) then
do ik = 1, num_kpts
do m = 1, num_bands
do n = 1, m
read (spn_in, *, err=110, end=110) s_real, s_img
spn_o(n, m, ik, 1) = cmplx(s_real, s_img, dp)
read (spn_in, *, err=110, end=110) s_real, s_img
spn_o(n, m, ik, 2) = cmplx(s_real, s_img, dp)
read (spn_in, *, err=110, end=110) s_real, s_img
spn_o(n, m, ik, 3) = cmplx(s_real, s_img, dp)
! Read upper-triangular part, now build the rest
spn_o(m, n, ik, 1) = conjg(spn_o(n, m, ik, 1))
spn_o(m, n, ik, 2) = conjg(spn_o(n, m, ik, 2))
spn_o(m, n, ik, 3) = conjg(spn_o(n, m, ik, 3))
end do
end do
end do
else
allocate (spn_temp(3, (num_bands*(num_bands + 1))/2), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating spm_temp in get_SHC_R', comm)
return
end if
do ik = 1, num_kpts
read (spn_in) ((spn_temp(s, m), s=1, 3), m=1, (num_bands*(num_bands + 1))/2)
counter = 0
do m = 1, num_bands
do n = 1, m
counter = counter + 1
spn_o(n, m, ik, 1) = spn_temp(1, counter)
spn_o(m, n, ik, 1) = conjg(spn_temp(1, counter))
spn_o(n, m, ik, 2) = spn_temp(2, counter)
spn_o(m, n, ik, 2) = conjg(spn_temp(2, counter))
spn_o(n, m, ik, 3) = spn_temp(3, counter)
spn_o(m, n, ik, 3) = conjg(spn_temp(3, counter))
end do
end do
end do
deallocate (spn_temp, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating spm_temp in get_SHC_R', comm)
return
end if
end if
close (spn_in)
end if !on_root
! end copying from get_SS_R, Junfeng Qiao
! start copying from get_HH_R, Junfeng Qiao
! Note this is different from get_HH_R, at here we need the
! original Hamiltonian to construct SHR_R, SH_R.
if (on_root) then
allocate (H_o(num_bands, num_bands, num_kpts))
H_o = cmplx_0
do ik = 1, num_kpts
do m = 1, num_bands
H_o(m, m, ik) = eigval(m, ik)
end do
! scissors shift applied to the original Hamiltonian
if (num_valence_bands > 0 .and. abs(scissors_shift) > 1.0e-7_dp) then
do m = num_valence_bands + 1, num_bands
H_o(m, m, ik) = H_o(m, m, ik) + scissors_shift
end do
else if (pw90_spin_hall%bandshift) then
do m = pw90_spin_hall%bandshift_firstband, num_bands
H_o(m, m, ik) = H_o(m, m, ik) + pw90_spin_hall%bandshift_energyshift
end do
end if
end do
end if !on_root
! end copying from get_HH_R, Junfeng Qiao
! start copying from get_AA_R, Junfeng Qiao
! read mmn file
!
if (on_root) then
allocate (SR_q(num_wann, num_wann, num_kpts, 3, 3))
allocate (SHR_q(num_wann, num_wann, num_kpts, 3, 3))
allocate (SH_q(num_wann, num_wann, num_kpts, 3))
allocate (S_o(num_bands, num_bands))
open (newunit=mmn_in, file=trim(seedname)//'.mmn', &
form='formatted', status='old', action='read', err=101)
write (stdout, '(/a)', advance='no') &
' Reading overlaps from '//trim(seedname)//'.mmn in get_SHC_R : '
! Read the comment line (header)
read (mmn_in, '(a)', err=102, end=102) header
write (stdout, '(a)') trim(header)
! Read the number of bands, k-points and nearest neighbours
read (mmn_in, *, err=102, end=102) nb_tmp, nkp_tmp, nntot_tmp
! Checks
if (nb_tmp .ne. num_bands) then
call set_error_fatal(error, trim(seedname)//'.mmn has wrong number of bands', comm)
return
end if
if (nkp_tmp .ne. num_kpts) then
call set_error_fatal(error, trim(seedname)//'.mmn has wrong number of k-points', comm)
return
end if
if (nntot_tmp .ne. kmesh_info%nntot) then
call set_error_fatal(error, trim(seedname)//'.mmn has wrong number of nearest neighbours', comm)
return
end if
SR_q = cmplx_0
SHR_q = cmplx_0
SH_q = cmplx_0
ik_prev = 0
! QZYZ18 Eq.(48)
allocate (SH_o(num_bands, num_bands, num_kpts, 3))
SH_o = cmplx_0
do ik = 1, num_kpts
do is = 1, 3
SH_o(:, :, ik, is) = matmul(spn_o(:, :, ik, is), H_o(:, :, ik))
call get_gauge_overlap_matrix(num_bands, num_wann, eigval, v_matrix, dis_manifold, &
ik, num_states(ik), ik, num_states(ik), &
SH_o(:, :, ik, is), have_disentangled, SH_q(:, :, ik, is))
end do
end do
! Composite loop over k-points ik (outer loop) and neighbors ik2 (inner)
do ncount = 1, num_kpts*kmesh_info%nntot
!
!Read from .mmn file the original overlap matrix
! S_o=<u_ik|u_ik2> between ab initio eigenstates
!
read (mmn_in, *, err=102, end=102) ik, ik2, nnl, nnm, nnn
do n = 1, num_bands
do m = 1, num_bands
read (mmn_in, *, err=102, end=102) m_real, m_imag
S_o(m, n) = cmplx(m_real, m_imag, kind=dp)
end do
end do
!debug
!OK
!if(ik.ne.ik_prev .and.ik_prev.ne.0) then
! if(nn_count.ne.nntot)&
! write(stdout,*) 'something wrong in get_AA_R!'
!endif
!enddebug
if (ik .ne. ik_prev) nn_count = 0
nn = 0
nn_found = .false.
do inn = 1, kmesh_info%nntot
if ((ik2 .eq. kmesh_info%nnlist(ik, inn)) .and. &
(nnl .eq. kmesh_info%nncell(1, ik, inn)) .and. &
(nnm .eq. kmesh_info%nncell(2, ik, inn)) .and. &
(nnn .eq. kmesh_info%nncell(3, ik, inn))) then
if (.not. nn_found) then
nn_found = .true.
nn = inn
else
call set_error_fatal(error, 'Error reading '//trim(seedname)//'.mmn.&
& More than one matching nearest neighbour found', comm)
return
end if
end if
end do
if (nn .eq. 0) then
write (stdout, '(/a,i8,2i5,i4,2x,3i3)') ' Error reading '//trim(seedname)//'.mmn:', &
ncount, ik, ik2, nn, nnl, nnm, nnn
call set_error_fatal(error, 'Neighbour not found', comm)
return
end if
nn_count = nn_count + 1 !Check: can also be place after nn=inn (?)
SM_o = cmplx_0
SHM_o = cmplx_0
SS_q = cmplx_0
SM_q = cmplx_0
SHM_q = cmplx_0
do is = 1, 3
! QZYZ18 Eq.(50)
SM_o(:, :, is) = matmul(spn_o(:, :, ik, is), S_o(:, :))
! QZYZ18 Eq.(51)
SHM_o(:, :, is) = matmul(SH_o(:, :, ik, is), S_o(:, :))
! Transform to projected subspace, Wannier gauge
!
! QZYZ18 Eq.(50)
call get_gauge_overlap_matrix(num_bands, num_wann, eigval, v_matrix, dis_manifold, &
ik, num_states(ik), ik, num_states(ik), &
spn_o(:, :, ik, is), have_disentangled, SS_q(:, :, is))
! QZYZ18 Eq.(50)
call get_gauge_overlap_matrix(num_bands, num_wann, eigval, v_matrix, dis_manifold, &
ik, num_states(ik), kmesh_info%nnlist(ik, nn), &
num_states(kmesh_info%nnlist(ik, nn)), SM_o(:, :, is), &
have_disentangled, SM_q(:, :, is))
! QZYZ18 Eq.(51)
call get_gauge_overlap_matrix(num_bands, num_wann, eigval, v_matrix, dis_manifold, &
ik, num_states(ik), kmesh_info%nnlist(ik, nn), &
num_states(kmesh_info%nnlist(ik, nn)), SHM_o(:, :, is), &
have_disentangled, SHM_q(:, :, is))
! Assuming all neighbors of a given point are read in sequence!
!
do idir = 1, 3
! QZYZ18 Eq.(50)
SR_q(:, :, ik, is, idir) = SR_q(:, :, ik, is, idir) &
+ kmesh_info%wb(nn)*kmesh_info%bk(idir, nn, ik) &
*(SM_q(:, :, is) - SS_q(:, :, is))
! QZYZ18 Eq.(51)
SHR_q(:, :, ik, is, idir) = SHR_q(:, :, ik, is, idir) &
+ kmesh_info%wb(nn)*kmesh_info%bk(idir, nn, ik) &
*(SHM_q(:, :, is) - SH_q(:, :, ik, is))
end do
end do
ik_prev = ik
end do !ncount
close (mmn_in)
do is = 1, 3
! QZYZ18 Eq.(46)
call fourier_q_to_R(num_kpts, wigner_seitz%nrpts, wigner_seitz%irvec, kpt_latt, &
SH_q(:, :, :, is), SH_R_temp(:, :, :, is))
do idir = 1, 3
! QZYZ18 Eq.(44)
call fourier_q_to_R(num_kpts, wigner_seitz%nrpts, wigner_seitz%irvec, kpt_latt, &
SR_q(:, :, :, is, idir), SR_R_temp(:, :, :, is, idir))
! QZYZ18 Eq.(45)
call fourier_q_to_R(num_kpts, wigner_seitz%nrpts, wigner_seitz%irvec, kpt_latt, &
SHR_q(:, :, :, is, idir), SHR_R_temp(:, :, :, is, idir))
end do
end do
do is = 1, 3
! QZYZ18 Eq.(46)
call operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, &
SH_R_temp(:, :, :, is), SH_R(:, :, :, is))
do idir = 1, 3
! QZYZ18 Eq.(44)
call operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, &
SR_R_temp(:, :, :, is, idir), SR_R(:, :, :, is, idir))
! QZYZ18 Eq.(45)
call operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, &
SHR_R_temp(:, :, :, is, idir), SHR_R(:, :, :, is, idir))
end do
end do
SR_R = cmplx_i*SR_R
SHR_R = cmplx_i*SHR_R
end if !on_root
call comms_bcast(SH_R(1, 1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts_pw90*3, error, comm)
if (allocated(error)) return
call comms_bcast(SR_R(1, 1, 1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts_pw90*3*3, error, comm)
if (allocated(error)) return
call comms_bcast(SHR_R(1, 1, 1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts_pw90*3*3, error, comm)
if (allocated(error)) return
! end copying from get_AA_R, Junfeng Qiao
if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
call io_stopwatch_stop('get_oper: get_SHC_R', timer)
return
deallocate (SH_R_temp)
deallocate (SR_R_temp)
deallocate (SHR_R_temp)
101 call set_error_file(error, 'Error: Problem opening input file '//trim(seedname)//'.mmn', comm)
return
102 call set_error_file(error, 'Error: Problem reading input file '//trim(seedname)//'.mmn', comm)
return
109 call set_error_file(error, 'Error: Problem opening input file '//trim(seedname)//'.spn', comm)
return
110 call set_error_file(error, 'Error: Problem reading input file '//trim(seedname)//'.spn', comm)
return
end subroutine get_SHC_R
!================================================
subroutine get_SH_R(dis_manifold, kmesh_info, kpt_latt, print_output, pw90_oper_read, &
pw90_spin_hall, SH_R, v_matrix, eigval, scissors_shift, &
wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, num_wann, &
num_valence_bands, have_disentangled, seedname, stdout, timer, error, comm)
!================================================
!
!! Compute several matrices for spin Hall conductivity
!! SH_R = <0n|sigma_{x,y,z}.H|Rm>
!
!================================================
use w90_postw90_types, only: pw90_oper_read_type, pw90_spin_hall_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(dis_manifold_type), intent(in) :: dis_manifold
type(kmesh_info_type), intent(in) :: kmesh_info
type(pw90_oper_read_type), intent(in) :: pw90_oper_read
type(print_output_type), intent(in) :: print_output
type(pw90_spin_hall_type), intent(in) :: pw90_spin_hall
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(timer_list_type), intent(inout) :: timer
type(w90_comm_type), intent(in) :: comm
type(w90_error_type), allocatable, intent(out) :: error
integer, intent(in) :: stdout, num_bands, num_kpts, num_wann, num_valence_bands
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(inout) :: SH_R(:, :, :, :) ! <0n|sigma_x,y,z.H|Rm>
complex(kind=dp), allocatable :: SH_R_temp(:, :, :, :)
character(len=50), intent(in) :: seedname
logical, intent(in) :: have_disentangled
! local variables
complex(kind=dp), allocatable :: SH_q(:, :, :, :)
complex(kind=dp), allocatable :: S_o(:, :)
complex(kind=dp), allocatable :: spn_o(:, :, :, :), spn_temp(:, :)
complex(kind=dp), allocatable :: H_o(:, :, :)
complex(kind=dp), allocatable :: SH_o(:, :, :, :)
real(kind=dp) :: s_real, s_img
integer :: spn_in, counter, ierr, s, is
integer :: n, m, ik, idir, nb_tmp, nkp_tmp
integer, allocatable :: num_states(:)
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_SH_R', timer)
allocate (SH_R_temp(num_wann, num_wann, wigner_seitz%nrpts, 3))
if (.not. allocated(SH_R)) then
allocate (SH_R(num_wann, num_wann, wigner_seitz%nrpts_pw90, 3))
else
if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
call io_stopwatch_stop('get_oper: get_SH_R', timer)
return
end if
! start copying from get_SS_R, Junfeng Qiao
! read spn file
if (on_root) then
allocate (spn_o(num_bands, num_bands, num_kpts, 3))
allocate (num_states(num_kpts))
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
! Read from .spn file the original spin matrices <psi_nk|sigma_i|psi_mk>
! (sigma_i = Pauli matrix) between ab initio eigenstates
!
if (pw90_oper_read%spn_formatted) then
open (newunit=spn_in, file=trim(seedname)//'.spn', form='formatted', &
status='old', err=109)
write (stdout, '(/a)', advance='no') &
' Reading spin matrices from '//trim(seedname)//'.spn in get_SH_R : '
read (spn_in, *, err=110, end=110) header
write (stdout, '(a)') trim(header)
read (spn_in, *, err=110, end=110) nb_tmp, nkp_tmp
else
open (newunit=spn_in, file=trim(seedname)//'.spn', form='unformatted', &
status='old', err=109)
write (stdout, '(/a)', advance='no') &
' Reading spin matrices from '//trim(seedname)//'.spn in get_SH_R : '
read (spn_in, err=110, end=110) header
write (stdout, '(a)') trim(header)
read (spn_in, err=110, end=110) nb_tmp, nkp_tmp
end if
if (nb_tmp .ne. num_bands) then
call set_error_fatal(error, trim(seedname)//'.spn has wrong number of bands', comm)
return
end if
if (nkp_tmp .ne. num_kpts) then
call set_error_fatal(error, trim(seedname)//'.spn has wrong number of k-points', comm)
return
end if
if (pw90_oper_read%spn_formatted) then
do ik = 1, num_kpts
do m = 1, num_bands
do n = 1, m
read (spn_in, *, err=110, end=110) s_real, s_img
spn_o(n, m, ik, 1) = cmplx(s_real, s_img, dp)
read (spn_in, *, err=110, end=110) s_real, s_img
spn_o(n, m, ik, 2) = cmplx(s_real, s_img, dp)
read (spn_in, *, err=110, end=110) s_real, s_img
spn_o(n, m, ik, 3) = cmplx(s_real, s_img, dp)
! Read upper-triangular part, now build the rest
spn_o(m, n, ik, 1) = conjg(spn_o(n, m, ik, 1))
spn_o(m, n, ik, 2) = conjg(spn_o(n, m, ik, 2))
spn_o(m, n, ik, 3) = conjg(spn_o(n, m, ik, 3))
end do
end do
end do
else
allocate (spn_temp(3, (num_bands*(num_bands + 1))/2), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating spm_temp in get_SH_R', comm)
return
end if
do ik = 1, num_kpts
read (spn_in) ((spn_temp(s, m), s=1, 3), m=1, (num_bands*(num_bands + 1))/2)
counter = 0
do m = 1, num_bands
do n = 1, m
counter = counter + 1
spn_o(n, m, ik, 1) = spn_temp(1, counter)
spn_o(m, n, ik, 1) = conjg(spn_temp(1, counter))
spn_o(n, m, ik, 2) = spn_temp(2, counter)
spn_o(m, n, ik, 2) = conjg(spn_temp(2, counter))
spn_o(n, m, ik, 3) = spn_temp(3, counter)
spn_o(m, n, ik, 3) = conjg(spn_temp(3, counter))
end do
end do
end do
deallocate (spn_temp, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating spm_temp in get_SH_R', comm)
return
end if
end if
close (spn_in)
end if !on_root
! end copying from get_SS_R, Junfeng Qiao
! start copying from get_HH_R, Junfeng Qiao
! Note this is different from get_HH_R, at here we need the
! original Hamiltonian to construct SHR_R, SH_R.
if (on_root) then
allocate (H_o(num_bands, num_bands, num_kpts))
H_o = cmplx_0
do ik = 1, num_kpts
do m = 1, num_bands
H_o(m, m, ik) = eigval(m, ik)
end do
! scissors shift applied to the original Hamiltonian
if (num_valence_bands > 0 .and. abs(scissors_shift) > 1.0e-7_dp) then
do m = num_valence_bands + 1, num_bands
H_o(m, m, ik) = H_o(m, m, ik) + scissors_shift
end do
else if (pw90_spin_hall%bandshift) then
do m = pw90_spin_hall%bandshift_firstband, num_bands
H_o(m, m, ik) = H_o(m, m, ik) + pw90_spin_hall%bandshift_energyshift
end do
end if
end do
end if !on_root
! end copying from get_HH_R, Junfeng Qiao
! start copying from get_AA_R, Junfeng Qiao
! read mmn file
!
if (on_root) then
allocate (SH_q(num_wann, num_wann, num_kpts, 3))
SH_q = cmplx_0
! QZYZ18 Eq.(48)
allocate (SH_o(num_bands, num_bands, num_kpts, 3))
SH_o = cmplx_0
do ik = 1, num_kpts
do is = 1, 3
SH_o(:, :, ik, is) = matmul(spn_o(:, :, ik, is), H_o(:, :, ik))
call get_gauge_overlap_matrix(num_bands, num_wann, eigval, v_matrix, dis_manifold, &
ik, num_states(ik), ik, num_states(ik), &
SH_o(:, :, ik, is), have_disentangled, SH_q(:, :, ik, is))
end do
end do
do is = 1, 3
! QZYZ18 Eq.(46)
call fourier_q_to_R(num_kpts, wigner_seitz%nrpts, wigner_seitz%irvec, kpt_latt, &
SH_q(:, :, :, is), SH_R_temp(:, :, :, is))
end do
do is = 1, 3
! QZYZ18 Eq.(46)
call operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, &
SH_R_temp(:, :, :, is), SH_R(:, :, :, is))
end do
end if !on_root
call comms_bcast(SH_R(1, 1, 1, 1), num_wann*num_wann*wigner_seitz%nrpts_pw90*3, error, comm)
if (allocated(error)) return
! end copying from get_AA_R, Junfeng Qiao
if (print_output%timing_level > 1 .and. print_output%iprint > 0) &
call io_stopwatch_stop('get_oper: get_SH_R', timer)
return
deallocate (SH_R_temp)
101 call set_error_file(error, 'Error: Problem opening input file '//trim(seedname)//'.mmn', comm)
return
102 call set_error_file(error, 'Error: Problem reading input file '//trim(seedname)//'.mmn', comm)
return
109 call set_error_file(error, 'Error: Problem opening input file '//trim(seedname)//'.spn', comm)
return
110 call set_error_file(error, 'Error: Problem reading input file '//trim(seedname)//'.spn', comm)
return
end subroutine get_SH_R
!================================================
subroutine get_SBB_R(pw90_berry, dis_manifold, kmesh_info, kpt_latt, print_output, SH_R, SBB_R, &
v_matrix, scissors_shift, wigner_seitz, ws_distance, ws_region, num_bands, &
num_kpts, num_wann, have_disentangled, seedname, stdout, timer, error, comm)
!================================================!
!
! SBB_ab(R) = <0|s_a.H.(r-R)_b|R> is the Fourier transform of
! SBB_ab(k) = <u|s_a.H|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) :: SH_R(:, :, :, :)
complex(kind=dp), allocatable, intent(inout) :: SBB_R(:, :, :, :, :) ! <0n|sigma_x,y,z.H.(r-R)_alpha|Rm>
complex(kind=dp), allocatable :: SBB_R_temp(:, :, :)
complex(kind=dp), allocatable :: SBB_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, sHu_in, qb2, winmin_q, winmin_qb2, w
integer :: ipol
integer, allocatable :: num_states(:)
integer, allocatable :: counts(:), displs(:)
complex(kind=dp), allocatable :: SBB_q_b(:, :, :, :, :, :)
complex(kind=dp), allocatable :: SBB_q(:, :, :, :, :)
complex(kind=dp), allocatable :: SBB_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_SBB_R', timer)
if (.not. allocated(SBB_R)) then
allocate (SBB_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_SBB_R', timer)
return
end if
if (on_root) then
allocate (SBB_q_b(num_wann, num_wann, num_kpts, kmesh_info%nntot, 3, 3))
SBB_R = cmplx_0
else
allocate (SBB_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 SBB_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=sHu_in, file=trim(seedname)//".sHu", form='unformatted', &
status='old', action='read', err=111)
write (stdout, '(/a)', advance='no') &
' Reading sHu overlaps from '//trim(seedname)//'.sHu in get_SBB_R: '
read (sHu_in, err=112, end=112) header
write (stdout, '(a)') trim(header)
read (sHu_in, err=112, end=112) nb_tmp, nkp_tmp, nntot_tmp
if (nb_tmp .ne. num_bands) then
call set_error_fatal(error, trim(seedname)//'.sHu 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)//'.sHu 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)//'.sHu has not the right number of nearest neighbours', comm)
return
end if
SBB_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 .sHu file the matrices <u_q|s_a H_q|u_{q+b2}>
! between the original ab initio eigenstates
!
read (sHu_in, err=112, end=112) &
((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)
SBB_q_b(:, :, ik, nn2o, ipol, b) = SBB_q_b(:, :, ik, nn2o, ipol, b) + &
cmplx_i*phase1(:, :)*kmesh_info%wb(nn2)*kmesh_info%bk(b, nn2, ik)*H_q_qb2(:, :)
end do
end do !ipol
end do !nn2
end do !ik
close (sHu_in)
deallocate (phase1)
end if !on_root
if (pw90_berry%transl_inv_full) then
if (.not. allocated(SH_R)) then
call set_error_fatal(error, 'transl_inv_full=T for SBB_R needs SH_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 (SBB_q_loc(num_wann, num_wann, counts(mpirank(comm))))
allocate (SBB_R_temp(num_wann, num_wann, wigner_seitz%nrpts))
if (on_root) then
allocate (SBB_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(SBB_q_loc, w*counts(mpirank(comm)), SBB_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, SBB_q_loc, SBB_R_temp)
call comms_reduce(SBB_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, SBB_R_temp, SBB_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)
SBB_R(:, :, :, ipol, b) = SBB_R(:, :, :, ipol, b) + &
SBB_R_b(:, :, :, ipol, b)*spread(spread(phase2, 1, num_wann), 1, num_wann)
end if
end do
end do
end do
deallocate (SBB_q_loc)
deallocate (SBB_R_temp)
if (on_root) then
deallocate (phase2)
deallocate (SBB_q_b)
deallocate (SBB_R_b)
do b = 1, 3
do ipol = 1, 3
do ir = 1, wigner_seitz%nrpts_pw90
SBB_R(:, :, ir, ipol, b) = SBB_R(:, :, ir, ipol, b) + &
(r0(:, :, b) - 0.5_dp*wigner_seitz%crvec_pw90(b, ir))*SH_R(:, :, ir, ipol)
end do
end do
end do
end if
else
allocate (SBB_q(num_wann, num_wann, num_kpts, 3, 3))
if (on_root) then
SBB_q = sum(SBB_q_b, 4)
deallocate (SBB_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 (SBB_q_loc(num_wann, num_wann, counts(mpirank(comm))))
allocate (SBB_R_temp(num_wann, num_wann, wigner_seitz%nrpts))
do b = 1, 3
do ipol = 1, 3
call comms_scatterv(SBB_q_loc, w*counts(mpirank(comm)), SBB_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, SBB_q_loc, SBB_R_temp)
call comms_reduce(SBB_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, SBB_R_temp, SBB_R(:, :, :, ipol, b))
end if
end do
end do
deallocate (SBB_q_loc)
deallocate (SBB_q)
deallocate (SBB_R_temp)
end if
call comms_bcast(SBB_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_SBB_R', timer)
return
deallocate (SBB_R_temp)
111 call set_error_file(error, 'Error: Problem opening input file '//trim(seedname)//'.sHu', comm)
return
112 call set_error_file(error, 'Error: Problem reading input file '//trim(seedname)//'.sHu', comm)
return
end subroutine get_SBB_R
!================================================
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
!================================================!
! PRIVATE PROCEDURES
!================================================!
!================================================!
subroutine fourier_q_to_R(num_kpts, nrpts, irvec, kpt_latt, op_q, op_R)
!================================================
!
!! Fourier transforms Wannier-gauge representation
!! of a given operator O from q-space to R-space:
!!
!! O_ij(q) --> O_ij(R) = (1/N_kpts) sum_q e^{-iqR} O_ij(q)
!
!================================================
implicit none
! Arguments
real(kind=dp), intent(in) :: kpt_latt(:, :)
integer, intent(in) :: num_kpts, nrpts, irvec(:, :)
complex(kind=dp), intent(in) :: op_q(:, :, :) !! Operator in q-space
complex(kind=dp), intent(out) :: op_R(:, :, :) !! Operator in R-space
! local variables
integer :: ir, ik
real(kind=dp) :: rdotq
complex(kind=dp) :: phase_fac
op_R = cmplx_0
do ir = 1, nrpts
do ik = 1, num_kpts
rdotq = twopi*dot_product(kpt_latt(:, ik), irvec(:, ir))
phase_fac = exp(-cmplx_i*rdotq)
op_R(:, :, ir) = op_R(:, :, ir) + phase_fac*op_q(:, :, ik)
end do
end do
op_R = op_R/real(num_kpts, dp)
end subroutine fourier_q_to_R
!================================================!
subroutine fourier_loc_q_to_R(num_kpts, counts, displs, rank, nrpts, irvec, kpt_latt, op_q, op_R)
!================================================
!
!! Fourier transforms Wannier-gauge representation
!! of a given operator O from q-space to R-space:
!!
!! O_ij(q) --> O_ij(R) = (1/N_kpts) sum_q e^{-iqR} O_ij(q)
!
!================================================
implicit none
! Arguments
real(kind=dp), intent(in) :: kpt_latt(:, :)
integer, intent(in) :: num_kpts, rank, counts(0:), displs(0:), nrpts, irvec(:, :)
complex(kind=dp), intent(in) :: op_q(:, :, :) !! Operator in q-space
complex(kind=dp), intent(out) :: op_R(:, :, :) !! Operator in R-space
! local variables
integer :: ir, ik, ik_start, ik_end
real(kind=dp) :: rdotq
complex(kind=dp) :: phase_fac
ik_start = displs(rank) + 1
ik_end = displs(rank) + counts(rank)
op_R = cmplx_0
do ir = 1, nrpts
do ik = ik_start, ik_end
rdotq = twopi*dot_product(kpt_latt(:, ik), irvec(:, ir))
phase_fac = exp(-cmplx_i*rdotq)
op_R(:, :, ir) = op_R(:, :, ir) + phase_fac*op_q(:, :, ik - ik_start + 1)
end do
end do
op_R = op_R/real(num_kpts, dp)
end subroutine fourier_loc_q_to_R
!================================================
subroutine get_win_min(num_bands, dis_manifold, ik, win_min, have_disentangled)
!================================================
!
!! Find the lower bound (band index) of the
!! outer energy window at the specified k-point
!
!================================================
use w90_types, only: dis_manifold_type, print_output_type
implicit none
! Arguments
integer, intent(in) :: num_bands, ik !! Index of the required k-point
integer, intent(out) :: win_min !! Index of the lower band of the outer energy window
type(dis_manifold_type), intent(in) :: dis_manifold
logical, intent(in) :: have_disentangled
! local variables
integer :: j
if (.not. have_disentangled) then
win_min = 1
return
end if
do j = 1, num_bands
if (dis_manifold%lwindow(j, ik)) then
win_min = j
exit
end if
end do
end subroutine get_win_min
!================================================
subroutine get_gauge_overlap_matrix(num_bands, num_wann, eigval, v_matrix, dis_manifold, ik_a, &
ns_a, ik_b, ns_b, S_o, have_disentangled, S, H)
!================================================
!
! Wannier-gauge overlap matrix S in the projected subspace
!
! TODO: Update this documentation of this routine and
! possibliy give it a better name. The routine has been
! generalized multiple times.
!
!================================================
use w90_types, only: dis_manifold_type
use w90_utility, only: utility_zgemmm
implicit none
! arguments
type(dis_manifold_type), intent(in) :: dis_manifold
real(kind=dp), intent(in) :: eigval(:, :)
complex(kind=dp), intent(in) :: S_o(:, :), v_matrix(:, :, :)
integer, intent(in) :: num_wann, num_bands, ik_a, ns_a, ik_b, ns_b
logical, intent(in) :: have_disentangled
complex(kind=dp), intent(out), optional :: S(:, :), H(:, :)
! local variables
integer :: wm_a, wm_b
call get_win_min(num_bands, dis_manifold, ik_a, wm_a, have_disentangled)
call get_win_min(num_bands, dis_manifold, ik_b, wm_b, have_disentangled)
call utility_zgemmm(v_matrix(1:ns_a, 1:num_wann, ik_a), 'C', &
S_o(wm_a:wm_a + ns_a - 1, wm_b:wm_b + ns_b - 1), 'N', &
v_matrix(1:ns_b, 1:num_wann, ik_b), 'N', &
S, eigval(wm_a:wm_a + ns_a - 1, ik_a), H)
end subroutine get_gauge_overlap_matrix
!============================================================================
subroutine operator_wigner_setup(ws_distance, ws_region, wigner_seitz, num_wann, op_R, op_R_opt_ws)
!==========================================================================
!
! Also, divide real-space matrix elements with the degeneracy factor.
! For use_ws_distance = true, reorder the real-space grid index
! using ir_ind_ws_to_pw90.
!
! After this routine, irvec_pw90, crvec_pw90, and nrpts_pw90 can be
! used in the fourier_R_to_k routines, irrespective of use_ws_distance.
!
!==========================================================================
use w90_constants, only: dp, cmplx_0
use w90_types, only: ws_region_type, ws_distance_type
use w90_postw90_types, only: wigner_seitz_type
type(ws_distance_type), intent(in) :: ws_distance
type(ws_region_type), intent(in) :: ws_region
type(wigner_seitz_type), intent(in) :: wigner_seitz
integer, intent(in) :: num_wann
complex(kind=dp), intent(in) :: op_R(num_wann, num_wann, wigner_seitz%nrpts)
!! operator in real-space grid, before applying ndegen
complex(kind=dp), intent(inout) :: op_R_opt_ws(num_wann, num_wann, wigner_seitz%nrpts_pw90)
!! operator in real-space grid, after applying ndegen
integer :: ir, jr, i, j, ideg
op_R_opt_ws = cmplx_0
if (ws_region%use_ws_distance) then
do ir = 1, wigner_seitz%nrpts
do j = 1, num_wann
do i = 1, num_wann
do ideg = 1, ws_distance%ndeg(i, j, ir)
jr = wigner_seitz%ir_ind_ws_to_pw90(ideg, i, j, ir)
op_R_opt_ws(i, j, jr) = op_R_opt_ws(i, j, jr) &
+ op_R(i, j, ir)/real(wigner_seitz%ndegen(ir)* &
ws_distance%ndeg(i, j, ir), dp)
end do
end do
end do
end do
else ! .not. use_ws_distance
! Note that nrpts_pw90 == nrpts if use_ws_distance == .false.
do ir = 1, wigner_seitz%nrpts
op_R_opt_ws(:, :, ir) = op_R(:, :, ir)/real(wigner_seitz%ndegen(ir), dp)
end do
end if ! use_ws_distance
end subroutine operator_wigner_setup
end module w90_get_oper
