AA_a(R) = <0|r_a|R> is the Fourier transform
of the Berrry connection AA_a(k) = i
subroutine get_AA_R
!==================================================
!
!! 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_constants, only: dp, cmplx_0, cmplx_i
use w90_parameters, only: num_kpts, nntot, num_wann, wb, bk, timing_level, &
num_bands, ndimwin, nnlist, have_disentangled, &
transl_inv, nncell, effective_model
use w90_postw90_common, only: nrpts
use w90_io, only: stdout, io_file_unit, io_error, io_stopwatch, &
seedname
use w90_comms, only: on_root, comms_bcast
complex(kind=dp), allocatable :: AA_q(:, :, :, :)
complex(kind=dp), allocatable :: AA_q_diag(:, :)
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, ivdum(3), ivdum_old(3)
integer, allocatable :: num_states(:)
real(kind=dp) :: m_real, m_imag, rdum1_real, rdum1_imag, &
rdum2_real, rdum2_imag, rdum3_real, rdum3_imag
logical :: nn_found
character(len=60) :: header
if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_AA_R', 1)
if (.not. allocated(AA_R)) then
allocate (AA_R(num_wann, num_wann, nrpts, 3))
else
if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_AA_R', 2)
return
end if
! Real-space position matrix elements read from file
!
if (effective_model) then
if (.not. allocated(HH_R)) call io_error( &
'Error in get_AA_R: Must read file'//trim(seedname)//'_HH_R.dat first')
AA_R = cmplx_0
if (on_root) then
write (stdout, '(/a)') ' Reading position matrix elements from file ' &
//trim(seedname)//'_AA_R.dat'
file_unit = io_file_unit()
open (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 io_error('Error in get_AA_R: orbital indices out of bounds')
endif
if (n > 1) then
if (ivdum(1) /= ivdum_old(1) .or. ivdum(2) /= ivdum_old(2) .or. &
ivdum(3) /= ivdum_old(3)) ir = ir + 1
endif
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
enddo
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 io_error('Error in get_AA_R: inconsistent nrpts values')
endif
endif
call comms_bcast(AA_R(1, 1, 1, 1), num_wann*num_wann*nrpts*3)
if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_AA_R', 2)
return
endif
! 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(num_wann, num_wann, num_kpts, 3))
allocate (AA_q_diag(num_wann, 3))
allocate (S_o(num_bands, num_bands))
allocate (S(num_wann, num_wann))
allocate (num_states(num_kpts))
do ik = 1, num_kpts
if (have_disentangled) then
num_states(ik) = ndimwin(ik)
else
num_states(ik) = num_wann
endif
enddo
mmn_in = io_file_unit()
open (unit=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) &
call io_error(trim(seedname)//'.mmn has wrong number of bands')
if (nkp_tmp .ne. num_kpts) &
call io_error(trim(seedname)//'.mmn has wrong number of k-points')
if (nntot_tmp .ne. nntot) &
call io_error &
(trim(seedname)//'.mmn has wrong number of nearest neighbours')
AA_q = cmplx_0
ik_prev = 0
! Composite loop over k-points ik (outer loop) and neighbors ik2 (inner)
do ncount = 1, num_kpts*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)
enddo
enddo
!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, nntot
if ((ik2 .eq. nnlist(ik, inn)) .and. &
(nnl .eq. nncell(1, ik, inn)) .and. &
(nnm .eq. nncell(2, ik, inn)) .and. &
(nnn .eq. nncell(3, ik, inn))) then
if (.not. nn_found) then
nn_found = .true.
nn = inn
else
call io_error('Error reading '//trim(seedname)//'.mmn.&
& More than one matching nearest neighbour found')
endif
endif
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 io_error('Neighbour not found')
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( &
ik, num_states(ik), &
nnlist(ik, nn), num_states(nnlist(ik, nn)), &
S_o, S)
! Berry connection matrix
! Assuming all neighbors of a given point are read in sequence!
!
if (transl_inv .and. ik .ne. ik_prev) AA_q_diag(:, :) = cmplx_0
do idir = 1, 3
AA_q(:, :, ik, idir) = AA_q(:, :, ik, idir) &
+ cmplx_i*wb(nn)*bk(idir, nn, ik)*S(:, :)
if (transl_inv) then
!
! Rewrite band-diagonal elements a la Eq.(31) of MV97
!
do i = 1, num_wann
AA_q_diag(i, idir) = AA_q_diag(i, idir) &
- wb(nn)*bk(idir, nn, ik)*aimag(log(S(i, i)))
enddo
endif
end do
! Assuming all neighbors of a given point are read in sequence!
if (nn_count == nntot) then !looped over all neighbors
do idir = 1, 3
if (transl_inv) then
do n = 1, num_wann
AA_q(n, n, ik, idir) = AA_q_diag(n, idir)
enddo
endif
!
! 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)
!
AA_q(:, :, ik, idir) = &
0.5_dp*(AA_q(:, :, ik, idir) &
+ conjg(transpose(AA_q(:, :, ik, idir))))
enddo
end if
ik_prev = ik
enddo !ncount
close (mmn_in)
call fourier_q_to_R(AA_q(:, :, :, 1), AA_R(:, :, :, 1))
call fourier_q_to_R(AA_q(:, :, :, 2), AA_R(:, :, :, 2))
call fourier_q_to_R(AA_q(:, :, :, 3), AA_R(:, :, :, 3))
endif !on_root
call comms_bcast(AA_R(1, 1, 1, 1), num_wann*num_wann*nrpts*3)
if (timing_level > 1 .and. on_root) call io_stopwatch('get_oper: get_AA_R', 2)
return
101 call io_error &
('Error: Problem opening input file '//trim(seedname)//'.mmn')
102 call io_error &
('Error: Problem reading input file '//trim(seedname)//'.mmn')
103 call io_error('Error in get_AA_R: problem opening file '// &
trim(seedname)//'_AA_R.dat')
end subroutine get_AA_R