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