subroutine calc_double_occ_from_rdm(rdm, rdm_trace, inst_occ)
! also write a routine which calculates the double occupancy from the
! 2-rdm, if it has been calculated!
type(rdm_list_t), intent(inout) :: rdm
real(dp), intent(in) :: rdm_trace(rdm%sign_length)
real(dp), intent(out), optional :: inst_occ
character(*), parameter :: this_routine = "calc_double_occ_from_rdm"
integer :: ielem, ij, kl, i, j, k, l, p, q, r, s, iproc, irdm, ierr, hash_val
integer(int_rdm) :: ijkl
real(dp) :: rdm_sign(rdm%sign_length)
real(dp) :: double_occ(rdm%sign_length), all_double_occ(rdm%sign_length)
real(dp), allocatable :: spatial_double_occ(:), all_spatial_double_occ(:)
logical :: tSuccess
double_occ = 0.0_dp
! just a quick addition to calculate spatially resolved double
! occupancy
if (t_spin_measurements) then
allocate(spatial_double_occ(nBasis / 2))
allocate(all_spatial_double_occ(nBasis / 2))
spatial_double_occ = 0.0_dp
all_spatial_double_occ = 0.0_dp
end if
! todo: find out about the flags to ensure the rdm was actually
! calculated!
! seperately on all processors do this summation and then
! communicate the results in the end..
do ielem = 1, rdm%nelements
ijkl = rdm%elements(0, ielem)
call calc_separate_rdm_labels(ijkl, ij, kl, i, j, k, l)
call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign)
! normalise
rdm_sign = rdm_sign / rdm_trace
! convert to spatial orbitals:
p = spatial(i)
q = spatial(j)
r = spatial(k)
s = spatial(l)
! only consider the diagonal elements!
if (p == q .and. p == r .and. p == s) then
ASSERT(.not. same_spin(i, j))
ASSERT(.not. same_spin(k, l))
! add up all the diagonal contributions
double_occ = double_occ + rdm_sign
if (t_spin_measurements) spatial_double_occ(p) = rdm_sign(1)
end if
end do
! at the end average over the spatial orbitals
double_occ = 2.0_dp * double_occ / real(nbasis, dp)
! MPI communicate:
call MPISumAll(double_occ, all_double_occ)
if (present(inst_occ)) then
inst_occ = double_occ(1)
end if
if (t_spin_measurements) then
call MPIAllreduce(spatial_double_occ, MPI_SUM, all_spatial_double_occ)
end if
if (iProcIndex == root) then
print *, "======"
print *, "Double occupancy from RDM: ", all_double_occ
print *, "======"
if (t_spin_measurements) then
print *, "======"
print *, "spatially resolved double occupancy from RDM: "
print *, all_spatial_double_occ
print *, "======"
end if
end if
end subroutine calc_double_occ_from_rdm