attempt_die_normal Function

public function attempt_die_normal(DetCurr, Kii, RealwSign, WalkExcitLevel, DetPosition) result(ndie)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: DetCurr(nel)
real(kind=dp), intent(in) :: Kii
real(kind=dp), intent(in), dimension(lenof_sign) :: RealwSign
integer, intent(in) :: WalkExcitLevel
integer, intent(in), optional :: DetPosition

Return Value real(kind=dp), dimension(lenof_sign)


Contents

Source Code


Source Code

    function attempt_die_normal(DetCurr, Kii, realwSign, WalkExcitLevel, DetPosition) result(ndie)

        ! Should we kill the particle at determinant DetCurr.
        ! The function allows multiple births (if +ve shift), or deaths from
        ! the same particle. The returned number is the number of deaths if
        ! positive, and the
        !
        ! In:  DetCurr - The determinant to consider
        !      Kii     - The diagonal matrix element of DetCurr (-Ecore)
        !      wSign   - The sign of the determinant being considered. If
        !                |wSign| > 1, attempt to die multiple particles at
        !                once (multiply probability of death by |wSign|)
        ! Ret: ndie    - The number of deaths (if +ve), or births (If -ve).

        integer, intent(in) :: DetCurr(nel)
        real(dp), dimension(lenof_sign), intent(in) :: RealwSign
        real(dp), intent(in) :: Kii
        real(dp), dimension(lenof_sign) :: ndie
        integer, intent(in) :: WalkExcitLevel
        integer, intent(in), optional :: DetPosition
        character(*), parameter :: t_r = 'attempt_die_normal'

        integer(kind=n_int) :: iLutnI(0:niftot)
        integer :: N_double

        real(dp) :: probsign, r
        real(dp), dimension(inum_runs) :: fac
        integer :: i, run
#ifdef CMPLX_
        real(dp) :: rat(2)
#else
        real(dp) :: rat(1)
#endif
        real(dp) :: shift
        real(dp) :: relShiftOffset ! ShiftOffset relative to Hii

        unused_var(DetCurr)

        do i = 1, inum_runs
            if (t_cepa_shift) then

                fac(i) = tau * (Kii - (DiagSft(i) - cepa_shift(i, WalkExcitLevel)))
                call log_death_magnitude(Kii - (DiagSft(i) - cepa_shift(i, WalkExcitLevel)))

            else if (tSkipRef(i) .and. all(DetCurr == projEdet(:, i))) then
                !If we are fixing the population of reference det, skip death/birth
                fac(i) = 0.0
            else
                ! rescale the shift. It is better to wait till the shift starts varying
                if (tAdaptiveShift .and. .not. tSinglePartPhase(i)) then
                    if (tAS_Offset) then
                        !Since DiagSft is relative to Hii, ShiftOffset should also be adjusted
                        relShiftOffset = ShiftOffset(i) - Hii
                    else
                        !Each replica's shift should be scaled using its own reference
                        relShiftOffset = proje_ref_energy_offsets(i) ! i.e. E(Ref(i)) - Hii
                    end if
                  shift = relShiftOffset + (DiagSft(i) - relShiftOffset) * shiftFactorFunction(DetPosition, i, mag_of_run(RealwSign, i))
                else
                    shift = DiagSft(i)
                end if

                fac(i) = tau * (Kii - shift)

                if (t_precond_hub .and. Kii > 10.0_dp) then
                    call log_death_magnitude(((Kii - shift) / (1.0_dp + Kii)))
                else
                    call log_death_magnitude(Kii - shift)
                end if
            end if

            if (t_precond_hub .and. Kii > 10.0_dp) then
                fac(i) = fac(i) / (1.0_dp + Kii)
            end if
        end do

        if (any(fac > 1.0_dp)) then
            if (any(fac > 2.0_dp)) then
                if ((tau_search_method /= possible_tau_search_methods%OFF) .or. t_scale_tau_to_death) then
                    ! If we are early in the calculation, and are using tau
                    ! searching, then this is not a big deal. Just let the
                    ! searching deal with it
                    write(stderr, '("** WARNING ** Death probability > 2: Algorithm unstable.")')
                    write(stderr, '("** WARNING ** Truncating spawn to ensure stability")')
                    do i = 1, inum_runs
                        fac(i) = min(2.0_dp, fac(i))
                    end do
                else
                    print *, "Diag sft", DiagSft
                    print *, "Death probability", fac
                    call stop_all(t_r, "Death probability > 2: Algorithm unstable. Reduce timestep.")
                end if
            else
                write (stderr, '("** WARNING ** Death probability > 1: Creating Antiparticles. "&
                    & //"Timestep errors possible: ")', advance='no')
                do i = 1, inum_runs
                    write(stderr, '(1X,f13.7)', advance='no') fac(i)
                end do
                write(stderr, '()')
            end if
        end if

        if ((tRealCoeffByExcitLevel .and. (WalkExcitLevel <= RealCoeffExcitThresh)) &
            .or. tAllRealCoeff) then
            do run = 1, inum_runs
                ndie(min_part_type(run)) = fac(run) * abs(realwSign(min_part_type(run)))
#ifdef CMPLX_
                ndie(max_part_type(run)) = fac(run) * abs(realwSign(max_part_type(run)))
#endif
            end do
        else
            do run = 1, inum_runs

                ! Subtract the current value of the shift, and multiply by tau.
                ! If there are multiple particles, scale the probability.

                rat(:) = fac(run) * abs(realwSign(min_part_type(run):max_part_type(run)))

                ndie(min_part_type(run):max_part_type(run)) = real(int(rat), dp)
                rat(:) = rat(:) - ndie(min_part_type(run):max_part_type(run))

                ! Choose to die or not stochastically
                r = genrand_real2_dSFMT()
                if (abs(rat(1)) > r) ndie(min_part_type(run)) = &
                    ndie(min_part_type(run)) + real(nint(sign(1.0_dp, rat(1))), dp)
#ifdef CMPLX_
                r = genrand_real2_dSFMT()
                if (abs(rat(2)) > r) ndie(max_part_type(run)) = &
                    ndie(max_part_type(run)) + real(nint(sign(1.0_dp, rat(2))), dp)
#endif
            end do
        end if
    end function attempt_die_normal