subroutine calcDoubleR2L_stochastic(ilut, csf_i, excitInfo, t, branch_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) :: branch_pgen
real(dp), intent(in) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
type(WeightObj_t), intent(in), optional :: opt_weight
character(*), parameter :: this_routine = "calcDoubleR2L_stochastic"
integer :: iOrb, start2, ende1, ende2, start1, switch
type(WeightObj_t) :: weights
real(dp) :: temp_pgen
HElement_t(dp) :: integral
ASSERT(.not. isThree(ilut, excitInfo%fullStart))
ASSERT(.not. isZero(ilut, excitInfo%secondStart))
ASSERT(.not. isZero(ilut, excitInfo%firstEnd))
ASSERT(.not. isThree(ilut, excitInfo%fullEnd))
start1 = excitInfo%fullStart
start2 = excitInfo%secondStart
ende1 = excitInfo%firstEnd
ende2 = excitInfo%fullEnd
! : create correct weights:
if (present(opt_weight)) then
weights = opt_weight
else
weights = init_fullDoubleWeight(csf_i, start2, ende1, ende2, negSwitches(start2), &
negSwitches(ende1), posSwitches(start2), posSwitches(ende1), &
csf_i%B_real(start2), csf_i%B_real(ende1))
end if
call createStochasticStart_single(ilut, csf_i, excitInfo, weights, posSwitches, &
negSwitches, t, branch_pgen)
! check validity
check_abort_excit(branch_pgen, t)
do iOrb = start1 + 1, start2 - 1
call singleStochasticUpdate(ilut, csf_i, iOrb, excitInfo, weights, posSwitches, &
negSwitches, t, temp_pgen)
branch_pgen = branch_pgen * temp_pgen
! check validity
check_abort_excit(branch_pgen, t)
end do
! change weights... maybe need both single and double type weights
! then do lowering semi start
weights = weights%ptr
call calcLoweringSemiStartStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
posSwitches, t, branch_pgen)
! check validity
check_abort_excit(branch_pgen, t)
do iOrb = start2 + 1, ende1 - 1
call doubleUpdateStochastic(ilut, csf_i, iOrb, excitInfo, weights, negSwitches, &
posSwitches, t, branch_pgen)
! check validity
check_abort_excit(branch_pgen, t)
end do
! then update weights and and to lowering semi-stop
weights = weights%ptr
call calcRaisingSemiStopStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
posSwitches, t, branch_pgen)
! check validity
check_abort_excit(branch_pgen, t)
excitInfo%currentGen = excitInfo%lastGen
do iOrb = ende1 + 1, ende2 - 1
call singleStochasticUpdate(ilut, csf_i, iOrb, excitInfo, weights, posSwitches, &
negSwitches, t, temp_pgen)
branch_pgen = branch_pgen * temp_pgen
! check validity
check_abort_excit(branch_pgen, t)
end do
! and finally to end step
call singleStochasticEnd(csf_i, excitInfo, t)
! if we do RDMs also store the x0 and x1 coupling coeffs
if (tFillingStochRDMOnFly) then
call encode_stochastic_rdm_info(GugaBits, t, rdm_ind= &
contract_2_rdm_ind(excitInfo%i, excitInfo%j, excitInfo%k, excitInfo%l, &
excit_lvl=2, excit_typ=excitInfo%typ), &
x0=extract_matrix_element(t, 1), &
x1=extract_matrix_element(t, 2))
end if
! for the additional contributing integrals:
! i have to consider that there might be a non-overlap excitation,
! which leads to the same excitation if there is no change in the
! stepvector in the overlap region!
! where i have to parts of the matrix elements, but which can be
! expressed in terms of the already calulated one i think..
! atleast in the case if R2L and L2R, as both the bottom and top
! parts are the same and only the semi-start ans semi-stop have to
! be modified..
! and also the x1 element has to be dismissed i guess..
! so i should change how i deal with the x1 elements and maybe
! keep them out here in these functions..
! determine if a switch happended:
switch = findFirstSwitch(ilut, t, start2 + 1, ende1)
! i have to rethink the matrix element contribution by the
! equivalent non-overlap excitations, since its possible that
! a +-2 excitations stays on the +-2 branch in the double overlap
! region all the time and then also can be produced by a non-overlap
! but then the matrix element is different then the calculation used
! right now..
! no!!! sinnce the excitations leading to such a csf wouldn be valid
! single excitations-.. since the deltaB values at the end would not
! match
if (switch > 0) then
! a switch happened and only mixed overlap contributes
! after determining how to deal with different x0 and x1 parts
! wait: if a switch happened i know that the x0-element is 0
! and only the x1-element of the overlap region counts!
integral = extract_matrix_element(t, 2) * (get_umat_el(start1, ende2, ende1, start2) + &
get_umat_el(ende2, start1, start2, ende1)) / 2.0_dp
if (near_zero(integral)) then
branch_pgen = 0.0_dp
t = 0
else
call encode_matrix_element(t, 0.0_dp, 2)
call encode_matrix_element(t, integral, 1)
end if
else
! no switch happened: so i have to think about the
! non-overlap contribution
! in the R2L and L2R case it is really easy if i keep the
! x0 and x1 elements seperate up until here.
! since the types of generators are the same and only the x0
! element(which does not get changed in the overlap region)
! are nececarry for the non-overlap excitation , it only gets
! a -t^2 = -1/2 factor...
! so to get the non-overlap version i need to multiply by -2.0!
integral = (-extract_matrix_element(t, 1) * (get_umat_el(start1, ende2, start2, ende1) + &
get_umat_el(ende2, start1, ende1, start2)) * 2.0_dp + (extract_matrix_element(t, 1) + &
extract_matrix_element(t, 2)) * (get_umat_el(start1, ende2, ende1, start2) + &
get_umat_el(ende2, start1, start2, ende1))) / 2.0_dp
if (near_zero(integral)) then
branch_pgen = 0.0_dp
t = 0_n_int
else
call encode_matrix_element(t, 0.0_dp, 2)
call encode_matrix_element(t, integral, 1)
end if
end if
end subroutine calcDoubleR2L_stochastic