Add_Hist_ExplicitRDM_Contrib Subroutine

    subroutine Add_Hist_ExplicitRDM_Contrib(iLutnI, blank_det)

        ! This is the general routine for taking a particular determinant in
        ! the spawned list, D_i and adding it's contribution to the reduced
        ! density matrix.

        use LoggingData, only: RDMExcitLevel
        use Parallel_neci, only: nProcessors
        use rdm_data, only: Sing_ExcList, Doub_ExcList, Sing_InitExcSlots, Doub_InitExcSlots
        use rdm_data, only: Sing_ExcDjs, Doub_ExcDjs

        integer(n_int), intent(in) :: iLutnI(0:NIfTot)
        logical, intent(in) :: blank_det
        integer :: i

        ! Set up excitation arrays.
        ! These are blocked according to the processor the excitation would be
        ! on if occupied. In each block, the first entry is the sign of
        ! determinant D_i and the second the bit string of the determinant
        ! (these need to be sent along with the excitations). Each processor
        ! will have a different Di.

        Sing_ExcDjs(:, :) = 0
        Sing_ExcList(:) = 0
        Sing_ExcList(:) = Sing_InitExcSlots(:)

        do i = 0, nProcessors - 1
            Sing_ExcDjs(:, Sing_ExcList(i)) = iLutnI(:)
            Sing_ExcList(i) = Sing_ExcList(i) + 1
        end do
        if (RDMExcitLevel /= 1) then
            Doub_ExcDjs(:, :) = 0
            Doub_ExcList(:) = 0
            Doub_ExcList(:) = Doub_InitExcSlots(:)

            do i = 0, nProcessors - 1
                Doub_ExcDjs(:, Doub_ExcList(i)) = iLutnI(:)
                Doub_ExcList(i) = Doub_ExcList(i) + 1
            end do
        end if

        ! Out of here we will get a filled ExcDjs array with all the single or
        ! double excitations  from Dj, this will be done for each proc.
        if (.not. blank_det) call Gen_Hist_ExcDjs(iLutnI)

        ! We then need to send the excitations to the relevant processors.
        ! This routine then calls SearchOccDets which takes each excitation
        ! and and binary searches the occupied determinants for this. If found,
        ! we re-find the orbitals and parity involved in the excitation, and add
        ! the c_i*c_j contributions to the corresponding matrix element.
        call Send_Hist_ProcExcDjs()

    end subroutine Add_Hist_ExplicitRDM_Contrib