subroutine sitesym_dis_extract_symmetry(ik, n, zmat, lambda, umat)
!==================================================================!
! !
! minimize Omega_I by steepest descendent !
! !
! delta U_{mu I}(k) = Z_{mu mu'}*U_{mu' I}(k) !
! - \sum_{J} lambda_{JI} U_{mu J}(k) !
! lambda_{JI}=U^{*}_{mu J} Z_{mu mu'} U_{mu' I} !
! !
!==================================================================!
use w90_parameters, only: num_bands, num_wann
implicit none
integer, intent(in) :: ik, n
complex(kind=dp), intent(in) :: zmat(num_bands, num_bands)
complex(kind=dp), intent(out) :: lambda(num_wann, num_wann)
complex(kind=dp), intent(inout) :: umat(num_bands, num_wann)
complex(kind=dp) :: umatnew(num_bands, num_wann)
complex(kind=dp) :: ZU(num_bands, num_wann)
complex(kind=dp) :: deltaU(num_bands, num_wann), carr(num_bands)
integer :: i, m, INFO, IFAIL(2), IWORK(5*2)
complex(kind=dp) :: HP(3), SP(3), V(2, 2), CWORK(2*2)
real(kind=dp) :: W(2), RWORK(7*2), sp3
integer :: iter
integer, parameter :: niter = 50
do iter = 1, niter
! Z*U
call zgemm('N', 'N', n, num_wann, n, cmplx_1, &
zmat, num_bands, umat, num_bands, &
cmplx_0, ZU, num_bands)
! lambda = U^{+}*Z*U
call zgemm('C', 'N', num_wann, num_wann, n, cmplx_1, &
umat, num_bands, ZU, num_bands, &
cmplx_0, lambda, num_wann)
deltaU(:, :) = ZU(:, :) - matmul(umat, lambda)
if (sum(abs(deltaU(:n, :))) .lt. 1e-10) return
! band-by-band minimization
do i = 1, num_wann
! diagonalize 2x2 matrix
HP(1) = real(dot_product(umat(1:n, i), ZU(1:n, i)), kind=dp)
HP(2) = dot_product(ZU(1:n, i), deltaU(1:n, i)) ! (1,2) matrix element
carr(1:n) = matmul(zmat(1:n, 1:n), deltaU(1:n, i))
HP(3) = real(dot_product(deltaU(1:n, i), carr(1:n)), kind=dp) ! (2,2)
SP(1) = real(dot_product(umat(1:n, i), umat(1:n, i)), kind=dp)
SP(2) = dot_product(umat(1:n, i), deltaU(1:n, i))
SP(3) = real(dot_product(deltaU(1:n, i), deltaU(1:n, i)), kind=dp)
sp3 = real(SP(3), kind=dp)
if (abs(sp3) .lt. 1e-10) then
umatnew(:, i) = umat(:, i)
cycle
endif
call ZHPGVX(1, 'V', 'A', 'U', 2, HP, SP, 0.0_dp, 0.0_dp, 0, 0, &
-1.0_dp, m, W, V, 2, CWORK, RWORK, IWORK, IFAIL, INFO)
if (INFO .ne. 0) then
write (stdout, *) 'error in sitesym_dis_extract_symmetry: INFO=', INFO
if (INFO .gt. 0) then
if (INFO .le. 2) then
write (stdout, *) INFO, ' eigenvectors failed to converge'
write (stdout, *) IFAIL(1:INFO)
else
write (stdout, *) ' S is not positive definite'
write (stdout, *) 'sp3=', sp3
endif
call io_error('error at sitesym_dis_extract_symmetry')
endif
endif
! choose the larger eigenstate
umatnew(:, i) = V(1, 2)*umat(:, i) + V(2, 2)*deltaU(:, i)
enddo ! i
call symmetrize_ukirr(ik2ir(ik), num_bands, umatnew, n)
umat(:, :) = umatnew(:, :)
enddo ! iter
return
end subroutine sitesym_dis_extract_symmetry