Add_RDM_From_IJ_Pair Subroutine

    subroutine Add_RDM_From_IJ_Pair(spawn, one_rdms, nI, nJ, realSignI, realSignJ)

        ! This routine takes a pair of different determinants Di and Dj, and
        ! figures out which type of elements need to be added in to the RDM.

        ! ------------------------ *IMPORTANT* --------------------------------
        ! nI refers to the ket determinant and nJ refers to the bra,
        ! determinant. Similarly, realSignI must be the sign corresponding to
        ! the ket determinant, and realSignJ to the bra, i.e. the element
        ! added in is equal to:
        !   realSignJ^* realSignI < nJ | ... | nI >
        ! Getting nI and nJ the wrong way around will result in RDM elements
        ! being added to the wrong side of the diagonal, which is a problem for
        ! transition RDMs which are *NOT* hermitian, so this is important.
        ! Similarly, getting nI and nJ right but realSignI and realSignJ the
        ! wrong way around will cause problems for complex RDMs, since the
        ! complex conjugate will be taken of the wrong sign.
        ! ---------------------------------------------------------------------

        use LoggingData, only: RDMExcitLevel
        use rdm_data, only: one_rdm_t
        use rdm_data_utils, only: add_to_rdm_spawn_t

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

        integer :: Ex(2, maxExcit), Ex_symm(2, maxExcit)
        logical :: tParity
        real(dp) :: full_sign(spawn%rdm_send%sign_length)
#ifdef DEBUG_
        character(*), parameter :: this_routine = "Add_RDM_From_IJ_Pair"
#endif
        Ex(:, :) = 0
        ! Maximum excitation level - we know they are connected by a double
        ! or a single excitation.
        Ex(1, 1) = 2
        tParity = .false.

        ! Ex(1,:) comes out as the orbital(s) excited from, i.e. i,j.
        ! Ex(2,:) comes out as the orbital(s) excited to, i.e. a,b.
        call GetExcitation(nI, nJ, nel, Ex, tParity)

        full_sign = 0.0_dp
        if (tParity) then
            full_sign = -realSignI * realSignJ
        else
            full_sign = realSignI * realSignJ
        end if

        if (any(ex < 0)) return

        if ((Ex(1, 2) == 0) .and. (Ex(2, 2) == 0)) then

            ! Di and Dj are separated by a single excitation.
            ! Add in the contribution from this pair into the 1-RDM.

            if (RDMExcitLevel == 1) then
                call fill_sings_1rdm(one_rdms, Ex, tParity, realSignI, realSignJ, .false.)
            else
                call fill_spawn_rdm_singles(spawn, nI, Ex, full_sign)
            end if

        else if (RDMExcitLevel /= 1) then

            ! Otherwise Di and Dj are connected by a double excitation.
            ! Add in this contribution to the 2-RDM, if being calculated.
            call add_to_rdm_spawn_t(spawn, Ex(2, 1), Ex(2, 2), Ex(1, 1), Ex(1, 2), full_sign, .false.)
        end if

    end subroutine Add_RDM_From_IJ_Pair