subroutine calc_mixed_start_contr_pgen(ilut, csf_i, t, excitInfo, branch_pgen, &
pgen)
integer(n_int), intent(in) :: ilut(0:nifguga), t(0:nifguga)
type(CSF_Info_t), intent(in) :: csf_i
type(ExcitationInformation_t), value :: excitInfo
real(dp), value :: branch_pgen
real(dp), intent(out) :: pgen
character(*), parameter :: this_routine = "calc_mixed_start_contr_pgen"
integer :: st, se, en, elecInd, holeInd, sw, step, i
real(dp) :: orb_pgen, zero_weight, switch_weight, stay_weight
real(dp) :: bot_cont, posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
real(dp) :: new_pgen, start_mat, stay_mat, start_weight
type(WeightObj_t) :: weights
logical :: below_flag
st = excitInfo%fullStart
se = excitInfo%firstEnd
en = excitInfo%fullEnd
! depending on the type of excitaiton, calculation of orbital pgens
! change
if (excitInfo%typ == excit_type%fullStart_L_to_R) then
elecInd = en
holeInd = se
else if (excitInfo%typ == excit_type%fullstart_R_to_L) then
elecInd = se
holeInd = en
else
call stop_all(this_routine, "should not be here!")
end if
sw = findFirstSwitch(ilut, t, st, se)
! what can i precalculate beforehand?
step = csf_i%stepvector(st)
! do i actually deal with the actual start orbital influence??
! fuck i don't think so.. wtf..
call calc_orbital_pgen_contrib_start(csf_i, [2 * st, 2 * elecInd], &
holeInd, orb_pgen)
pgen = orb_pgen * branch_pgen
! since weights only depend on the number of switches at the
! semistop and semistop and full-end index i can calculate
! it beforehand for all?
excitInfo%fullStart = 1
excitInfo%secondStart = 1
call calcRemainingSwitches_excitInfo_double(csf_i, excitInfo, posSwitches, negSwitches)
weights = init_fullStartWeight(csf_i, se, en, negSwitches(se), &
posSwitches(se), csf_i%B_real(se))
! determine the original starting weight
zero_weight = weights%proc%zero(negSwitches(st), posSwitches(st), &
csf_i%B_real(st), weights%dat)
if (step == 1) then
switch_weight = weights%proc%plus(posSwitches(st), csf_i%B_real(st), &
weights%dat)
if (isOne(t, st)) then
bot_cont = Root2 * sqrt((csf_i%B_real(st) - 1.0_dp) / &
(csf_i%B_real(st) + 1.0_dp))
stay_weight = calcStayingProb(zero_weight, switch_weight, &
csf_i%B_real(st))
start_weight = zero_weight / (zero_weight + switch_weight)
else
bot_cont = -sqrt(2.0_dp / ((csf_i%B_real(st) - 1.0_dp) * &
(csf_i%B_real(st) + 1.0_dp)))
stay_weight = 1.0_dp - calcStayingProb(zero_weight, switch_weight, &
csf_i%B_real(st))
start_weight = switch_weight / (zero_weight + switch_weight)
end if
else
switch_weight = weights%proc%minus(negSwitches(st), &
csf_i%B_real(st), weights%dat)
if (isOne(t, st)) then
bot_cont = -sqrt(2.0_dp / ((csf_i%B_real(st) + 1.0_dp) * &
(csf_i%B_real(st) + 3.0_dp)))
stay_weight = 1.0_dp - calcStayingProb(zero_weight, switch_weight, &
csf_i%B_real(st))
start_weight = switch_weight / (zero_weight + switch_weight)
else
bot_cont = -Root2 * sqrt((csf_i%B_real(st) + 3.0_dp) / &
(csf_i%B_real(st) + 1.0_dp))
stay_weight = calcStayingProb(zero_weight, switch_weight, &
csf_i%B_real(st))
start_weight = zero_weight / (zero_weight + switch_weight)
end if
end if
ASSERT(.not. near_zero(start_weight))
! update the pgen stumbs here to reuse start_weight variable
new_pgen = stay_weight * branch_pgen / start_weight
! divide out the original starting weight:
branch_pgen = branch_pgen / start_weight
! loop from start backwards so i can abort at a d=1 & b=1 stepvalue
! also consider if bot_cont < EPS to avoid unnecarry calculations
if (.not. near_zero(bot_cont)) then
below_flag = .false.
do i = st - 1, 1, -1
if (csf_i%Occ_int(i) /= 1) cycle
! then check if thats the last stepvalue to consider
if (csf_i%stepvector(i) == 1 .and. csf_i%B_int(i) == 1) then
below_flag = .true.
end if
! then i need to calculate the orbital probability
! from the fact that this is a lowering into raising fullstart
! i know, more about the restrictions...
! and that fullend is the electron eg.
! depening on the type of excitation (r2l or l2r) the electron
! orbitals change here
call calc_orbital_pgen_contrib_start(csf_i, [2*i, 2*elecInd], &
holeInd, orb_pgen)
! then deal with the matrix element and branching probabilities
step = csf_i%stepvector(i)
! get both start and staying matrix elements -> and update
! matrix element contributions on the fly to avoid second loop!
call getDoubleMatrixElement(step, step, -1, gen_type%R, gen_type%L, &
csf_i%B_real(i), 1.0_dp, x1_element=start_mat)
call getDoubleMatrixElement(step, step, 0, gen_type%R, gen_type%L, &
csf_i%B_real(i), 1.0_dp, x1_element=stay_mat)
! another check.. although this should not happen
! except the other d = 1 & b = 1 condition is already met
! above, to not continue:
if (near_zero(stay_mat)) below_flag = .true.
zero_weight = weights%proc%zero(negSwitches(i), posSwitches(i), &
csf_i%B_real(i), weights%dat)
if (step == 1) then
switch_weight = weights%proc%plus(posSwitches(i), &
csf_i%B_real(i), weights%dat)
else
switch_weight = weights%proc%minus(negSwitches(i), &
csf_i%B_real(i), weights%dat)
end if
start_weight = zero_weight / (zero_weight + switch_weight)
stay_weight = calcStayingProb(zero_weight, switch_weight, &
csf_i%B_real(i))
! i think i could avoid the second loop over j
! if i express everything in terms of already calculated
! quantities!
! "normally" matrix element shouldnt be 0 anymore... still check
if (.not. near_zero(start_mat)) then
pgen = pgen + orb_pgen * start_weight * new_pgen
end if
if (below_flag) exit
if (t_trunc_guga_pgen .or. &
(t_trunc_guga_pgen_noninits .and. .not. is_init_guga)) then
if (new_pgen < trunc_guga_pgen) then
new_pgen = 0.0_dp
end if
end if
! update new_pgen for next cycle
new_pgen = stay_weight * new_pgen
end do
end if
! start to switch loop: here matrix elements are not 0!
! and its only db = 0 branch and no stepvalue change!
! if the start is the switch nothing happens
step = csf_i%stepvector(st)
! calculate the necarry values needed to formulate everything in terms
! of the already calculated quantities:
call getDoubleMatrixElement(step, step, -1, gen_type%L, gen_type%R, &
csf_i%B_real(st), 1.0_dp, x1_element=start_mat)
if (.not. near_zero(abs(start_mat))) then
do i = st + 1, sw - 1
! the good thing here is, i do not need to loop a second time,
! since i can recalc. the matrix elements and pgens on-the fly
! here the matrix elements should not be 0 or otherwise the
! excitation wouldnt have happended anyways
if (csf_i%Occ_int(i) /= 1) cycle
! calculate orbitals pgen first and cycle if 0
call calc_orbital_pgen_contrib_start(csf_i, [2 * i, 2 * elecInd], &
holeInd, orb_pgen)
step = csf_i%stepvector(i)
! update inverse product
#ifdef DEBUG_
call getDoubleMatrixElement(step, step, 0, gen_type%L, gen_type%R, &
csf_i%B_real(i), 1.0_dp, x1_element=stay_mat)
ASSERT(.not. near_zero(stay_mat))
#endif
! and update pgens also
zero_weight = weights%proc%zero(negSwitches(i), posSwitches(i), &
csf_i%B_real(i), weights%dat)
if (step == 1) then
switch_weight = weights%proc%plus(posSwitches(i), &
csf_i%B_real(i), weights%dat)
else
switch_weight = weights%proc%minus(negSwitches(i), &
csf_i%B_real(i), weights%dat)
end if
stay_weight = calcStayingProb(zero_weight, switch_weight, &
csf_i%B_real(i))
start_weight = zero_weight / (zero_weight + switch_weight)
branch_pgen = branch_pgen / stay_weight
if (t_trunc_guga_pgen .or. &
(t_trunc_guga_pgen_noninits .and. .not. is_init_guga)) then
if (branch_pgen * start_weight > trunc_guga_pgen) then
pgen = pgen + orb_pgen * branch_pgen * start_weight
end if
else
! and add up correctly
pgen = pgen + orb_pgen * branch_pgen * start_weight
end if
end do
! handle switch seperately (but only if switch > start)
if (sw > st) then
! check orb_pgen otherwise no influencce
call calc_orbital_pgen_contrib_start(csf_i, [2 * sw, 2 * elecInd], &
holeInd, orb_pgen)
if (.not. near_zero(orb_pgen)) then
step = csf_i%stepvector(sw)
zero_weight = weights%proc%zero(negSwitches(sw), posSwitches(sw), &
csf_i%B_real(sw), weights%dat)
! on the switch the original probability is:
if (step == 1) then
switch_weight = weights%proc%plus(posSwitches(sw), &
csf_i%B_real(sw), weights%dat)
#ifdef DEBUG_
call getDoubleMatrixElement(2, 1, 0, gen_type%L, gen_type%R, &
csf_i%B_real(sw), 1.0_dp, x1_element=stay_mat)
#endif
else
switch_weight = weights%proc%minus(negSwitches(sw), &
csf_i%B_real(sw), weights%dat)
#ifdef DEBUG_
call getDoubleMatrixElement(1, 2, 0, gen_type%L, gen_type%R, &
csf_i%B_real(sw), 1.0_dp, x1_element=stay_mat)
#endif
end if
! update inverse product
! and also get starting contribution
ASSERT(.not. near_zero(stay_mat))
stay_weight = 1.0_dp - calcStayingProb(zero_weight, switch_weight, &
csf_i%B_real(sw))
! and the new startProb is also the non-b=0 branch
start_weight = switch_weight / (zero_weight + switch_weight)
if (t_trunc_guga_pgen .or. &
(t_trunc_guga_pgen_noninits .and. .not. is_init_guga)) then
if (branch_pgen * start_weight / start_weight > trunc_guga_pgen) then
pgen = pgen + orb_pgen * branch_pgen * start_weight / stay_weight
end if
else
pgen = pgen + orb_pgen * branch_pgen * start_weight / stay_weight
end if
end if
end if
end if
! i also need to consider the electron pair picking probability..
pgen = pgen / real(ElecPairs, dp)
! and if the second electron is in a double occupied orbital I have
! to modify it with 2
if (csf_i%stepvector(elecInd) == 3) pgen = pgen * 2.0_dp
end subroutine calc_mixed_start_contr_pgen