subroutine calcFullStartFullStopMixedStochastic(ilut, csf_i, excitInfo, t, pgen, &
posSwitches, negSwitches, opt_weight)
integer(n_int), intent(in) :: ilut(0:nifguga)
type(CSF_Info_t), intent(in) :: csf_i
type(ExcitationInformation_t), intent(inout) :: excitInfo
integer(n_int), intent(out) :: t(0:nifguga)
real(dp), intent(out) :: pgen
real(dp), intent(in) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
type(WeightObj_t), intent(in), optional :: opt_weight
character(*), parameter :: this_routine = "calcFullStartFullStopMixedStochastic"
type(WeightObj_t) :: weights
real(dp) :: branch_pgen, temp_pgen, above_cpt, below_cpt, rdm_mat, p_orig
HElement_t(dp) :: integral
integer :: iOrb, i, j, k, l, typ
ASSERT(.not. isZero(ilut, excitInfo%fullStart))
ASSERT(.not. isThree(ilut, excitInfo%fullStart))
ASSERT(.not. isZero(ilut, excitInfo%fullEnd))
ASSERT(.not. isThree(ilut, excitInfo%fullEnd))
if (present(opt_weight)) then
weights = opt_weight
else
if (t_approx_exchange .or. (t_approx_exchange_noninits .and. &
(.not. is_init_guga))) then
weights = init_forced_end_exchange_weight(csf_i, excitInfo%fullEnd)
else
weights = init_doubleWeight(csf_i, excitInfo%fullEnd)
end if
end if
if (t_approx_exchange .or. (t_approx_exchange_noninits .and. &
(.not. is_init_guga))) then
call forced_mixed_start(ilut, csf_i, excitInfo, t, branch_pgen)
else
call mixedFullStartStochastic(ilut, csf_i, excitInfo, weights, posSwitches, &
negSwitches, t, branch_pgen)
end if
! set pgen to 0 in case of early exit
pgen = 0.0_dp
! in this excitation there has to be a switch somewhere in the double
! overlap region, so only x1 matrix element counts essentially
! and this can be 0 at the full-start already -> so check that here
! and in case abort
! should i do that in the mixedFullStartStochastic routine
! do that x1 matrix element in the routine and only check probWeight here
check_abort_excit(branch_pgen, t)
temp_pgen = 1.0_dp
do iOrb = excitInfo%fullStart + 1, excitInfo%fullEnd - 1
call doubleUpdateStochastic(ilut, csf_i, iOrb, excitInfo, weights, negSwitches, &
posSwitches, t, branch_pgen)
! zero x1 - elements can also happen in the double update
if (near_zero(extract_matrix_element(t, 2)) .or. near_zero(branch_pgen)) then
t = 0_n_int
return
end if
end do
call mixedFullStopStochastic(ilut, csf_i, excitInfo, t)
! check if there was a change in the stepvector in the double
! overlap region
if (.not. near_zero(extract_matrix_element(t, 1))) then
t = 0_n_int
return
end if
! if we do RDMs also store the x0 and x1 coupling coeffs
! and I need to do it before the routines below since excitInfo
! gets changed there
if (tFillingStochRDMOnFly) then
! i need to unbias against the total pgen later on in the
! RDM sampling otherwise the rdm-bias factor is not correct!
! encode the necessary information in the rdm-matele!
i = excitInfo%i
j = excitInfo%j
k = excitInfo%k
l = excitInfo%l
typ = excitInfo%typ
rdm_mat = extract_matrix_element(t, 2)
call calc_orbital_pgen_contr(csf_i, [2 * i, 2 * j], above_cpt, &
below_cpt)
p_orig = (above_cpt + below_cpt) * branch_pgen
if (.not. (t_heisenberg_model .or. t_tJ_model)) then
p_orig = p_orig / real(ElecPairs, dp)
end if
end if
global_excitInfo = excitInfo
if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
if (getDeltaB(t) == 0) then
t = 0_n_int
return
end if
call calc_mixed_contr_integral(ilut, csf_i, t, excitInfo%fullStart, &
excitInfo%fullEnd, integral)
pgen = branch_pgen
! just to be save that a switch always happens at the end
! print that out for now
else
call calc_mixed_contr_integral(ilut, csf_i, t, excitInfo%fullStart, &
excitInfo%fullEnd, integral)
if (.not. near_zero(integral)) then
call calc_mixed_contr_pgen(ilut, csf_i, t, excitInfo, pgen)
end if
! call calc_mixed_contr_sym(ilut, csf_i, t, excitInfo, pgen, integral)
end if
if (near_zero(integral)) then
t = 0_n_int
pgen = 0.0_dp
return
end if
if (tFillingStochRDMOnFly) then
if (.not. near_zero(p_orig)) then
call encode_stochastic_rdm_info(GugaBits, t, rdm_ind= &
contract_2_rdm_ind(i, j, k, l, excit_lvl = 2, &
excit_typ=typ), x0 = 0.0_dp, x1 = rdm_mat * pgen / p_orig)
end if
end if
call encode_matrix_element(t, 0.0_dp, 2)
call encode_matrix_element(t, integral, 1)
end subroutine calcFullStartFullStopMixedStochastic