subroutine symmetrize_ukirr(ir, ndim, umat, n)
!==================================================================!
! !
! calculate u~(k)=1/N_{R'} \sum_{R'} d^{+}(R',k) u(k) D(R',k) !
! where R'k = k !
! and orthonormalize it !
! !
!==================================================================!
use w90_parameters, only: num_wann, num_bands, symmetrize_eps
implicit none
integer, intent(in) :: ir, ndim
complex(kind=dp), intent(inout) :: umat(ndim, num_wann)
integer, optional, intent(in) :: n
integer :: isym, ngk, i, iter, ntmp
integer, parameter :: niter = 100
real(kind=dp) :: diff
complex(kind=dp) :: usum(ndim, num_wann)
complex(kind=dp) :: cmat_sub(ndim, num_wann)
complex(kind=dp) :: cmat(ndim, num_wann)
complex(kind=dp) :: cmat2(num_wann, num_wann)
!write(stdout,"(a)") '-- symmetrize_ukirr --'
if (present(n)) then
if (ndim .ne. num_bands) call io_error('ndim!=num_bands')
ntmp = n
else
if (ndim .ne. num_wann) call io_error('ndim!=num_wann')
ntmp = ndim
endif
ngk = count(kptsym(:, ir) .eq. ir2ik(ir))
if (ngk .eq. 1) then
call orthogonalize_u(ndim, num_wann, umat, ntmp)
return
endif
do iter = 1, niter
usum(:, :) = 0
cmat2(:, :) = 0
do i = 1, num_wann
cmat2(i, i) = cmat2(i, i) + ngk
enddo
do isym = 1, nsymmetry
if (kptsym(isym, ir) .ne. ir2ik(ir)) cycle
!
! cmat = d^{+}(R,k) U(k) D(R,k)
! size of umat: umat(ndim,num_wann)
!
! cmat_sub = U(k) D(R,k)
call zgemm('N', 'N', ntmp, num_wann, num_wann, cmplx_1, &
umat, ndim, d_matrix_wann(:, :, isym, ir), num_wann, &
cmplx_0, cmat_sub, ndim)
! cmat = d^{+}(R,k) * cmat_sub
call zgemm('C', 'N', ntmp, num_wann, ntmp, cmplx_1, &
d_matrix_band(:, :, isym, ir), ndim, cmat_sub, ndim, &
cmplx_0, cmat, ndim)
usum(:, :) = usum(:, :) + cmat(:, :)
! check
cmat2(:, :) = cmat2(:, :) - &
matmul(conjg(transpose(umat(:ntmp, :))), cmat(:ntmp, :))
enddo ! isym
diff = sum(abs(cmat2))
if (diff .lt. symmetrize_eps) exit
if (iter .eq. niter) then
write (stdout, "(a)") 'Error in symmetrize_u: not converged'
write (stdout, "(a)") 'Either eps is too small or specified irreps is not'
write (stdout, "(a)") ' compatible with the bands'
write (stdout, "(a,2e20.10)") 'diff,eps=', diff, symmetrize_eps
call io_error('symmetrize_ukirr: not converged')
endif
usum = usum/ngk
call orthogonalize_u(ndim, num_wann, usum, ntmp)
umat(:, :) = usum
enddo ! iter
return
end subroutine symmetrize_ukirr