Uses guiding centres to pick phases which give a consistent choice of branch cut for the spread definition
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| complex(kind=dp), | intent(out) | :: | csheet(:,:,:) | Choice of phase |
||
| real(kind=dp), | intent(out) | :: | sheet(:,:,:) | Choice of branch cut |
||
| real(kind=dp), | intent(inout) | :: | rguide(:,:) | Guiding centres |
||
| integer, | intent(in) | :: | irguide | Zero if first call to this routine |
||
| real(kind=dp), | intent(in), | optional | :: | m_w(:,:,:) | Used in the Gamma point routines as an optimisation |
subroutine wann_phases(csheet, sheet, rguide, irguide, m_w)
!==================================================================!
!! Uses guiding centres to pick phases which give a
!! consistent choice of branch cut for the spread definition
! !
!===================================================================
use w90_constants, only: eps6
use w90_parameters, only: num_wann, nntot, neigh, &
nnh, bk, bka, num_kpts, timing_level, m_matrix, gamma_only
use w90_io, only: io_stopwatch
use w90_utility, only: utility_inv3
implicit none
complex(kind=dp), intent(out) :: csheet(:, :, :)
!! Choice of phase
real(kind=dp), intent(out) :: sheet(:, :, :)
!! Choice of branch cut
real(kind=dp), intent(inout) :: rguide(:, :)
!! Guiding centres
integer, intent(in) :: irguide
!! Zero if first call to this routine
real(kind=dp), intent(in), optional :: m_w(:, :, :)
!! Used in the Gamma point routines as an optimisation
!local
complex(kind=dp) :: csum(nnh)
real(kind=dp) :: xx(nnh)
real(kind=dp) :: smat(3, 3), svec(3), sinv(3, 3)
real(kind=dp) :: xx0, det, brn
complex(kind=dp) :: csumt
integer :: loop_wann, na, nkp, i, j, nn, ind, m, nkp_loc
if (timing_level > 1 .and. on_root) call io_stopwatch('wann: phases', 1)
csum = cmplx_0; xx = 0.0_dp
! report problem to solve
! for each band, csum is determined and then its appropriate
! guiding center rguide(3,nwann)
do loop_wann = 1, num_wann
if (.not. present(m_w)) then
! get average phase for each unique bk direction
if (gamma_only) then
do na = 1, nnh
csum(na) = cmplx_0
do nkp_loc = 1, counts(my_node_id)
nkp = nkp_loc + displs(my_node_id)
nn = neigh(nkp, na)
csum(na) = csum(na) + m_matrix(loop_wann, loop_wann, nn, nkp_loc)
enddo
enddo
else
do na = 1, nnh
csum(na) = cmplx_0
do nkp_loc = 1, counts(my_node_id)
nkp = nkp_loc + displs(my_node_id)
nn = neigh(nkp, na)
csum(na) = csum(na) + m_matrix_loc(loop_wann, loop_wann, nn, nkp_loc)
enddo
enddo
endif
else
do na = 1, nnh
csum(na) = cmplx_0
do nkp_loc = 1, counts(my_node_id)
nkp = nkp_loc + displs(my_node_id)
nn = neigh(nkp, na)
csum(na) = csum(na) &
+ cmplx(m_w(loop_wann, loop_wann, 2*nn - 1), m_w(loop_wann, loop_wann, 2*nn), dp)
enddo
enddo
end if
call comms_allreduce(csum(1), nnh, 'SUM')
! now analyze that information to get good guess at
! wannier center
! write(*,*)
! do na=1,nnh
! write(*,'a,3f10.5,a,2f10.5)')
! & ' bka=',(bka(j,na),j=1,3),' csum=',csum(na)
! end do
! problem is to find a real-space 3-vector rguide such that
! phase of csum(nn) ~= phase of exp[ -i bka(nn) dot rguide ]
! or, letting
! xx(nn) = - Im ln csum(nn) (modulo 2*pi)
! then
! bka(nn) dot rguide ~= xx(nn)
!
! we take an arbitrary branch cut for first three xx(nn)
! and determine rguide from these; then for each additional bka
! vector, we first determine the most consistent branch cut,
! and then update rguide
!
! rguide is obtained by minimizing
! sum_nn [ bka(nn) dot rguide - xx(nn) ] ^2
! or, setting the derivative with respect to rcenter to zero,
! sum_i smat(j,i) * rguide(i,nwann) = svec(j)
! where
! smat(j,i) = sum_nn bka(j,nn) * bka(i,nn)
! svec(j) = sum_nn bka(j,nn) * xx(nn)
! initialize smat and svec
smat = 0.0_dp
svec = 0.0_dp
do nn = 1, nnh
if (nn .le. 3) then
! obtain xx with arbitrary branch cut choice
xx(nn) = -aimag(log(csum(nn)))
else
! obtain xx with branch cut choice guided by rguide
xx0 = 0.0_dp
do j = 1, 3
xx0 = xx0 + bka(j, nn)*rguide(j, loop_wann)
enddo
! xx0 is expected value for xx
! csumt = exp (ci * xx0)
csumt = exp(cmplx_i*xx0)
! csumt has opposite of expected phase of csum(nn)
xx(nn) = xx0 - aimag(log(csum(nn)*csumt))
endif
! write(*,'(a,i5,3f7.3,2f10.5)') 'nn, bka, xx, mag =',
! 1 nn,(bka(j,nn),j=1,3),xx(nn),abs(csum(nn))/float(num_kpts)
! update smat and svec
do j = 1, 3
do i = 1, 3
smat(j, i) = smat(j, i) + bka(j, nn)*bka(i, nn)
enddo
svec(j) = svec(j) + bka(j, nn)*xx(nn)
enddo
if (nn .ge. 3) then
! determine rguide
call utility_inv3(smat, sinv, det)
! the inverse of smat is sinv/det
if (abs(det) .gt. eps6) then
! to check that the first nn bka vectors are not
! linearly dependent - this is a change from original code
if (irguide .ne. 0) then
do j = 1, 3
rguide(j, loop_wann) = 0.0_dp
do i = 1, 3
rguide(j, loop_wann) = rguide(j, loop_wann) + sinv(j, i) &
*svec(i)/det
enddo
enddo
endif
endif
endif
enddo
enddo
! obtain branch cut choice guided by rguid
sheet = 0.0_dp
do nkp = 1, num_kpts
do nn = 1, nntot
do loop_wann = 1, num_wann
! sheet (loop_wann, nn, nkp) = 0.d0
do j = 1, 3
sheet(loop_wann, nn, nkp) = sheet(loop_wann, nn, nkp) &
+ bk(j, nn, nkp)*rguide(j, loop_wann)
enddo
! csheet (loop_wann, nn, nkp) = exp (ci * sheet (loop_wann, nn, nkp) )
enddo
enddo
enddo
csheet = exp(cmplx_i*sheet)
! now check that we picked the proper sheet for the log
! of m_matrix. criterion: q_n^{k,b}=Im(ln(M_nn^{k,b})) + b \cdot r_n are
! circa 0 for a good solution, circa multiples of 2 pi for a bad one.
! I use the guiding center, instead of r_n, to understand which could be
! right sheet
rnkb = 0.0_dp
do nkp = 1, num_kpts
do nn = 1, nntot
do m = 1, num_wann
! rnkb (m, nn, nkp) = 0.0_dp
brn = 0.0_dp
do ind = 1, 3
brn = brn + bk(ind, nn, nkp)*rguide(ind, m)
enddo
rnkb(m, nn, nkp) = rnkb(m, nn, nkp) + brn
enddo
enddo
enddo
! write ( stdout , * ) ' '
! write ( stdout , * ) ' PHASES ARE SET USING THE GUIDING CENTERS'
! write ( stdout , * ) ' '
! do nkp = 1, num_kpts
! do n = 1, num_wann
! do nn = 1, nntot
! pherr = aimag(log(csheet(n,nn,nkp)*m_matrix(n,n,nn,nkp))) &
! - sheet(n,nn,nkp)+rnkb(n,nn,nkp)-aimag(log(m_matrix(n,n,nn,nkp)))
! enddo
! enddo
! enddo
if (timing_level > 1 .and. on_root) call io_stopwatch('wann: phases', 2)
return
end subroutine wann_phases