subroutine sitesym_symmetrize_gradient(sitesym, grad, imode, num_kpts, num_wann, error, comm)
!================================================!
use w90_error, only: w90_error_type, set_error_fatal
use w90_utility, only: utility_zgemm
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) :: imode, num_wann, num_kpts
complex(kind=dp), intent(inout) :: grad(:, :, :) !(num_wann, num_wann, num_kpts)
! local variables
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, sitesym%nkptirr
ik = sitesym%ir2ik(ir)
grad_total = grad(:, :, ik)
lfound(ik) = .true.
do isym = 2, sitesym%nsymmetry
irk = sitesym%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, sitesym%d_matrix_wann(:, :, isym, ir), 'C', grad(:, :, irk), 'N', num_wann)
call utility_zgemm(cmat1, cmat2, 'N', sitesym%d_matrix_wann(:, :, isym, ir), 'N', num_wann)
grad_total = grad_total + cmat1
end do
grad(:, :, ik) = grad_total
end do
do ik = 1, num_kpts
if (sitesym%ir2ik(sitesym%ik2ir(ik)) .ne. ik) grad(:, :, ik) = 0
end do
elseif (imode .eq. 2) then
! JJ, 20 July 2022, note:
! previously the following algorithm was *also applied* after the above for "mode 1"
! changed such that two algorithms are mutually exclusive.
! old results (test case testw90_disentanglement_sawfs) require mode 1 followed by mode 2
! see call in wannierise's wann_domega() routine
! surely the two modes do the same thing??? if not, then replace elseif with endif as before --JJ
!
! grad -> 1/N_{R'} \sum_{R'} D^{+}(R',k) grad D(R',k)
! where R' k = k
!
do ir = 1, sitesym%nkptirr
ik = sitesym%ir2ik(ir)
ngk = count(sitesym%kptsym(:, ir) .eq. ik)
if (ngk .eq. 1) cycle
grad_total = grad(:, :, ik)
do isym = 2, sitesym%nsymmetry
if (sitesym%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', sitesym%d_matrix_wann(:, :, isym, ir), 'N', num_wann)
! step 2: cmat1 = D^{+}(R,k) * cmat2
call utility_zgemm(cmat1, sitesym%d_matrix_wann(:, :, isym, ir), 'C', cmat2, 'N', num_wann)
grad_total = grad_total + cmat1
end do
grad(:, :, ik) = grad_total/ngk
end do
else
call set_error_fatal(error, 'unknown mode argument in sitesym_symmetrize_gradient', comm)
return
end if
return
end subroutine sitesym_symmetrize_gradient