subroutine calcFullStartL2R_stochastic(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 = "calcFullStartL2R_stochastic"
integer :: i, st, en, se, gen, j, k, l, elecInd, holeInd
type(WeightObj_t) :: weights
real(dp) :: branch_pgen, temp_pgen, p_orig, orb_pgen, rdm_mat
HElement_t(dp) :: integral
procedure(calc_pgen_general), pointer :: calc_pgen_yix_start
! the integral contributions to this mixed excitaiton full starts are
! more involved then i have imagined. since all the orbitals from below
! the first stepvector change(which has to happen at some point of the
! excitation, or else it is single-like,(which is also still unclear
! how i should implement that)) contribute to the same excitation.
! similar to the full-start full-stop mixed generator case.
! to ensure a stepvector change i could write some kind of restricted
! double excitiation update function, which accumulates switch
! probability with decreasing switch possibilities...
! and then during the excitation creation i have to save the first
! stepvector change and depending on that calculate the remaining
! integral contribution. todo
! and same with the fullstartr2l, and fullstoprl2/l2r
! for both full starts and full stop routines, the contributions are
! the same atleast.
ASSERT(isProperCSF_ilut(ilut))
ASSERT(.not. isZero(ilut, excitInfo%fullStart))
ASSERT(.not. isThree(ilut, excitInfo%fullStart))
! create the fullStart
st = excitInfo%fullStart
en = excitInfo%fullEnd
se = excitInfo%firstEnd
gen = excitInfo%lastGen
! create correct weights:
if (present(opt_weight)) then
weights = opt_weight
else
weights = init_fullStartWeight(csf_i, se, en, negSwitches(se), &
posSwitches(se), csf_i%B_real(se))
end if
! in the case of the approximate exchange excitations I need to
! force a switch at the beginning
if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
! do the switch
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
! in case of early exit pgen should be set to 0
pgen = 0.0_dp
check_abort_excit(branch_pgen, t)
! in the mixed fullstart case there has to be a switch at some point
! in the double overlap region, or else it would be a single-like
! excitation -> so check if x1 matrix element is 0
! then for the overlap region i need a double update routine, which
! somehow garuantees a switch happens at some point, to avoid
! single excitations
do i = st + 1, se - 1
call doubleUpdateStochastic(ilut, csf_i, i, excitInfo, &
weights, negSwitches, posSwitches, t, branch_pgen)
! to keep it general, i cant only check weights in doubleUpdate
if (near_zero(extract_matrix_element(t, 2)) .or. near_zero(branch_pgen)) then
t = 0_n_int
return
end if
end do
! then deal with specific semi-stop
! and update weights here
weights = weights%ptr
call calcLoweringSemiStopStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
posSwitches, t, branch_pgen)
! check validity
check_abort_excit(branch_pgen, t)
excitInfo%currentGen = excitInfo%lastGen
do i = se + 1, en - 1
call singleStochasticUpdate(ilut, csf_i, i, excitInfo, weights, posSwitches, &
negSwitches, t, temp_pgen)
branch_pgen = branch_pgen * temp_pgen
if (t_trunc_guga_pgen .or. &
(t_trunc_guga_pgen_noninits .and. .not. is_init_guga)) then
if (branch_pgen < trunc_guga_pgen) then
t = 0_n_int
return
end if
end if
! check validity
check_abort_excit(branch_pgen, t)
end do
call singleStochasticEnd(csf_i, excitInfo, t)
! 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
elecInd = en
holeInd = se
rdm_mat = extract_matrix_element(t, 2)
call calc_orbital_pgen_contrib_start(&
csf_i, [2 * st, 2 * elecInd], holeInd, orb_pgen)
p_orig = orb_pgen * branch_pgen / real(ElecPairs, dp)
if (csf_i%stepvector(elecInd) == 3) p_orig = p_orig * 2.0_dp
end if
! put everything in first entry
call encode_matrix_element(t, extract_matrix_element(t, 1) + &
extract_matrix_element(t, 2), 1)
call encode_matrix_element(t, 0.0_dp, 2)
! now have to check if there has be a change in the double overlap
! region, or else its just a single-like excitation:
! todo write a routine for that
! switchFlag = checkForSwitch(ilut, t, st, se-1)
! if a switch happened i have to calculated the additional contributing
! matrix elements from all indices below the fullstart
! can i formulate that in terms of the already calulated matrix
! element, as in the diagonal matrix element calculation?
! probably... but todo!
! update: since i can formulate everything in terms of the already
! calculated matrix element i can abort here if it is zero
if (near_zero(extract_matrix_element(t, 1))) then
t = 0_n_int
return
end if
global_excitInfo = excitInfo
if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
call calc_mixed_start_contr_approx(t, csf_i, excitInfo, integral)
pgen = branch_pgen
else
call calc_mixed_start_contr_integral(ilut, csf_i, t, excitInfo, &
integral)
call calc_mixed_start_contr_pgen(ilut, csf_i, t, excitInfo, &
branch_pgen, pgen)
! call calc_mixed_start_contr_sym(ilut, csf_i, t, excitInfo, branch_pgen, pgen, integral)
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=excitInfo%typ), x0=0.0_dp, &
x1=rdm_mat * pgen / p_orig)
end if
end if
! and finally update the matrix element with all contributions
call update_matrix_element(t, (get_umat_el(st, se, en, st) + &
get_umat_el(se, st, st, en)) / 2.0_dp + integral, 1)
end subroutine calcFullStartL2R_stochastic