store_parent_with_spawned Subroutine

    subroutine store_parent_with_spawned(RDMBiasFacCurr, WalkerNumber, iLutI, &
                                         DetSpawningAttempts, iLutJ, procJ)

        ! We are spawning from iLutI to SpawnedParts(:,ValidSpawnedList(proc)).
        ! This routine stores the parent (D_i) with the spawned child (D_j) so
        ! that we can add in Ci.Cj to the RDM later on. The parent is nifd
        ! integers long, and stored in the second part of the SpawnedParts array
        ! from NIfTot+1 -> NIfTot+1 + nifd.

        use DetBitOps, only: DetBitEQ
        use FciMCData, only: SpawnedParts, ValidSpawnedList, TempSpawnedParts, &
                             TempSpawnedPartsInd
        use CalcData, only: tNonInitsForRDMs

        real(dp), intent(in) :: RDMBiasFacCurr
        integer, intent(in) :: WalkerNumber, procJ
        integer, intent(in) :: DetSpawningAttempts
        integer(n_int), intent(in) :: iLutI(0:niftot), iLutJ(0:niftot)
        logical :: tRDMStoreParent
        integer :: j

        if (abs(RDMBiasFacCurr) < 1.0e-12_dp) then
            ! If RDMBiasFacCurr is exactly zero, any contribution from Ci.Cj
            ! will be zero so it is not worth carrying on.
            call zero_parent(SpawnedParts(:, ValidSpawnedList(procJ)))
            ! in the case of non-variational rdms, we also require the parent
            ! to be an initiator -> only take non-initiators when the
            ! init-agnostic rdms are requested (tNonInitsForRDMs)
            return
        end if

        if (tNonInitsForRDMs .or. all_runs_are_initiator(ilutI)) then

            ! First we want to check if this Di.Dj pair has already been accounted for.
            ! This means searching the Dj's that have already been spawned from this Di, to make sure
            ! the new Dj being spawned on here is not the same.
            ! The Dj children spawned by the current Di are being stored in the array TempSpawnedParts,
            ! so that the reaccurance of a Di.Dj pair may be monitored.

            ! Store the Di parent with the spawned child, unless we find this Dj has already been spawned on.
            tRDMStoreParent = .true.

            ! Run through the Dj walkers that have already been spawned from this particular Di.
            ! If this is the first to be spawned from Di, TempSpawnedPartsInd will be zero, so we
            ! just wont run over anything.

            do j = 1, TempSpawnedPartsInd
                if (DetBitEQ(iLutJ(0:nifd), TempSpawnedParts(0:nifd, j), nifd)) then
                    ! If this Dj is found, we do not want to store the parent with this spawned walker.
                    tRDMStoreParent = .false.
                    exit
                end if
            end do

            if (tRDMStoreParent) then
                ! This is a new Dj that has been spawned from this Di.
                ! We want to store it in the temporary list of spawned parts which have come from this Di.
                if (WalkerNumber /= DetSpawningAttempts) then
                    ! Don't bother storing these if we're on the last walker, or if we only have one
                    ! walker on Di.
                    TempSpawnedPartsInd = TempSpawnedPartsInd + 1
                    TempSpawnedParts(0:nifd, TempSpawnedPartsInd) = iLutJ(0:nifd)
                end if

                ! We also want to make sure the parent Di is stored with this Dj.

                ! We need to carry with the child (and the parent), the sign of the parent.
                ! In actual fact this is the sign of the parent divided by the probability of generating
                ! that pair Di and Dj, to account for the
                ! fact that Di and Dj are not always added to the RDM, but only when Di spawns on Dj.
                ! This RDMBiasFacCurr factor is turned into an integer to pass around to the relevant processors.
                call encode_parent(SpawnedParts(:, ValidSpawnedList(procJ)), &
                                   ilutI, RDMBiasFacCurr)

            else
                ! This Di has already spawned on this Dj - don't store the Di parent with this child,
                ! so that the pair is not double counted.
                ! We are using the probability that Di spawns onto Dj *at least once*, so we don't want to
                ! double count this pair.
                call zero_parent(SpawnedParts(:, ValidSpawnedList(procJ)))
            end if
        else
            ! This Di is a non-initiator, and we specified to take only initiators into account
            call zero_parent(SpawnedParts(:, ValidSpawnedList(procJ)))
        end if

    end subroutine store_parent_with_spawned