subroutine sitesym_symmetrize_gradient(imode, grad)
!==================================================================!
use w90_parameters, only: num_wann, num_kpts
use w90_utility, only: utility_zgemm
implicit none
integer, intent(in) :: imode
complex(kind=dp), intent(inout) :: grad(num_wann, num_wann, num_kpts)
integer :: ik, ir, isym, irk, ngk
complex(kind=dp) :: grad_total(num_wann, num_wann)
complex(kind=dp) :: cmat1(num_wann, num_wann)
complex(kind=dp) :: cmat2(num_wann, num_wann)
logical :: lfound(num_kpts)
if (imode .eq. 1) then
lfound = .false.
do ir = 1, nkptirr
ik = ir2ik(ir)
grad_total = grad(:, :, ik)
lfound(ik) = .true.
do isym = 2, nsymmetry
irk = kptsym(isym, ir)
if (lfound(irk)) cycle
lfound(irk) = .true.
!
! cmat1 = D(R,k)^{+} G(Rk) D(R,k)
! cmat2 = D(R,k)^{\dagger} G(Rk)
!
call utility_zgemm(cmat2, d_matrix_wann(:, :, isym, ir), 'C', &
grad(:, :, irk), 'N', num_wann)
call utility_zgemm(cmat1, cmat2, 'N', &
d_matrix_wann(:, :, isym, ir), 'N', num_wann)
grad_total = grad_total + cmat1
enddo
grad(:, :, ik) = grad_total
enddo
do ik = 1, num_kpts
if (ir2ik(ik2ir(ik)) .ne. ik) grad(:, :, ik) = 0
enddo
endif ! if (imode.eq.1)
!
! grad -> 1/N_{R'} \sum_{R'} D^{+}(R',k) grad D(R',k)
! where R' k = k
!
do ir = 1, nkptirr
ik = ir2ik(ir)
ngk = count(kptsym(:, ir) .eq. ik)
if (ngk .eq. 1) cycle
grad_total = grad(:, :, ik)
do isym = 2, nsymmetry
if (kptsym(isym, ir) .ne. ik) cycle
!
! calculate cmat1 = D^{+}(R,k) G(Rk) D(R,k)
!
! step 1: cmat2 = G(Rk) D(R,k)
call utility_zgemm(cmat2, grad(:, :, ik), 'N', &
d_matrix_wann(:, :, isym, ir), 'N', num_wann)
! step 2: cmat1 = D^{+}(R,k) * cmat2
call utility_zgemm(cmat1, d_matrix_wann(:, :, isym, ir), 'C', &
cmat2, 'N', num_wann)
grad_total = grad_total + cmat1
enddo
grad(:, :, ik) = grad_total/ngk
enddo
return
end subroutine sitesym_symmetrize_gradient