subroutine plot_write_r2mn(num_kpts, num_wann, kmesh_info, m_matrix, seedname, dist_k, error, comm)
!================================================!
!
! Write seedname.r2mn file
!
!================================================!
use w90_comms, only: w90_comm_type
use w90_constants, only: dp
use w90_error, only: w90_error_type, set_error_file
use w90_types, only: kmesh_info_type
implicit none
type(kmesh_info_type), intent(in) :: kmesh_info
type(w90_error_type), allocatable, intent(out) :: error
type(w90_comm_type), intent(in) :: comm
integer, intent(in) :: num_kpts, num_wann
integer, intent(in) :: dist_k(:) ! MPI k-point distribution
complex(kind=dp), intent(in) :: m_matrix(:, :, :, :)
character(len=50), intent(in) :: seedname
integer :: r2mnunit, nw1, nw2, nkp, nn, ierr
integer :: nkp_rank, my_node_id
real(kind=dp) :: r2ave_mn, delta
logical :: on_root = .false.
! note that here I use formulas analogue to Eq. 23, and not to the
! shift-invariant Eq. 32 .
my_node_id = mpirank(comm)
if (my_node_id == 0) on_root = .true.
if (on_root) then
open (newunit=r2mnunit, file=trim(seedname)//'.r2mn', form='formatted', iostat=ierr)
if (ierr /= 0) then
call set_error_file(error, 'Error opening file '//trim(seedname)//'.r2mn in plot_write_r2mn', comm)
return
end if
end if
do nw1 = 1, num_wann
do nw2 = 1, num_wann
r2ave_mn = 0.0_dp
delta = 0.0_dp
if (nw1 .eq. nw2) delta = 1.0_dp
nkp_rank = 1
do nkp = 1, num_kpts
if (dist_k(nkp) /= my_node_id) cycle
do nn = 1, kmesh_info%nntot
! [GP-begin, Apr13, 2012: corrected sign inside "real"]
r2ave_mn = r2ave_mn + kmesh_info%wb(nn)* &
(2.0_dp*delta - real(m_matrix(nw1, nw2, nn, nkp_rank) + &
conjg(m_matrix(nw2, nw1, nn, nkp_rank)), dp))
end do
nkp_rank = nkp_rank + 1
end do ! global k list
call comms_reduce(r2ave_mn, 1, 'SUM', error, comm)
r2ave_mn = r2ave_mn/real(num_kpts, dp)
if (on_root) write (r2mnunit, '(2i6,f20.12)') nw1, nw2, r2ave_mn
end do
end do
if (on_root) close (r2mnunit)
end subroutine plot_write_r2mn