function attempt_create_realtime(DetCurr, iLutCurr, RealwSign, nJ, iLutnJ, &
prob, HElGen, ic, ex, tParity, walkExcitLevel, part_type, AvSignCurr, &
AvExPerWalker, RDMBiasFacCurr, precond_fac) result(child)
! create a specific attempt_create function for the real-time fciqmc
! to avoid preprocessor flag jungle..
integer, intent(in) :: DetCurr(nel), nJ(nel)
integer, intent(in) :: part_type ! 1 = Real parent particle, 2 = Imag parent particle
integer(kind=n_int), intent(in) :: iLutCurr(0:NIfTot)
integer(kind=n_int), intent(inout) :: iLutnJ(0:niftot)
integer, intent(in) :: ic, ex(2, ic), walkExcitLevel
real(dp), dimension(lenof_sign), intent(in) :: RealwSign
logical, intent(in) :: tParity
real(dp), intent(inout) :: prob
real(dp), dimension(lenof_sign) :: child
real(dp), dimension(lenof_sign), intent(in) :: AvSignCurr
real(dp), intent(in) :: AvExPerWalker
real(dp), intent(out) :: RDMBiasFacCurr
real(dp), intent(in) :: precond_fac
HElement_t(dp), intent(inout) :: HElGen
character(*), parameter :: this_routine = 'attempt_create_realtime'
real(dp) :: walkerweight, pSpawn, nSpawn, MatEl, p_spawn_rdmfac, &
sepSign, fac_unused
integer :: extracreate, tgt_cpt, component, iUnused
integer :: TargetExcitLevel, tmp_ex(2, ic)
logical :: tRealSpawning
real(dp) :: rh_imag
HElement_t(dp) :: rh_used
unused_var(precond_fac)
unused_var(AvSignCurr)
! This is crucial
child = 0.0_dp
! If each walker does not have exactly one spawning attempt
! (if AvMCExcits /= 1.0_dp) then the probability of an excitation
! having been chosen, prob, must be altered accordingly.
prob = prob * AvExPerWalker
! In the case of using HPHF, and when tGenMatHEl is on, the matrix
! element is calculated at the time of the excitation generation,
! and returned in HElGen. In this case, get_spawn_helement simply
! returns HElGen, rather than recomputing the matrix element.
tmp_ex(1, :) = ex(2, :)
tmp_ex(2, :) = ex(1, :)
rh_used = get_spawn_helement(nJ, DetCurr, iLutnJ, ilutCurr, ic, tmp_ex, &
tParity, HElGen)
tRealSpawning = .false.
if (tAllRealCoeff) then
tRealSpawning = .true.
else if (tRealCoeffByExcitLevel) then
TargetExcitLevel = FindBitExcitLevel(iLutRef, iLutnJ)
if (TargetExcitLevel <= RealCoeffExcitThresh) &
tRealSpawning = .true.
end if
! do the whole real-time shabang here, since it also depends on the
! fact if it is a pure real-hamiltonian(or atleast we can save some effort)
! change in code here for the real-time fciqmc
! for now itsonly implemented for pure real hamiltonians
! but the -i * H in the diff. eq. makes it kind of all imaginary
! for the walker dynamics(except no influence on the conjugation
! of H for H_ij -> H_ji*
! todo: the case of a paritally imaginary Hamiltonian as input
! has also to be considered later on!
! here i have then contributions from both real and imaginary
! walker populations: when writing the Hamiltonian as: H = H + iJ
! and a determinant weight as: n + ic:
! n_j + ic_j = -i(H_ij + iJ_ij)^+ (n_i + ic_i)
! n_j + ic_j = -i(H_ji - iJ_ji) (n_i + ic_i) =>
! leads to =>
! n_j = H_ji c_i - J_ji n_i
! c_j = -H_ji n_i - J_ji c_i
! ... but this should also be with the "normal" complex Hamiltonians
! or? talk to ali..
! yes this is exactly what is done already.. now implement that
! also for the additional -i in the real-time walker dynamics
if (.not. t_complex_ints .and. .not. t_rotated_time) then
! if it is a pure real-hamiltonian there is only spawing from
! real to complex walkers and v.v.
tgt_cpt = rotate_part(part_type)
walkerweight = sign(1.0_dp, RealwSign(part_type))
MatEl = real(rh_used, dp)
! spawn from real-parent to imaginary child: no sign change
! from imaginary to real -> sign change
if (mod(part_type, 2) == 0) walkerweight = -walkerweight
nSpawn = -tau * MatEl * walkerweight / prob
! n.b. if we ever end up with |walkerweight| /= 1, then this
! will need to ffed further through.
if ((tau_search_method == possible_tau_search_methods%CONVENTIONAL) .and. (.not. tFillingStochRDMonFly)) then
call log_spawn_magnitude(ic, ex, matel, prob)
end if
! Keep track of the biggest spawn this cycle
max_cyc_spawn = max(abs(nSpawn), max_cyc_spawn)
if (tRealSpawning) then
! Continuous spawning. Add in acceptance probabilities.
if (tRealSpawnCutoff .and. &
abs(nSpawn) < RealSpawnCutoff) then
p_spawn_rdmfac = abs(nSpawn) / RealSpawnCutoff
nSpawn = RealSpawnCutoff &
* stochastic_round(nSpawn / RealSpawnCutoff)
else
p_spawn_rdmfac = 1.0_dp !The acceptance probability of some kind of child was equal to 1
end if
else
if (abs(nSpawn) >= 1) then
p_spawn_rdmfac = 1.0_dp !We were certain to create a child here.
! This is the special case whereby if P_spawn(j | i) > 1,
! then we will definitely spawn from i->j.
! I.e. the pair Di,Dj will definitely be in the SpawnedParts list.
! We don't care about multiple spawns - if it's in the list, an RDM contribution will result
! regardless of the number spawned - so if P_spawn(j | i) > 1, we treat it as = 1.
else
p_spawn_rdmfac = abs(nSpawn)
end if
! How many children should we spawn?
! And round this to an integer in the usual way
! HACK: To use the same number of random numbers for the tests.
nSpawn = real(stochastic_round(nSpawn), dp)
end if
! And create the parcticles
child(tgt_cpt) = nSpawn
else
! have to loop over the tgt_cpt similar to the complex impl
! if the Hamiltonian has real and imaginary components do it
! similarily to complex implementation with H <-> J switched
! rmneci_setup: adjusted for multirun, fixed complex -> real spawns
do component = 1, (lenof_sign / inum_runs)
tgt_cpt = min_part_type(part_type_to_run(part_type)) - 1 + component
! keep track of the sign due to the kind of spawn event
sepSign = 1.0_dp
! if (part_type == 2 .and. inum_runs == 1) component = 3 - tgt_cpt !?
walkerweight = sign(1.0_dp, RealwSign(part_type))
if (mod(part_type, 2) == 0 .and. component == 1) &
sepSign = (-1.0_dp)
if (t_real_time_fciqmc) then
#ifdef CMPLX_
rh_imag = real(aimag(rh_used), dp)
#else
rh_imag = 0.0_dp
#endif
! part_type is given as input, for that part_type, the real part of
! the HElement is used if rotation occurs and the imaginary part if not
if (mod(component, 2) == mod(part_type, 2)) then
! spawn part_type -> part_type
MatEl = -rh_imag * tau_real - real(rh_used, dp) * tau_imag
else
! spawn part_type -> rotate_part(part_type)
MatEl = real(rh_used, dp) * tau_real - rh_imag * tau_imag
end if
end if
nSpawn = -sepSign * MatEl * walkerweight / prob
! n.b. if we ever end up with |walkerweight| /= 1, then this
! will need to ffed further through.
if ((tau_search_method == possible_tau_search_methods%CONVENTIONAL) .and. (.not. tFillingStochRDMonFly)) then
call log_spawn_magnitude(ic, ex, matel, prob)
end if
! Keep track of the biggest spawn this cycle
max_cyc_spawn = max(abs(nSpawn), max_cyc_spawn)
if (tRealSpawning) then
! Continuous spawning. Add in acceptance probabilities.
if (tRealSpawnCutoff .and. &
abs(nSpawn) < RealSpawnCutoff) then
p_spawn_rdmfac = abs(nSpawn) / RealSpawnCutoff
nSpawn = RealSpawnCutoff &
* stochastic_round(nSpawn / RealSpawnCutoff)
else
p_spawn_rdmfac = 1.0_dp !The acceptance probability of some kind of child was equal to 1
end if
else
if (abs(nSpawn) >= 1) then
p_spawn_rdmfac = 1.0_dp !We were certain to create a child here.
! This is the special case whereby if P_spawn(j | i) > 1,
! then we will definitely spawn from i->j.
! I.e. the pair Di,Dj will definitely be in the SpawnedParts list.
! We don't care about multiple spawns - if it's in the list, an RDM contribution will result
! regardless of the number spawned - so if P_spawn(j | i) > 1, we treat it as = 1.
else
p_spawn_rdmfac = abs(nSpawn)
end if
! How many children should we spawn?
! And round this to an integer in the usual way
! HACK: To use the same number of random numbers for the tests.
nSpawn = real(stochastic_round(nSpawn), dp)
end if
! And create the parcticles (in the correct run)
child(tgt_cpt) = nSpawn
end do
end if
if (tFillingStochRDMonFly) then
if (.not. near_zero(child(part_type))) then
!Only add in contributions for spawning events within population 1
!(Otherwise it becomes tricky in annihilation as spawnedparents doesn't tell you which population
!the event came from at present)
call calc_rdmbiasfac(p_spawn_rdmfac, prob, realwSign(part_type), RDMBiasFacCurr)
else
RDMBiasFacCurr = 0.0_dp
end if
else
! Not filling the RDM stochastically, bias is zero.
RDMBiasFacCurr = 0.0_dp
end if
! Avoid compiler warnings
iUnused = walkExcitLevel
end function attempt_create_realtime