Doub_SearchOccDets Subroutine

    subroutine Doub_SearchOccDets(recvcounts, recvdisps)

        ! We now have arrays SingExcDjs2 which contain all the single excitations
        ! from each processor. These number sent from processor i is
        ! recvcounts(i), and the first 2 have information  about the determinant
        ! Di from which the Dj's are single excitations (and it's sign).

        use FciMCData, only: CurrentDets, TotWalkers
        use Parallel_neci, only: nProcessors, MPIArg
        use rdm_data, only: Doub_ExcDjs, Doub_ExcDjs2, two_rdm_spawn
        use rdm_data_utils, only: add_to_rdm_spawn_t
        use searching, only: BinSearchParts_rdm
        use SystemData, only: nel
#ifdef DEBUG_
        character(*), parameter :: this_routine = "Doub_SearchOccDets"
#endif
        integer(MPIArg), intent(in) :: recvcounts(nProcessors), recvdisps(nProcessors)
        integer(n_int) :: iLutnJ(0:NIfTot)
        real(dp) :: SignDi(lenof_sign), SignDj(lenof_sign), full_sign(1)
        integer :: i, j, NoDets, StartDets, PartInd
        integer :: nI(NEl), nJ(NEl), Ex(2, maxExcit), FlagsDi, FlagsDj
        logical :: tDetFound, tParity

        ! Take each Dj, and binary search CurrentDets to see if it is occupied.
        do i = 1, nProcessors
            ! Doing determinants from each processor separately because each
            ! has a different D_i it was excited from.
            NoDets = recvcounts(i) / (NIfTot + 1)
            StartDets = (recvdisps(i) / (NIfTot + 1)) + 1

            if (NoDets > 1) then
                call extract_bit_rep(Doub_ExcDjs2(:, StartDets), nI, SignDi, FlagsDi)

                do j = StartDets + 1, (NoDets + StartDets - 1)
                    ! D_i is in the first spot - start from the second.
                    iLutnJ(:) = Doub_ExcDjs2(:, j)

                    ! This binary searches CurrentDets between 1 and TotWalkers
                    ! for determinant iLutnJ. If found, tDetFound will be true,
                    ! and PartInd the index in CurrentDets where the determinant is.
                    call BinSearchParts_rdm(iLutnJ, 1, int(TotWalkers), PartInd, tDetFound)

                    if (tDetFound) then
                        ! Determinant occupied; add c_i*c_j to the relevant
                        ! element of nElRDM. Need to first find the orbitals
                        ! involved in the excitation from D_i -> D_j and
                        ! the parity.
                        Ex(:, :) = 0
                        ! Ex(1,1) goes in as the max number of excitations -
                        ! we know this is an excitation of level RDMExcitLevel.
                        Ex(1, 1) = 2
                        tParity = .false.

                        call extract_bit_rep(CurrentDets(:, PartInd), nJ, SignDj, FlagsDj)

                        ! Ex(1,:) comes out as the orbital(s) excited from,
                        ! Ex(2,:) comes out as the orbital(s) excited to.
                        call GetExcitation(nI, nJ, NEl, Ex, tParity)

                        if (Ex(1, 1) <= 0) call Stop_All('SearchOccDets', &
                                                         'nJ is not the correct excitation of nI.')

                        ASSERT(.not. t_3_body_excits)
                        if (tParity) then
                            full_sign = -SignDi(1) * SignDj(lenof_sign)
                        else
                            full_sign = SignDi(1) * SignDj(lenof_sign)
                        end if

                        call add_to_rdm_spawn_t(two_rdm_spawn, Ex(2, 1), Ex(2, 2), &
                            Ex(1, 1), Ex(1, 2), full_sign, .false.)
                        ! Add in symmetric contribution.
                        call add_to_rdm_spawn_t(two_rdm_spawn, Ex(1, 1), Ex(1, 2), &
                            Ex(2, 1), Ex(2, 2), full_sign, .false.)
                    end if
                end do
            end if
        end do

    end subroutine Doub_SearchOccDets