calc_2rdm_estimates_wrapper Subroutine

    subroutine calc_2rdm_estimates_wrapper(rdm_defs, est, rdm, en_pert)

        ! Calculate the estimates for the 2-RDM stored in rdm. The full estimates
        ! are stored using this object, and also instantaneous estimates. The
        ! instantaneous estimates are calculated by subtracting the previous
        ! stored values from the newly calculated ones. This works so long as
        ! the estimates are linear functions of the RDMs, which they will be
        ! for any observable.

        use CalcData, only: tEN2
        use Parallel_neci, only: MPISumAll
        use rdm_data, only: rdm_estimates_t, rdm_list_t, rdm_definitions_t
        use rdm_data, only: en_pert_t
        use SystemData, only: nel, ecore
        use OneEInts, only: PropCore
        use LoggingData, only: iNumPropToEst, tCalcPropEst

        type(rdm_definitions_t), intent(in) :: rdm_defs
        type(rdm_estimates_t), intent(inout) :: est
        type(rdm_list_t), intent(in) :: rdm
        type(en_pert_t), intent(in) :: en_pert

        integer :: irdm, iprop
        real(dp) :: rdm_trace(est%nrdms), rdm_norm(est%nrdms)
        real(dp) :: rdm_energy_1(est%nrdms), rdm_energy_2(est%nrdms)
        real(dp) :: rdm_spin(est%nrdms)
        real(dp) :: rdm_prop(iNumProptoEst, est%nrdms)
        real(dp) :: energy_pert(est%nrdms_standard)

        ! Use the _inst variables as temporary variables to store the current
        ! total values. These are updated at the end of this routine.
        est%trace_inst = est%trace
        est%norm_inst = est%norm
        est%energy_1_num_inst = est%energy_1_num
        est%energy_2_num_inst = est%energy_2_num
        est%energy_num_inst = est%energy_num
        if (.not. tGUGA) est%spin_num_inst = est%spin_num
        if (tCalcPropEst) est%property_inst = est%property

        ! Calculate the new total values.

        call calc_rdm_trace(rdm, rdm_trace)
        call MPISumAll(rdm_trace, est%trace)

        ! RDMs are normalised so that their trace is nel*(nel-1)/2.
        if (tGUGA) then
            rdm_norm = rdm_trace * 1.0_dp / (nel * (nel - 1))
            est%norm = est%trace * 1.0_dp / (nel * (nel - 1))
        else
            rdm_norm = rdm_trace * 2.0_dp / (nel * (nel - 1))
            est%norm = est%trace * 2.0_dp / (nel * (nel - 1))
        end if

        ! For the transition RDMs, we want to calculate the norms using the
        ! non-transition RDMs.
        do irdm = est%nrdms_standard + 1, est%nrdms
            est%norm(irdm) = sqrt(est%norm(rdm_defs%state_labels(1, irdm)) * est%norm(rdm_defs%state_labels(2, irdm)))
        end do

        ! The 1- and 2- electron operator contributions to the RDM energy.
        if (tGUGA) then
            call calc_rdm_energy_guga(rdm, rdm_energy_1, rdm_energy_2)
        else
            call calc_rdm_energy(rdm, rdm_energy_1, rdm_energy_2)
        end if

        call MPISumAll(rdm_energy_1, est%energy_1_num)
        call MPISumAll(rdm_energy_2, est%energy_2_num)
        ! The *total* energy, including the core contribution.
        est%energy_num = est%energy_1_num + est%energy_2_num + ecore * est%norm

        ! Estimate of the expectation value of the spin squared operator
        ! (equal to S(S+1) for spin quantum number S).
        if (.not. tGUGA) then
            call calc_rdm_spin(rdm, rdm_norm, rdm_spin)
            call MPISumAll(rdm_spin, est%spin_num)
        end if

        if (tCalcPropEst) then
            ! Estimate of the properties using different property integrals
            ! and all the standard and transition rdms.
            call calc_rdm_prop(rdm, rdm_prop)
            call MPISumAll(rdm_prop, est%property)
            ! Add the contribution from the core (zero body) part of the perturbation.
            do iprop = 1, iNumPropToEst
                est%property(iprop, :) = est%property(iprop, :) + est%norm * PropCore(iprop)
            end do
        end if

        ! Calculate the instantaneous values by subtracting the old total
        ! values from the new total ones.
        est%trace_inst = est%trace - est%trace_inst
        est%norm_inst = est%norm - est%norm_inst
        est%energy_1_num_inst = est%energy_1_num - est%energy_1_num_inst
        est%energy_2_num_inst = est%energy_2_num - est%energy_2_num_inst
        est%energy_num_inst = est%energy_num - est%energy_num_inst
        if (.not. tGUGA) est%spin_num_inst = est%spin_num - est%spin_num_inst
        if (tCalcPropEst) est%property_inst = est%property - est%property_inst

        ! For the EN Perturbation terms, we clear them at the start of
        ! every RDM averaging cycle, so they're treated a bit differently.
        if (tEN2) then
            call calc_en_pert_energy(en_pert, est%energy_num_inst(1:en_pert%sign_length), &
                                     est%norm_inst(1:en_pert%sign_length), energy_pert)

            call MPISumAll(energy_pert, est%energy_pert_inst)

            est%energy_pert = est%energy_pert + est%energy_pert_inst
        end if

    end subroutine calc_2rdm_estimates_wrapper