fill_diag_1rdm Subroutine

    subroutine fill_diag_1rdm(one_rdms, nI, contrib_sign, tCoreSpaceDetIn, RDMItersIn, tLagrCorr)

        ! Add the contribution to the diagonal elements of the 1-RDM from
        ! determinant nI, with the corresponding sign(s) from an FCIQMC
        ! simulation gievn by contrib_sign.

        ! The (spinned) 1-RDM is defined by
        !
        ! \gamma_{i,j} = < \Psi | a^{\dagger}_i a_j | \Psi >.
        !
        ! If we have a closed shell system then we instead add to the spin-free
        ! 1-RDM, obtained by tracing over the spin component.

        use rdm_data, only: one_rdm_t, tOpenShell
        use LoggingData, only: ThreshOccRDM, tThreshOccRDMDiag, tExplicitAllRDM
        use RotateOrbsData, only: SymLabelListInv_rot
        use UMatCache, only: gtID

        type(one_rdm_t), intent(inout) :: one_rdms(:)
        integer, intent(in) :: nI(:)
        real(dp), intent(in) :: contrib_sign(:)
        logical, optional, intent(in) :: tCoreSpaceDetIn
        integer, optional, intent(in) :: RDMItersIn(:)
        logical, optional, intent(in) :: tLagrCorr

        integer :: i, ind, irdm
        integer :: RDMIters(size(one_rdms))
        real(dp) :: ScaleContribFac, final_contrib(size(one_rdms))
        logical :: tCoreSpaceDet, tLC

        if (present(tLagrCorr)) then
            tLC = tLagrCorr
        else
            tLC = .true.
        end if

        ScaleContribFac = 1.0_dp

        RDMIters = 1
        if (present(RDMItersIn)) RDMIters = RDMItersIn

        tCoreSpaceDet = .false.
        if (present(tCoreSpaceDetIn)) tCoreSpaceDet = tCoreSpaceDetIn

        ! This is the single-run cutoff being applied (do not use in DR mode):
        if (.not. tCoreSpaceDet) then
            ! Dets in the core space are never removed from main list, so
            ! strictly do not require corrections.
            if (tThreshOccRDMDiag .and. (abs(contrib_sign(1)) <= ThreshOccRDM)) ScaleContribFac = 0.0_dp
        end if

        do i = 1, size(nI)
            ! The SymLabelListInv_rot array is used to index the 1-RDM so that
            ! it will be filled in block diagonal order, where each block holds
            ! one symmetry sector of the 1-RDM (i.e., orbitals are ordered
            ! first by symmetry, then by the standard order).
            if (tOpenShell) then
                ind = SymLabelListInv_rot(nI(i))
            else
                ind = SymLabelListInv_rot(gtID(nI(i)))
            end if
            if (size(contrib_sign) == 1) then
                final_contrib = contrib_sign**2 * RDMIters * ScaleContribFac
            else
                final_contrib = contrib_sign(1::2)*contrib_sign(2::2) * RDMIters * ScaleContribFac
            end if
            ! in adaptive shift mode, the reference contribution is rescaled
            ! we assume that projEDet is the same on all runs, else there is no point
            if (tAdaptiveShift .and. all(nI == projEDet(:, 1)) &
                .and. tNonInitsForRDMs .and. .not. tInitsRDMRef .and. tLC) &
                final_contrib = final_contrib + RDMIters * ScaleContribFac * rdmCorrectionFactor
            do irdm = 1, size(one_rdms)
                one_rdms(irdm)%matrix(ind, ind) = one_rdms(irdm)%matrix(ind, ind) + final_contrib(irdm)
            end do
        end do

    end subroutine fill_diag_1rdm