wann_phases Subroutine

  subroutine wann_phases(csheet, sheet, rguide, irguide, num_wann, kmesh_info, num_kpts, &
                         use_ss_functional, m_matrix_loc, timing_level, iprint, timer, &
                         nkrank, global_k, error, comm, 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, cmplx_0, cmplx_i
    use w90_io, only: io_stopwatch_start, io_stopwatch_stop
    use w90_utility, only: utility_inv3
    use w90_comms, only: comms_allreduce, w90_comm_type
    use w90_types, only: kmesh_info_type, timer_list_type

    implicit none

    ! arguments
    type(w90_comm_type), intent(in) :: comm
    type(kmesh_info_type), intent(in) :: kmesh_info
    type(timer_list_type), intent(inout) :: timer
    type(w90_error_type), allocatable, intent(out) :: error

    integer, intent(in) :: timing_level
    integer, intent(in) :: num_wann
    integer, intent(in) :: num_kpts
    integer, intent(in) :: irguide !! Zero if first call to this routine
    integer, intent(in) :: iprint
    integer, intent(in) :: nkrank
    integer, intent(in) :: global_k(:)

    real(kind=dp), intent(out) :: sheet(:, :, :) !! Choice of branch cut
    real(kind=dp), intent(inout) :: rguide(:, :) !! Guiding centres
    real(kind=dp), intent(in), optional :: m_w(:, :, :)

    complex(kind=dp), intent(out) :: csheet(:, :, :) !! Choice of phase
    complex(kind=dp), intent(in) :: m_matrix_loc(:, :, :, :)
    logical, intent(in) :: use_ss_functional

    !local
    complex(kind=dp) :: csum(kmesh_info%nnh)
    real(kind=dp) :: xx(kmesh_info%nnh)
    real(kind=dp) :: smat(3, 3), svec(3), sinv(3, 3)
    real(kind=dp) :: xx0, det
    complex(kind=dp) :: csumt
    integer :: loop_wann, na, nkp, i, j, nn, nkp_loc

    if (timing_level > 1 .and. iprint > 0) call io_stopwatch_start('wann: phases', timer)

    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
        do na = 1, kmesh_info%nnh
          csum(na) = cmplx_0
          do nkp_loc = 1, nkrank
            nkp = global_k(nkp_loc)
            nn = kmesh_info%neigh(nkp, na)
            csum(na) = csum(na) + m_matrix_loc(loop_wann, loop_wann, nn, nkp_loc)
          end do
        end do
      else
        do na = 1, kmesh_info%nnh
          csum(na) = cmplx_0
          do nkp_loc = 1, nkrank
            nkp = global_k(nkp_loc)
            nn = kmesh_info%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)
          end do
        end do
      end if ! m_w present

      call comms_allreduce(csum(1), kmesh_info%nnh, 'SUM', error, comm)
      if (allocated(error)) return

      ! 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, kmesh_info%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 + kmesh_info%bka(j, nn)*rguide(j, loop_wann)
          end do
          ! 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))
        end if

        ! update smat and svec
        do j = 1, 3
          do i = 1, 3
            smat(j, i) = smat(j, i) + kmesh_info%bka(j, nn)*kmesh_info%bka(i, nn)
          end do
          svec(j) = svec(j) + kmesh_info%bka(j, nn)*xx(nn)
        end do

        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
                end do
              end do
            end if
          end if
        end if

      end do !nnh
    end do !loop_wann

    ! obtain branch cut choice guided by rguide
    sheet = 0.0_dp
    if (use_ss_functional) then
      do nn = 1, kmesh_info%nntot
        do loop_wann = 1, num_wann
          do j = 1, 3
            sheet(loop_wann, nn, 1) = sheet(loop_wann, nn, 1) &
                                      + kmesh_info%bk(j, nn, 1)*rguide(j, loop_wann)
          end do
        end do
      end do
    else
      do nkp = 1, num_kpts
        do nn = 1, kmesh_info%nntot
          do loop_wann = 1, num_wann
            do j = 1, 3
              sheet(loop_wann, nn, nkp) = sheet(loop_wann, nn, nkp) &
                                          + kmesh_info%bk(j, nn, nkp)*rguide(j, loop_wann)
            end do
          end do
        end do
      end do
    end if
    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, kmesh_info%nntot
    !     do m = 1, num_wann
    !       !           rnkb (m, nn, nkp) = 0.0_dp
    !       brn = 0.0_dp
    !       do ind = 1, 3
    !         brn = brn + kmesh_info%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. iprint > 0) call io_stopwatch_stop('wann: phases', timer)
    return
  end subroutine wann_phases