Fill_Spin_Coupled_RDM Subroutine

    subroutine Fill_Spin_Coupled_RDM(spawn, one_rdms, iLutnI, iLutnJ, nI, nJ, realSignI, realSignJ)

        ! It takes to HPHF functions, and calculates what needs to be summed
        ! into the RDMs.

        ! If the two HPHF determinants we're considering consist of I + I' and
        ! J + J', where X' is the spin coupled (all spins flipped) version of X,
        ! then we have already considered the I -> J excitation. And if I and
        ! J are connected by a double excitation, tDoubleConnection is true
        ! and we have also considered I' -> J'. But we need to also account
        ! for I -> J' and I' -> J.

        use DetBitOps, only: TestClosedShellDet, FindBitExcitLevel
        use hphf_integrals, only: hphf_sign
        use HPHFRandExcitMod, only: FindExcitBitDetSym, FindDetSpinSym
        use rdm_data, only: one_rdm_t
        use SystemData, only: nel

        type(rdm_spawn_t), intent(inout) :: spawn
        type(one_rdm_t), intent(inout) :: one_rdms(:)
        integer(n_int), intent(in) :: iLutnI(0:NIfTot), iLutnJ(0:NIfTot)
        real(dp), intent(in) :: realSignI(:), realSignJ(:)
        integer, intent(in) :: nI(nel), nJ(nel)

        integer(n_int) :: iLutnI2(0:NIfTot)
        integer :: nI2(nel), nJ2(nel)
        real(dp) :: NewSignJ(size(realSignJ)), NewSignI(size(realSignI))
        real(dp) :: PermSignJ(size(realSignJ)), PermSignI(size(realSignI))
        integer :: I_J_ExcLevel, ICoup_J_ExcLevel

        if (TestClosedShellDet(iLutnI)) then

            if (TestClosedShellDet(iLutnJ)) then
                ! Closed shell -> Closed shell - just as in determinant case
                call Add_RDM_From_IJ_Pair(spawn, one_rdms, nI, nJ, realSignI, realSignJ)
            else
                ! Closed shell -> open shell.
                call FindDetSpinSym(nJ, nJ2, nel)
                NewSignJ = realSignJ / Root2
                call Add_RDM_From_IJ_Pair(spawn, one_rdms, nI, nJ, realSignI, NewSignJ)
                ! What is the permutation between Di and Dj'
                NewSignJ = NewSignJ * hphf_sign(iLutnJ)
                call Add_RDM_From_IJ_Pair(spawn, one_rdms, nI, nJ2, realSignI, NewSignJ)
            end if

        else if (TestClosedShellDet(iLutnJ)) then
            ! Open shell -> closed shell
            call FindDetSpinSym(nI, nI2, nel)
            NewSignI = realSignI / Root2
            call Add_RDM_From_IJ_Pair(spawn, one_rdms, nI, nJ, NewSignI, realSignJ)
            ! What is the permutation between Di' and Dj?
            NewSignI = NewSignI * hphf_sign(iLutnI)
            call Add_RDM_From_IJ_Pair(spawn, one_rdms, nI2, nJ, NewSignI, realSignJ)

        else
            ! Open shell -> open shell
            NewSignI = realSignI / Root2
            NewSignJ = realSignJ / Root2
            PermSignJ = NewSignJ * real(hphf_sign(iLutnJ), dp)
            PermSignI = NewSignI * real(hphf_sign(iLutnI), dp)
            call FindExcitBitDetSym(iLutnI, iLutnI2)
            call FindDetSpinSym(nI, nI2, nel)
            call FindDetSpinSym(nJ, nJ2, nel)
            I_J_ExcLevel = FindBitExcitLevel(iLutnI, iLutnJ, 2)
            ICoup_J_ExcLevel = FindBitExcitLevel(iLutnI2, iLutnJ, 2)

            if (I_J_ExcLevel <= 2) then
                call Add_RDM_From_IJ_Pair(spawn, one_rdms, nI, nJ, NewSignI, NewSignJ)
                ! Di -> Dj
                call Add_RDM_From_IJ_Pair(spawn, one_rdms, nI2, nJ2, PermSignI, PermSignJ)
                ! Di' -> Dj'  (both permuted sign)
            end if

            if (ICoup_J_ExcLevel <= 2) then
                call Add_RDM_From_IJ_Pair(spawn, one_rdms, nI2, nJ, PermSignI, NewSignJ)
                ! Di' -> Dj  (i permuted sign)
                call Add_RDM_From_IJ_Pair(spawn, one_rdms, nI, nJ2, NewSignI, PermSignJ)
                ! Di  -> Dj'  (j permuted sign)
            end if
        end if

    end subroutine Fill_Spin_Coupled_RDM