create_particle Subroutine

public subroutine create_particle(nJ, iLutJ, child, part_type, hdiag_spawn, err, ilutI, SignCurr, WalkerNo, RDMBiasFacCurr, WalkersToSpawn, matel, ParentPos)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nJ(nel)
integer(kind=n_int), intent(in) :: iLutJ(0:niftot)
real(kind=dp), intent(in) :: child(lenof_sign)
integer, intent(in) :: part_type
real(kind=dp), intent(in) :: hdiag_spawn
integer, intent(out) :: err
integer(kind=n_int), intent(in), optional :: ilutI(0:niftot)
real(kind=dp), intent(in), optional :: SignCurr(lenof_sign)
integer, intent(in), optional :: WalkerNo
real(kind=dp), intent(in), optional :: RDMBiasFacCurr
integer, intent(in), optional :: WalkersToSpawn
real(kind=dp), intent(in), optional :: matel
integer, intent(in), optional :: ParentPos

Contents

Source Code


Source Code

    subroutine create_particle(nJ, iLutJ, child, part_type, hdiag_spawn, err, ilutI, SignCurr, &
                               WalkerNo, RDMBiasFacCurr, WalkersToSpawn, matel, ParentPos)
        ! Create a child in the spawned particles arrays. We spawn particles
        ! into a separate array, but non-contiguously. The processor that the
        ! newly-spawned particle is going to be sent to has to be determined,
        ! and then it will be put into the appropriate element determined by
        ! ValidSpawnedList

        ! 'type' of the particle - i.e. real/imag
        integer, intent(in) :: nJ(nel), part_type
        integer(n_int), intent(in) :: iLutJ(0:niftot)
        real(dp), intent(in) :: child(lenof_sign)
        integer, intent(out) :: err
        HElement_t(dp), intent(in) :: hdiag_spawn
        integer(n_int), intent(in), optional :: ilutI(0:niftot)
        real(dp), intent(in), optional :: SignCurr(lenof_sign)
        integer, intent(in), optional :: WalkerNo
        real(dp), intent(in), optional :: RDMBiasFacCurr
        integer, intent(in), optional :: WalkersToSpawn
        real(dp), intent(in), optional :: matel
        integer, intent(in), optional :: ParentPos
        integer :: proc, j, run
        real(dp) :: r
        integer, parameter :: flags = 0
        logical :: list_full, allowed_child
        character(*), parameter :: this_routine = 'create_particle'

        logical :: parent_init
        real(dp)  :: weight_acc, weight_rej, weight_rev, weight_den, weight_den2

        ! error flag to indicate if the attempt was successful
        err = 0
        !Ensure no cross spawning between runs - run of child same as run of
        !parent
        run = part_type_to_run(part_type)
        ASSERT(sum(abs(child)) - sum(abs(child(min_part_type(run):max_part_type(run)))) < 1.0e-12_dp)
        ! Determine which processor the particle should end up on in the
        ! DirectAnnihilation algorithm.
        proc = DetermineDetNode(nel, nJ, 0)    ! (0 -> nNodes-1)

        if (checkValidSpawnedList(proc)) then
            err = 1
            return
        end if
        !We initially encode no flags
        call encode_bit_rep(SpawnedParts(:, ValidSpawnedList(proc)), iLutJ, &
                            child, flags)

        ! If the parent was an initiator then set the initiator flag for the
        ! child, to allow it to survive.
        if (tTruncInitiator) then
            allowed_child = .false.
            ! optionally: allow all spawns with a given sign
            if (allowedSpawnSign /= 0) then
                if (allowedSpawnSign * child(part_type) * SignCurr(part_type) > 0) &
                    allowed_child = .true.
            end if
            if (allowed_child .or. test_flag(ilutI, get_initiator_flag(part_type))) then
                call set_flag(SpawnedParts(:, ValidSpawnedList(proc)), get_initiator_flag(part_type))
            end if
        end if

        if (tAutoAdaptiveShift) then
            SpawnInfo(SpawnParentIdx, ValidSpawnedList(proc)) = ParentPos
            SpawnInfo(SpawnRun, ValidSpawnedList(proc)) = run
            if (tAAS_MatEle) then
                weight_acc = abs(matel)
                weight_rej = abs(matel)
            else if (tAAS_MatEle2) then
                weight_den = abs((get_diagonal_matel(nJ, ilutJ) - Hii) - DiagSft(run))
                if (weight_den < AAS_DenCut) then
                    weight_den = AAS_DenCut
                end if
                weight_acc = abs(matel) / weight_den
                weight_rej = abs(matel) / weight_den
            else if (tAAS_MatEle3) then
                weight_den = abs((get_diagonal_matel(nJ, ilutJ) - Hii) - DiagSft(run))
                if (weight_den < AAS_DenCut) then
                    weight_den = AAS_DenCut
                end if
                weight_acc = 1.0_dp
                weight_rej = abs(matel) / weight_den
            else if (tAAS_MatEle4) then
                weight_den = abs((get_diagonal_matel(nJ, ilutJ) - Hii) - DiagSft(run))
                if (weight_den < AAS_DenCut) then
                    weight_den = AAS_DenCut
                end if
                weight_den2 = abs((get_diagonal_matel(nJ, ilutJ) - Hii))
                if (weight_den2 < AAS_DenCut) then
                    weight_den2 = AAS_DenCut
                end if
                weight_rej = abs(matel) / weight_den
                weight_acc = abs(matel) / weight_den2
            else
                weight_acc = 1.0_dp
                weight_rej = 1.0_dp
            end if

            !Enocde weight, which is real, as an integer
            SpawnInfo(SpawnWeightAcc, ValidSpawnedList(proc)) = &
                transfer(weight_acc, SpawnInfo(SpawnWeightAcc, ValidSpawnedList(proc)))
            SpawnInfo(SpawnWeightRej, ValidSpawnedList(proc)) = &
                transfer(weight_rej, SpawnInfo(SpawnWeightRej, ValidSpawnedList(proc)))

        end if

        ! store global data - number of spawns
        ! Is using the pure initiator space option, then if this spawning
        ! occurs to within the defined initiator space (regardless of
        ! whether or not it is occupied already), then it should never be
        ! rejected by the initiator criterion, so set the initiator
        ! flag now if not done already.
        if (tPureInitiatorSpace) then
            if (.not. test_flag(SpawnedParts(:, ValidSpawnedList(proc)), &
                                get_initiator_flag(part_type))) then
                if (is_in_initiator_space(SpawnedParts(:, ValidSpawnedList(proc)), nJ)) then
                    call set_flag(SpawnedParts(:, ValidSpawnedList(proc)), &
                                  get_initiator_flag(part_type))
                end if
            end if
        end if
        if (tLogNumSpawns) call log_spawn(SpawnedParts(:, ValidSpawnedList(proc)))

        if (tFillingStochRDMonFly) then
            ! We are spawning from ilutI to
            ! SpawnedParts(:,ValidSpawnedList(proc)). We want to store the
            ! parent (D_i) with the spawned child (D_j) so that we can add in
            ! Ci.Cj to the RDM later.
            ! The parent is nifd integers long, and stored in the second
            ! part of the SpawnedParts array from NIfTot+1 --> NIfTot+1+nifd
            call store_parent_with_spawned(RDMBiasFacCurr, WalkerNo, &
                                           ilutI, WalkersToSpawn, ilutJ, proc)

            ! in the GUGA case I also need to store the rdm info in the
            ! spawned parts arrays here
            if (tGUGA) then
                call transfer_stochastic_rdm_info(ilutJ, &
                                                  SpawnedParts(:, ValidSpawnedList(proc)), &
                                                  BitIndex_from=IlutBits, BitIndex_to=IlutBits)
            end if
        end if

        if (tPreCond .or. tReplicaEstimates) then
            call encode_spawn_hdiag(SpawnedParts(:, ValidSpawnedList(proc)), hdiag_spawn)
        end if

        ValidSpawnedList(proc) = ValidSpawnedList(proc) + 1

        ! Sum the number of created children to use in acceptance ratio.
        ! Note that if child is an array, it should only have one non-zero
        ! element which has changed.
        ! rmneci_setup: clarified dependence of run on part_type
        run = part_type_to_run(part_type)
        acceptances(run) = acceptances(run) + &
                           sum(abs(child(min_part_type(run):max_part_type(run))))

    end subroutine create_particle