calc_double_occ_from_rdm Subroutine

    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