Sing_Hist_SearchOccDets Subroutine

    subroutine Sing_Hist_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 DetBitOps, only: FindBitExcitLevel
        use FciMCData, only: ilutHF_True, TotWalkers
        use hist, only: find_hist_coeff_explicit
        use hist_data, only: AllHistogram
        use LoggingData, only: RDMExcitLevel
        use Parallel_neci, only: nProcessors, MPIArg
        use rdm_data, only: Sing_ExcDjs, Sing_ExcDjs2, ExcNorm
        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

        integer(MPIArg), intent(in) :: recvcounts(nProcessors), recvdisps(nProcessors)
        integer(n_int) :: iLutnJ(0:NIfTot)
        integer, dimension(lenof_sign) :: HistPos
        real(dp), dimension(lenof_sign) :: RealHistPos

        integer :: i, j, NoDets, StartDets, PartInd, ExcitLevel
        integer :: nI(NEl), nJ(NEl), Ex(2, maxExcit), Ex_symm(2, maxExcit), FlagsDi, FlagsDj
        logical :: tDetFound, tParity
        real(dp) :: SignDi(lenof_sign), SignDj(lenof_sign), full_sign(1)

#ifdef DEBUG_
        character(*), parameter :: this_routine = "Sing_Hist_SearchOccDets"
#endif
        ! 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, RealHistPos, FlagsDi)

                HistPos = int(RealHistPos)

                SignDi = AllHistogram(1, HistPos(1)) / ExcNorm

                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)

                    ExcitLevel = FindBitExcitLevel(iLutHF_true, iLutnJ, NEl)
                    call find_hist_coeff_explicit(iLutnJ, ExcitLevel, 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 decode_bit_det(nJ, iLutnJ)

                        SignDj = AllHistogram(1, PartInd) / ExcNorm

                        ! 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_Hist_SearchOccDets