Sing_SearchOccDets Subroutine

    subroutine Sing_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 LoggingData, only: RDMExcitLevel
        use Parallel_neci, only: nProcessors, MPIArg
        use rdm_data, only: Sing_ExcDjs, Sing_ExcDjs2
        use rdm_data, only: one_rdms, two_rdm_spawn
        use rdm_filling, only: fill_sings_1rdm, fill_spawn_rdm_singles
        use searching, only: BinSearchParts_rdm
        use SystemData, only: nel
#ifdef DEBUG_
        character(*), parameter :: this_routine = "Sing_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), Ex_symm(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(Sing_ExcDjs2(:, StartDets), nI, SignDi, FlagsDi)

                do j = StartDets + 1, (NoDets + StartDets - 1)

                    ! D_i is in the first spot - start from the second.
                    iLutnJ(:) = Sing_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) = 1
                        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('Sing_SearchOccDets', 'nJ is not the correct excitation of nI.')

                        if (RDMExcitLevel == 1) then
                            call fill_sings_1rdm(one_rdms, Ex, tParity, SignDi, SignDj, .true.)
                        else
                            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 fill_spawn_rdm_singles(two_rdm_spawn, nI, Ex, full_sign)
                            ! Add in symmetric contribution.
                            Ex_symm(1, :) = Ex(2, :)
                            Ex_symm(2, :) = Ex(1, :)
                            call fill_spawn_rdm_singles(two_rdm_spawn, nI, Ex_symm, full_sign)
                        end if

                        ! No normalisation factor just yet - possibly need to revise.
                    end if

                end do
            end if
        end do

    end subroutine Sing_SearchOccDets