CompressSpawnedList Subroutine

public subroutine CompressSpawnedList(ValidSpawned, iter_data)


Type IntentOptional Attributes Name
integer :: ValidSpawned
type(fcimc_iter_data), intent(inout) :: iter_data


Source Code

Source Code

    subroutine CompressSpawnedList(ValidSpawned, iter_data)

        ! This sorts and compresses the spawned list to make it easier for the
        ! rest of the annihilation process. This is not essential, but should
        ! prove worthwhile.

        type(fcimc_iter_data), intent(inout) :: iter_data
        integer :: VecInd, ValidSpawned, DetsMerged, i, BeginningBlockDet, FirstInitIndex, CurrentBlockDet
        real(dp) :: SpawnedSign(lenof_sign), Temp_Sign(lenof_sign)
        integer :: EndBlockDet, part_type, Parent_Array_Ind
        integer :: No_Spawned_Parents
        integer(kind=n_int), pointer :: PointTemp(:, :)
        integer(n_int) :: cum_det(0:IlutBits%len_bcast), temp_det(0:IlutBits%len_bcast)
        character(len=*), parameter :: t_r = 'CompressSpawnedList'
        type(timer), save :: Sort_time
        ! We want to sort the list of newly spawned particles, in order for
        ! quicker binary searching later on. They should remain sorted after
        ! annihilation between spawned.

        if (.not. bNodeRoot) return

        Sort_time%timer_name = 'Compress Sort interface'
        call set_timer(Sort_time, 20)

        call sort(SpawnedParts(0:IlutBits%len_bcast, 1:ValidSpawned), ilut_lt, ilut_gt)

        call halt_timer(Sort_time)

        if (tHistSpawn) HistMinInd2(1:NEl) = FCIDetIndex(1:NEl)

        ! First, we compress the list of spawned particles, so that they are
        ! only specified at most once in each processors list. During this, we
        ! transfer the particles from SpawnedParts to SpawnedParts2. If we are
        ! working with complex walkers, we essentially do the same thing twice,
        ! annihilating real and imaginary particles seperately.

        ! This is the index in the SpawnedParts2 array to copy the compressed
        ! walkers into.
        VecInd = 1
        ! BeginningBlockDet will indicate the index of the first entry for a
        ! given determinant in SpawnedParts.
        BeginningBlockDet = 1
        DetsMerged = 0
        Parent_Array_Ind = 1
        Spawned_Parts_Zero = 0

        do while (BeginningBlockDet <= ValidSpawned)

            ! Loop in blocks of the same determinant to the end of the list of
            ! walkers.

            FirstInitIndex = 0
            CurrentBlockDet = BeginningBlockDet + 1

            do while (CurrentBlockDet <= ValidSpawned)
                if (.not. (DetBitEQ(SpawnedParts(:, BeginningBlockDet), &
                                    SpawnedParts(:, CurrentBlockDet)))) exit
                ! Loop over walkers on the same determinant in SpawnedParts.
                CurrentBlockDet = CurrentBlockDet + 1
            end do

            ! EndBlockDet indicates that we have reached the end of the block
            ! of similar dets.
            EndBlockDet = CurrentBlockDet - 1

            if (EndBlockDet == BeginningBlockDet) then
                ! Optimisation: This block only consists of one entry. Simply
                !               copy it across rather than explicitly searching
                !               the list.

                ! If this one entry has no amplitude then don't add it to the
                ! compressed list, but just cycle.
                call extract_sign(SpawnedParts(:, BeginningBlockDet), temp_sign)
                if ((sum(abs(temp_sign)) < 1.e-12_dp) .and. &
                    (.not. (tFillingStochRDMonFly .and. (.not. tNoNewRDMContrib)))) then
                    DetsMerged = DetsMerged + 1
                    BeginningBlockDet = CurrentBlockDet
                end if

                ! Transfer all info to the other array.
                SpawnedParts2(:, VecInd) = SpawnedParts(:, BeginningBlockDet)

                if (tFillingStochRDMonFly .and. (.not. tNoNewRDMContrib)) then
                    ! SpawnedParts contains the determinants spawned on (Dj),
                    ! and it's parent (Di) plus it's sign (Cj).
                    ! As in | Dj | Di | Ci |
                    ! We then compress multiple occurances of Dj, but these may
                    ! have come from different parents, and we want to keep
                    ! track of all Di's. As we compress SpawnedParts, we
                    ! therefore move all the parents (Di's) into Spawned_Parents.
                    ! If the compressed Dj is at position VecInd in SpawnedParts,
                    ! then Spawned_Parents_Index(1,VecInd) is the starting point
                    ! of it's parents (Di) in Spawned_Parents, and there are
                    ! Spawned_Parents_Index(2,VecInd) entries corresponding to
                    ! this Dj.

                    if (.not. bit_parent_zero(SpawnedParts(:, BeginningBlockDet))) then

                        ! If the parent determinant is null, the contribution to
                        ! the RDM is zero. No point in doing anything more with it.

                        Spawned_Parents(0:IlutBitsParent%ind_flag, Parent_Array_Ind) = &
                            SpawnedParts(IlutBits%ind_parent:IlutBits%ind_parent_flag, BeginningBlockDet)

                        call extract_sign(SpawnedParts(:, BeginningBlockDet), temp_sign)

                        ! Search to see which sign is non-zero, and therefore
                        ! find which simulation the spawning occured from and to.
                        ! NOTE: it is safe to compare against zero exactly here,
                        ! because all other components will have been set to zero
                        ! exactly and can't have changed at all.
                        Spawned_Parents(IlutBitsParent%ind_source, Parent_Array_Ind) = 0
                        do part_type = 1, lenof_sign
                            if (abs(temp_sign(part_type)) > 1.0e-12_dp) then
                                Spawned_Parents(IlutBitsParent%ind_source, Parent_Array_Ind) &
                                    = part_type
                            end if
                        end do

                        ! in the guga case we also need to transfer the rdm information
                        if (tGUGA) then
                            call transfer_stochastic_rdm_info( &
                                SpawnedParts(:, BeginningBlockDet), &
                                Spawned_Parents(:, Parent_Array_Ind), &
                                BitIndex_from=IlutBits, &
                        end if

                        ! The first nifd of the Spawned_Parents entry is the
                        ! parent determinant, the nifd + 1 entry is the Ci.
                        ! Parent_Array_Ind keeps track of the position in
                        ! Spawned_Parents.
                        Spawned_Parents_Index(1, VecInd) = Parent_Array_Ind
                        Spawned_Parents_Index(2, VecInd) = 1
                        ! In this case there is only one instance of Dj - so
                        ! therefore only 1 parent Di.
                        Parent_Array_Ind = Parent_Array_Ind + 1
                        Spawned_Parents_Index(1, VecInd) = Parent_Array_Ind
                        Spawned_Parents_Index(2, VecInd) = 0
                    end if

                    call extract_sign(SpawnedParts(:, BeginningBlockDet), temp_sign)
                    if (IsUnoccDet(temp_sign)) then
                        Spawned_Parts_Zero = Spawned_Parts_Zero + 1
                    end if
                end if

                VecInd = VecInd + 1
                ! Move onto the next block of determinants.
                BeginningBlockDet = CurrentBlockDet
                cycle ! Skip the rest of this block.
            end if

            ! Reset the cumulative determinant
            cum_det = 0_n_int
            cum_det(0:nifd) = SpawnedParts(0:nifd, BeginningBlockDet)

            if (tPreCond .or. tReplicaEstimates) then
                cum_det(IlutBits%ind_hdiag) = SpawnedParts(IlutBits%ind_hdiag, BeginningBlockDet)
            end if

            if (tFillingStochRDMonFly .and. (.not. tNoNewRDMContrib)) then
                ! This is the first Dj determinant - set the index for the
                ! beginning of where the parents for this Dj can be found in
                ! Spawned_Parents.
                Spawned_Parents_Index(1, VecInd) = Parent_Array_Ind

                ! In this case, multiple Dj's must be compressed, and therefore
                ! the Di's dealt with as  described above. We first just
                ! initialise the position in the Spawned_Parents array to enter
                ! the Di's.
                Spawned_Parents_Index(2, VecInd) = 0
            end if

            do i = BeginningBlockDet, EndBlockDet
                ! if logged, accumulate the number of spawn events
                if (tLogNumSpawns) then
                    cum_det(IlutBits%ind_spawn) = cum_det(IlutBits%ind_spawn) + &
                                                  SpawnedParts(IlutBits%ind_spawn, i)
                end if
                ! Annihilate in this block seperately for walkers of different types.
                do part_type = 1, lenof_sign
                    if (tHistSpawn) then
                        call extract_sign(SpawnedParts(:, i), SpawnedSign)
                        call extract_sign(SpawnedParts(:, BeginningBlockDet), &
                        call HistAnnihilEvent(SpawnedParts, SpawnedSign, &
                                              temp_sign, part_type)
                    end if

                    call FindResidualParticle(cum_det, SpawnedParts(:, i), &
                                              part_type, iter_data, VecInd, Parent_Array_Ind)

                end do
            end do ! Loop over particle type.

            ! Copy details into the final array.
            call extract_sign(cum_det, temp_sign)

            if ((sum(abs(temp_sign)) > 1.e-12_dp) .or. (tFillingStochRDMonFly .and. (.not. tNoNewRDMContrib))) then
                ! Transfer all info into the other array.
                ! Usually this is only done if the final sign on the compressed
                ! Dj is not equal to zero. But in the case of the stochastic RDM,
                ! we are concerned with the sign of Dj in the CurrentDets array,
                ! not the newly spawned sign.  We still want to check if Dj has
                ! a non-zero Cj in CurrentDets, so we need to carry this Dj
                ! through to the stage of checking CurrentDets regardless of
                ! the sign here.  Also getting rid of them here would make the
                ! biased sign of Ci slightly wrong.

                ! for the stochastic GUGA RDMs we lose the information about
                ! the RDMs here in the SpawnedParts array
                ! So we have to use the Spawned_Parents array..
                SpawnedParts2(0:NIfTot, VecInd) = cum_det(0:NIfTot)

                if (tPreCond .or. tReplicaEstimates) then
                    SpawnedParts2(IlutBits%ind_hdiag, VecInd) = cum_det(IlutBits%ind_hdiag)
                end if

                VecInd = VecInd + 1
                DetsMerged = DetsMerged + EndBlockDet - BeginningBlockDet

                ! Spawned_Parts_Zero is the number of spawned parts that are
                ! zero after compression of the spawned_parts list - and should
                ! have been removed from SpawnedParts if we weren't calculating
                ! the RDM - need this for a check later.
                if (IsUnoccDet(temp_sign)) Spawned_Parts_Zero = Spawned_Parts_Zero + 1
                ! All particles from block have been annihilated.
                DetsMerged = DetsMerged + EndBlockDet - BeginningBlockDet + 1
                ! This spawned entry will be removed - don't want to store any
                ! parents. Reset Parent_Array_Ind so that the parents will be
                ! written over.

                ! For semi-stochastic simulations, store this state so that we
                ! can set the flags of the next state to be the same, if it is
                ! the same state.
                temp_det(0:NIfTot) = cum_det(0:NIfTot)
            end if

            ! Move onto the next block of determinants.
            BeginningBlockDet = CurrentBlockDet

        end do

        if (tFillingStochRDMonFly .and. (.not. tNoNewRDMContrib)) No_Spawned_Parents = Parent_Array_Ind - 1
        ! This is the new number of unique spawned determinants on the processor.
        ValidSpawned = ValidSpawned - DetsMerged
        if (ValidSpawned /= (VecInd - 1)) then
            call stop_all(t_r, "Error in compression of spawned list")
        end if

        ! Want the compressed list in spawnedparts at the end of it - swap
        ! pointers around.
        PointTemp => SpawnedParts2
        SpawnedParts2 => SpawnedParts
        SpawnedParts => PointTemp

        if (tAutoAdaptiveShift) then
            PointTemp => SpawnInfo2
            SpawnInfo2 => SpawnInfo
            SpawnInfo => PointTemp
        end if
    end subroutine CompressSpawnedList