subroutine sitesym_symmetrize_zmatrix(sitesym, czmat, num_bands, num_kpts, lwindow_in)
!================================================!
!
! Z(k) <- \sum_{R} d^{+}(R,k) Z(Rk) d(R,k)
!
!================================================!
use w90_wannier90_types, only: sitesym_type
implicit none
! arguments
type(sitesym_type), intent(in) :: sitesym
integer, intent(in) :: num_bands
integer, intent(in) :: num_kpts
complex(kind=dp), intent(inout) :: czmat(num_bands, num_bands, num_kpts)
logical, intent(in) :: lwindow_in(num_bands, num_kpts)
! local variables
logical :: lfound(num_kpts)
integer :: ik, ir, isym, irk, nd
complex(kind=dp) :: cztmp(num_bands, num_bands)
complex(kind=dp) :: cmat1(num_bands, num_bands)
complex(kind=dp) :: cmat2(num_bands, num_bands)
lfound = .false.
do ir = 1, sitesym%nkptirr
ik = sitesym%ir2ik(ir)
nd = count(lwindow_in(:, ik))
lfound(ik) = .true.
do isym = 2, sitesym%nsymmetry
irk = sitesym%kptsym(isym, ir)
if (lfound(irk)) cycle
lfound(irk) = .true.
! cmat1 = Z(R,k)*d(R,k)
call zgemm('N', 'N', nd, nd, nd, cmplx_1, czmat(:, :, irk), num_bands, &
sitesym%d_matrix_band(:, :, isym, ir), num_bands, cmplx_0, cmat1, num_bands)
! cmat2 = d^{+}(R,k) Z(R,k) d(R,k) = d^{+}(R,k) cmat1
call zgemm('C', 'N', nd, nd, nd, cmplx_1, sitesym%d_matrix_band(:, :, isym, ir), &
num_bands, cmat1, num_bands, cmplx_0, cmat2, num_bands)
czmat(:, :, ik) = czmat(:, :, ik) + cmat2(:, :)
end do
cztmp(:, :) = czmat(:, :, ik)
do isym = 2, sitesym%nsymmetry
irk = sitesym%kptsym(isym, ir)
if (irk .ne. ik) cycle
call zgemm('N', 'N', nd, nd, nd, cmplx_1, cztmp, num_bands, &
sitesym%d_matrix_band(:, :, isym, ir), num_bands, cmplx_0, cmat1, num_bands)
! cmat2 = d^{+}(R,k) Z(R,k) d(R,k) = d^{+}(R,k) cmat1
call zgemm('C', 'N', nd, nd, nd, cmplx_1, sitesym%d_matrix_band(:, :, isym, ir), &
num_bands, cmat1, num_bands, cmplx_0, cmat2, num_bands)
czmat(:, :, ik) = czmat(:, :, ik) + cmat2(:, :)
end do
czmat(:, :, ik) = czmat(:, :, ik)/count(sitesym%kptsym(:, ir) .eq. ik)
end do
return
end subroutine sitesym_symmetrize_zmatrix