subroutine fill_diag_1rdm_guga(one_rdms, nI, contrib_sign, &
tCoreSpaceDetIn, RDMItersIn)
! I think I have to change these routines to also take into
! account the possible coupling coefficient of doubly occupied
! orbitals
type(one_rdm_t), intent(inout) :: one_rdms(:)
integer, intent(in) :: nI(nel)
real(dp), intent(in) :: contrib_sign(:)
logical, optional, intent(in) :: tCoreSpaceDetIn
integer, optional, intent(in) :: RDMItersIn(:)
integer :: i, ind, irdm, s, s_orbs(nel), inc
integer :: RDMIters(size(one_rdms))
real(dp) :: ScaleContribFac, final_contrib(size(one_rdms))
logical :: tCoreSpaceDet
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
s_orbs = gtID(nI)
i = 1
do while (i <= nel)
s = s_orbs(i)
ind = SymLabelListInv_rot(s)
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
if (isDouble(nI, i)) then
inc = 2
final_contrib = 2.0_dp * final_contrib
else
inc = 1
end if
if (tAdaptiveShift .and. all(nI == projEDet(:, 1))) &
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
i = i + inc
end do
end subroutine fill_diag_1rdm_guga