subroutine wann_domega(csheet, sheet, rave, num_wann, kmesh_info, num_kpts, wann_slwf, use_ss_functional, &
lsitesymmetry, ln_tmp_loc, m_matrix_loc, rnkb_loc, cdodq_loc, &
lambda_loc, timing_level, sitesym, timer, nkrank, global_k, error, comm, &
iprint, cdodq)
!================================================!
!
! Calculate the Gradient of the Wannier Function spread
!
! Modified by Valerio Vitale for the SLWF+C method (PRB 90, 165125)
! Jun 2018, based on previous work by Charles T. Johnson and
! Radu Miron at Imperial College London
!================================================
use w90_comms, only: comms_gatherv, comms_allreduce, w90_comm_type, mpirank
use w90_constants, only: cmplx_0
use w90_io, only: io_stopwatch_start, io_stopwatch_stop
use w90_sitesym, only: sitesym_symmetrize_gradient
use w90_types, only: kmesh_info_type, timer_list_type
use w90_wannier90_types, only: wann_slwf_type, sitesym_type
implicit none
! arguments
type(kmesh_info_type), intent(in) :: kmesh_info
type(sitesym_type), intent(in) :: sitesym
type(timer_list_type), intent(inout) :: timer
type(w90_comm_type), intent(in) :: comm
type(w90_error_type), allocatable, intent(out) :: error
type(wann_slwf_type), intent(inout) :: wann_slwf
integer, intent(in) :: num_wann
integer, intent(in) :: num_kpts
integer, intent(in) :: timing_level, iprint
integer, intent(in) :: nkrank
integer, intent(in) :: global_k(:)
logical, intent(in) :: use_ss_functional
real(kind=dp), intent(in) :: sheet(:, :, :)
real(kind=dp), intent(out) :: rave(:, :)
real(kind=dp), intent(inout) :: ln_tmp_loc(:, :, :)
real(kind=dp), intent(inout) :: rnkb_loc(:, :, :)
real(kind=dp), intent(in) :: lambda_loc
! as we work on the local cdodq, returning the full cdodq array is now made optional
complex(kind=dp), intent(out), optional :: cdodq(:, :, :)
complex(kind=dp), intent(in) :: csheet(:, :, :)
complex(kind=dp), intent(in) :: m_matrix_loc(:, :, :, :)
complex(kind=dp), intent(out) :: cdodq_loc(:, :, :)
logical, intent(in) :: lsitesymmetry
! local
complex(kind=dp), allocatable :: cr(:, :)
complex(kind=dp), allocatable :: crt(:, :)
real(kind=dp), allocatable :: r0kb(:, :, :)
complex(kind=dp), allocatable :: sum_mnn(:, :)
integer :: iw, ind, nkp, nn, m, n, ierr, nkp_loc, cnn, cnn2
complex(kind=dp) :: mnn
integer :: my_node_id
my_node_id = mpirank(comm)
if (timing_level > 1 .and. iprint > 0) call io_stopwatch_start('wann: domega', timer)
allocate (cr(num_wann, num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating cr in wann_domega', comm)
return
end if
allocate (crt(num_wann, num_wann), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating crt in wann_domega', comm)
return
end if
if (wann_slwf%selective_loc .and. wann_slwf%constrain) then
allocate (r0kb(num_wann, kmesh_info%nntot, num_kpts), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating r0kb in wann_domega', comm)
return
end if
end if
if (use_ss_functional) then
allocate (sum_mnn(num_wann, kmesh_info%nntot), stat=ierr)
if (ierr /= 0) then
call set_error_alloc(error, 'Error in allocating sum_mnn in wann_domega', comm)
return
end if
sum_mnn = 0.0_dp
do nkp_loc = 1, nkrank
nkp = global_k(nkp_loc)
do nn = 1, kmesh_info%nntot
cnn = kmesh_info%nnord(nn, nkp) ! enforce uniform order of bk vectors
do n = 1, num_wann
sum_mnn(n, nn) = sum_mnn(n, nn) + csheet(n, nn, 1)*m_matrix_loc(n, n, cnn, nkp_loc)
end do
end do
end do
call comms_allreduce(sum_mnn(1, 1), num_wann*kmesh_info%nntot, 'SUM', error, comm)
if (allocated(error)) return
sum_mnn = sum_mnn/real(num_kpts, dp)
! k-index is always 1 in SS method (k summation alread accomplished)
do nn = 1, kmesh_info%nntot
ln_tmp_loc(:, nn, 1) = kmesh_info%wb(nn)*(aimag(log(sum_mnn(:, nn))) - sheet(:, nn, 1))
end do
rave = 0.0_dp
do iw = 1, num_wann
do ind = 1, 3
do nn = 1, kmesh_info%nntot
rave(ind, iw) = rave(ind, iw) + kmesh_info%bk(ind, nn, 1)*ln_tmp_loc(iw, nn, 1)
end do
end do
end do
rave = -rave
rnkb_loc = 0.0_dp
do nkp_loc = 1, nkrank
do nn = 1, kmesh_info%nntot
do n = 1, num_wann
rnkb_loc(n, nn, nkp_loc) = sum(kmesh_info%bk(:, nn, 1)*rave(:, n))
end do
end do
end do
cdodq_loc = cmplx_0
do nkp_loc = 1, nkrank
do n = 1, num_wann
do m = 1, num_wann
do nn = 1, kmesh_info%nntot
nkp = global_k(nkp_loc)
cnn = kmesh_info%nnord(nn, nkp) ! enforce uniform order of bk vectors
cnn2 = kmesh_info%nnrev(nn, nkp) ! (b-vector nn2) is opposite to/negative of (b-vector nn)
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) + &
kmesh_info%wb(nn)*m_matrix_loc(m, n, cnn, nkp_loc)* &
conjg(sum_mnn(n, nn)/csheet(n, nn, 1))
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) - &
kmesh_info%wb(nn)*conjg(m_matrix_loc(n, m, cnn2, nkp_loc))* &
conjg(sum_mnn(m, nn)/csheet(m, nn, 1))
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) - &
kmesh_info%wb(nn)*conjg(m_matrix_loc(n, m, cnn, nkp_loc))* &
sum_mnn(m, nn)/csheet(m, nn, 1)
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) + &
kmesh_info%wb(nn)*m_matrix_loc(m, n, cnn2, nkp_loc)* &
sum_mnn(n, nn)/csheet(n, nn, 1)
end do
end do
end do
end do
cdodq_loc = cdodq_loc/real(num_kpts, dp)
if (present(cdodq)) then
! each process communicates its result to other processes
cdodq(:, :, :) = 0.0_dp
do nkp_loc = 1, nkrank
nkp = global_k(nkp_loc)
cdodq(:, :, nkp) = cdodq_loc(:, :, nkp_loc)
end do
call comms_allreduce(cdodq(1, 1, 1), num_wann*num_wann*num_kpts, 'SUM', error, comm)
if (allocated(error)) return
end if
deallocate (sum_mnn, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating sum_mnn in wann_domega', comm)
return
end if
else ! not Stengel-Spaldin
do nkp_loc = 1, nkrank
nkp = global_k(nkp_loc)
do nn = 1, kmesh_info%nntot
do n = 1, num_wann
! Note that this ln_tmp is defined differently wrt the one in wann_omega
ln_tmp_loc(n, nn, nkp_loc) = kmesh_info%wb(nn)*(aimag(log(csheet(n, nn, nkp) &
*m_matrix_loc(n, n, nn, nkp_loc))) - sheet(n, nn, nkp))
end do
end do
end do
! recalculate rave
rave = 0.0_dp
do iw = 1, num_wann
do ind = 1, 3
do nkp_loc = 1, nkrank
nkp = global_k(nkp_loc)
do nn = 1, kmesh_info%nntot
rave(ind, iw) = rave(ind, iw) + kmesh_info%bk(ind, nn, nkp) &
*ln_tmp_loc(iw, nn, nkp_loc)
end do
end do
end do
end do
rave = -rave/real(num_kpts, dp)
call comms_allreduce(rave(1, 1), num_wann*3, 'SUM', error, comm)
if (allocated(error)) return
! b.r_0n are calculated
if (wann_slwf%selective_loc .and. wann_slwf%constrain) then
r0kb = 0.0_dp
do nkp_loc = 1, nkrank
nkp = global_k(nkp_loc)
do nn = 1, kmesh_info%nntot
do n = 1, num_wann
r0kb(n, nn, nkp_loc) = sum(kmesh_info%bk(:, nn, nkp) &
*wann_slwf%centres(n, :))
end do
end do
end do
end if
rnkb_loc = 0.0_dp
do nkp_loc = 1, nkrank
nkp = global_k(nkp_loc)
do nn = 1, kmesh_info%nntot
do n = 1, num_wann
rnkb_loc(n, nn, nkp_loc) = sum(kmesh_info%bk(:, nn, nkp)*rave(:, n))
end do
end do
end do
! cd0dq(m,n,nkp) is calculated
cdodq_loc = cmplx_0
cr = cmplx_0
crt = cmplx_0
do nkp_loc = 1, nkrank
nkp = global_k(nkp_loc)
do nn = 1, kmesh_info%nntot
do n = 1, num_wann ! R^{k,b} and R~^{k,b} have columns of zeroes for the non-objective Wannier functions
mnn = m_matrix_loc(n, n, nn, nkp_loc)
crt(:, n) = m_matrix_loc(1:num_wann, n, nn, nkp_loc)/mnn !JJ potential for division by zero
cr(:, n) = m_matrix_loc(1:num_wann, n, nn, nkp_loc)*conjg(mnn)
end do
if (wann_slwf%selective_loc) then
do n = 1, num_wann
do m = 1, num_wann
if (m <= wann_slwf%slwf_num) then
if (n <= wann_slwf%slwf_num) then
! A[R^{k,b}]=(R-Rdag)/2
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) &
+ kmesh_info%wb(nn)*0.5_dp*(cr(m, n) - conjg(cr(n, m)))
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) &
- (crt(m, n)*ln_tmp_loc(n, nn, nkp_loc) &
+ conjg(crt(n, m)*ln_tmp_loc(m, nn, nkp_loc))) &
*cmplx(0.0_dp, -0.5_dp, kind=dp)
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) &
- (crt(m, n)*rnkb_loc(n, nn, nkp_loc) &
+ conjg(crt(n, m)*rnkb_loc(m, nn, nkp_loc))) &
*cmplx(0.0_dp, -0.5_dp, kind=dp)
if (wann_slwf%constrain) then
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) + lambda_loc &
*(crt(m, n)*ln_tmp_loc(n, nn, nkp_loc) &
+ conjg(crt(n, m)*ln_tmp_loc(m, nn, nkp_loc))) &
*cmplx(0.0_dp, -0.5_dp, kind=dp)
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) &
+ kmesh_info%wb(nn)*lambda_loc &
*(crt(m, n)*rnkb_loc(n, nn, nkp_loc) &
+ conjg(crt(n, m)*rnkb_loc(m, nn, nkp_loc))) &
*cmplx(0.0_dp, -0.5_dp, kind=dp)
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) - lambda_loc &
*(crt(m, n)*ln_tmp_loc(n, nn, nkp_loc) &
+ conjg(crt(n, m))*ln_tmp_loc(m, nn, nkp_loc)) &
*cmplx(0.0_dp, -0.5_dp, kind=dp)
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) &
- kmesh_info%wb(nn)*lambda_loc &
*(r0kb(n, nn, nkp_loc)*crt(m, n) &
+ r0kb(m, nn, nkp_loc)*conjg(crt(n, m))) &
*cmplx(0.0_dp, -0.5_dp, kind=dp)
end if
else
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) - kmesh_info%wb(nn) &
*0.5_dp*conjg(cr(n, m)) &
- conjg(crt(n, m)*(ln_tmp_loc(m, nn, nkp_loc) &
+ kmesh_info%wb(nn)*rnkb_loc(m, nn, nkp_loc))) &
*cmplx(0.0_dp, -0.5_dp, kind=dp)
if (wann_slwf%constrain) then
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) + lambda_loc &
*conjg(crt(n, m)*(ln_tmp_loc(m, nn, nkp_loc) &
+ kmesh_info%wb(nn)*rnkb_loc(m, nn, nkp_loc))) &
*cmplx(0.0_dp, -0.5_dp, kind=dp) &
- lambda_loc*(conjg(crt(n, m)) &
*ln_tmp_loc(m, nn, nkp_loc)) &
*cmplx(0.0_dp, -0.5_dp, kind=dp)
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) &
- kmesh_info%wb(nn)*lambda_loc &
*r0kb(m, nn, nkp_loc)*conjg(crt(n, m)) &
*cmplx(0.0_dp, -0.5_dp, kind=dp)
end if
end if
else if (n <= wann_slwf%slwf_num) then
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) &
+ kmesh_info%wb(nn)*cr(m, n)*0.5_dp &
- crt(m, n)*(ln_tmp_loc(n, nn, nkp_loc) &
+ kmesh_info%wb(nn)*rnkb_loc(n, nn, nkp_loc)) &
*cmplx(0.0_dp, -0.5_dp, kind=dp)
if (wann_slwf%constrain) then
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) + lambda_loc &
*crt(m, n)*(ln_tmp_loc(n, nn, nkp_loc) &
+ kmesh_info%wb(nn)*rnkb_loc(n, nn, nkp_loc)) &
*cmplx(0.0_dp, -0.5_dp, kind=dp) &
- lambda_loc*crt(m, n)*ln_tmp_loc(n, nn, nkp_loc) &
*cmplx(0.0_dp, -0.5_dp, kind=dp)
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) - kmesh_info%wb(nn) &
*lambda_loc &
*r0kb(n, nn, nkp_loc)*crt(m, n) &
*cmplx(0.0_dp, -0.5_dp, kind=dp)
end if
else
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc)
end if
end do
end do
else
do n = 1, num_wann
do m = 1, num_wann
! A[R^{k,b}]=(R-Rdag)/2
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) &
+ kmesh_info%wb(nn)*0.5_dp &
*(cr(m, n) - conjg(cr(n, m)))
! -S[T^{k,b}]=-(T+Tdag)/2i ; T_mn = Rt_mn q_n
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) - &
(crt(m, n)*ln_tmp_loc(n, nn, nkp_loc) &
+ conjg(crt(n, m)*ln_tmp_loc(m, nn, nkp_loc))) &
*cmplx(0.0_dp, -0.5_dp, kind=dp)
cdodq_loc(m, n, nkp_loc) = cdodq_loc(m, n, nkp_loc) - kmesh_info%wb(nn) &
*(crt(m, n)*rnkb_loc(n, nn, nkp_loc) &
+ conjg(crt(n, m)*rnkb_loc(m, nn, nkp_loc))) &
*cmplx(0.0_dp, -0.5_dp, kind=dp)
end do
end do
end if
end do
end do
cdodq_loc = cdodq_loc/real(num_kpts, dp)*4.0_dp
if (present(cdodq)) then
! each process communicates its result to other processes
cdodq(:, :, :) = 0.0_dp
do nkp_loc = 1, nkrank
nkp = global_k(nkp_loc)
cdodq(:, :, nkp) = cdodq_loc(:, :, nkp_loc)
end do
call comms_allreduce(cdodq(1, 1, 1), num_wann*num_wann*num_kpts, 'SUM', error, comm)
if (allocated(error)) return
if (lsitesymmetry) then
! correct behaviour is reproduced if algorithm 1 (mode 1) is followed by algorithm 2
call sitesym_symmetrize_gradient(sitesym, cdodq, 1, num_kpts, num_wann, error, comm)
call sitesym_symmetrize_gradient(sitesym, cdodq, 2, num_kpts, num_wann, error, comm)
do nkp_loc = 1, nkrank
nkp = global_k(nkp_loc)
cdodq_loc(:, :, nkp_loc) = cdodq(:, :, nkp)
end do
end if
end if
end if
deallocate (cr, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating cr in wann_domega', comm)
return
end if
deallocate (crt, stat=ierr)
if (ierr /= 0) then
call set_error_dealloc(error, 'Error in deallocating crt in wann_domega', comm)
return
end if
if (timing_level > 1 .and. iprint > 0) call io_stopwatch_stop('wann: domega', timer)
return
end subroutine wann_domega
