subroutine sitesym_dis_extract_symmetry(sitesym, lambda, umat, zmat, ik, n, num_bands, num_wann, &
stdout, error, comm)
!================================================!
!
! 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_wannier90_types, only: sitesym_type
implicit none
! arguments
type(sitesym_type), intent(in) :: sitesym
type(w90_error_type), allocatable, intent(out) :: error
type(w90_comm_type), intent(in) :: comm
integer, intent(in) :: num_bands
integer, intent(in) :: stdout
integer, intent(in) :: num_wann
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)
! local variables
complex(kind=dp), allocatable :: umatnew(:, :) !(num_bands, num_wann)
complex(kind=dp), allocatable :: ZU(:, :) !(num_bands, num_wann)
complex(kind=dp), allocatable :: deltaU(:, :) !(num_bands, num_wann)
complex(kind=dp), allocatable :: 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, ierr
integer, parameter :: niter = 50
allocate (umatnew(num_bands, num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating umatnew in sitesym_dis_extract_symmetry', comm)
return
end if
allocate (ZU(num_bands, num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating ZU in sitesym_dis_extract_symmetry', comm)
return
end if
allocate (deltaU(num_bands, num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating deltaU in sitesym_dis_extract_symmetry', comm)
return
end if
allocate (carr(num_bands), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating carr in sitesym_dis_extract_symmetry', comm)
return
end if
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
end if
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
end if
call set_error_fatal(error, 'error at sitesym_dis_extract_symmetry', comm)
return
end if
end if
! choose the larger eigenstate
umatnew(:, i) = V(1, 2)*umat(:, i) + V(2, 2)*deltaU(:, i)
end do ! i
call symmetrize_ukirr(num_wann, num_bands, sitesym%ik2ir(ik), num_bands, umatnew, sitesym, &
stdout, error, comm, n)
if (allocated(error)) return
umat(:, :) = umatnew(:, :)
end do ! iter
deallocate (umatnew, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating umatnew in sitesym_dis_extract_symmetry', comm)
return
end if
deallocate (ZU, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating ZU in sitesym_dis_extract_symmetry', comm)
return
end if
deallocate (deltaU, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating deltaU in sitesym_dis_extract_symmetry', comm)
return
end if
deallocate (carr, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating carr in sitesym_dis_extract_symmetry', comm)
return
end if
return
end subroutine sitesym_dis_extract_symmetry
