subroutine calc_mixed_end_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_end_contr_pgen"
integer :: st, se, en, step, sw, elecInd, holeInd, i, j
real(dp) :: top_cont, stay_mat, end_mat, orb_pgen, new_pgen, &
posSwitches(nSpatOrbs), negSwitches(nSpatOrbs), &
tmp_pos(nSpatOrbs), tmp_neg(nSpatOrbs)
logical :: above_flag
type(BranchWeightArr_t) :: weight_funcs(nSpatOrbs)
type(WeightObj_t) :: weights
! do as much stuff as possible beforehand
st = excitInfo%fullStart
se = excitInfo%secondStart
en = excitInfo%fullEnd
if (excitInfo%typ == excit_type%fullstop_L_to_R) then
elecInd = st
holeInd = se
else if (excitInfo%typ == excit_type%fullstop_R_to_L) then
elecInd = se
holeInd = st
else
call stop_all(this_routine, "should not be here!")
end if
! also here i didn't consider the actual end contribution or? ...
call calc_orbital_pgen_contrib_end(csf_i, [2 * elecInd, 2 * en], &
holeInd, orb_pgen)
pgen = orb_pgen * branch_pgen
step = csf_i%stepvector(en)
sw = findLastSwitch(ilut, t, se, en)
call calcRemainingSwitches_excitInfo_double(csf_i, excitInfo, posSwitches, &
negSwitches)
! need temporary switch arrays for more efficiently recalcing
! weights
tmp_pos = posSwitches
tmp_neg = negSwitches
! after last switch only dB = 0 branches! consider that
call setup_weight_funcs(t, csf_i, st, se, weight_funcs)
if (en < nSpatOrbs) then
select case (step)
case (1)
if (isOne(t, en)) then
top_cont = -Root2 * sqrt((csf_i%B_real(en) + 2.0_dp) / &
csf_i%B_real(en))
else
top_cont = -Root2 / sqrt(csf_i%B_real(en) * (csf_i%B_real(en) + 2.0_dp))
end if
case (2)
if (isOne(t, en)) then
top_cont = -Root2 / sqrt(csf_i%B_real(en) * (csf_i%B_real(en) + 2.0_dp))
else
top_cont = Root2 * sqrt(csf_i%B_real(en) / &
(csf_i%B_real(en) + 2.0_dp))
end if
case default
call stop_all(this_routine, "wrong stepvalues!")
end select
if (.not. near_zero(top_cont)) then
above_flag = .false.
! to avoid to recalc. remaining switches all the time
! just increment them correctly
if (step == 1) then
tmp_neg(se:en - 1) = tmp_neg(se:en - 1) + 1.0_dp
else
tmp_pos(se:en - 1) = tmp_pos(se:en - 1) + 1.0_dp
end if
do i = en + 1, nSpatOrbs
if (csf_i%Occ_int(i) /= 1) cycle
! then check if thats the last step
if (csf_i%stepvector(i) == 2 .and. csf_i%B_int(i) == 0) then
above_flag = .true.
end if
! then calc. orbital probability
call calc_orbital_pgen_contrib_end(csf_i, [2 * elecInd, 2 * i], &
holeInd, orb_pgen)
! should be able to do that without second loop too!
! figure out!
step = csf_i%stepvector(i)
call getDoubleMatrixElement(step, step, 0, gen_type%L, gen_type%R, csf_i%B_real(i), &
1.0_dp, x1_element=stay_mat)
call getMixedFullStop(step, step, 0, csf_i%B_real(i), &
x1_element=end_mat)
! this check should never be true, but just to be sure
if (near_zero(stay_mat)) above_flag = .true.
if (.not. near_zero(end_mat)) then
! also only recalc. pgen if matrix element is not 0
excitInfo%fullEnd = i
excitInfo%firstEnd = i
weights = init_semiStartWeight(csf_i, se, i, tmp_neg(se), &
tmp_pos(se), csf_i%B_real(se))
new_pgen = 1.0_dp
! deal with the start and semi-start seperately
if (csf_i%Occ_int(st) /= 1) then
new_pgen = new_pgen * weight_funcs(st)%ptr(weights, &
csf_i%B_real(st), tmp_neg(st), tmp_pos(st))
end if
do j = st + 1, se - 1
! can and do i have to cycle here if its not
! singly occupied??
if (csf_i%Occ_int(j) /= 1) cycle
new_pgen = new_pgen * weight_funcs(j)%ptr(weights, &
csf_i%B_real(j), tmp_neg(j), tmp_pos(j))
end do
! then need to reinit double weight
weights = weights%ptr
! and also with the semi-start
if (csf_i%Occ_int(se) /= 1) then
new_pgen = new_pgen * weight_funcs(se)%ptr(weights, &
csf_i%B_real(se), tmp_neg(se), tmp_pos(se))
end if
do j = se + 1, i - 1
if (csf_i%Occ_int(j) /= 1) cycle
new_pgen = new_pgen * weight_funcs(j)%ptr(weights, &
csf_i%B_real(j), tmp_neg(j), tmp_pos(j))
end do
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
pgen = pgen + new_pgen * orb_pgen
end if
if (above_flag) exit
! update the switches
if (csf_i%stepvector(i) == 1) then
tmp_neg(se:i - 1) = tmp_neg(se:i - 1) + 1.0_dp
else
tmp_pos(se:i - 1) = tmp_pos(se:i - 1) + 1.0_dp
end if
end do
end if
end if
if (sw < en) then
step = csf_i%stepvector(en)
! have to change the switches before the first cycle:
! but for cycling backwards, thats not so easy.. need todo
do i = en - 1, sw + 1, -1
if (csf_i%Occ_int(i) /= 1) cycle
! get orbital pgen
call calc_orbital_pgen_contrib_end(csf_i, [2 * elecInd, 2 * i], &
holeInd, orb_pgen)
if (csf_i%stepvector(i) == 1) then
! by looping in this direction i have to reduce
! the number of switches at the beginning
! but only to the left or??
! i think i have to rethink that.. thats not so easy..
negSwitches(se:i - 1) = negSwitches(se:i - 1) - 1.0_dp
else
posSwitches(se:i - 1) = posSwitches(se:i - 1) - 1.0_dp
end if
step = csf_i%stepvector(i)
call getMixedFullStop(step, step, 0, csf_i%B_real(i), x1_element=end_mat)
if (.not. near_zero(end_mat)) then
! only recalc. pgen if matrix element is not 0
excitInfo%fullEnd = i
excitInfo%firstEnd = i
weights = init_semiStartWeight(csf_i, se, i, negSwitches(se), &
posSwitches(se), csf_i%B_real(se))
new_pgen = 1.0_dp
! deal with the start and semi-start seperately
if (csf_i%Occ_int(st) /= 1) then
new_pgen = new_pgen * weight_funcs(st)%ptr(weights, &
csf_i%B_real(st), negSwitches(st), posSwitches(st))
end if
do j = st + 1, se - 1
if (csf_i%Occ_int(j) /= 1) cycle
new_pgen = new_pgen * weight_funcs(j)%ptr(weights, &
csf_i%B_real(j), negSwitches(j), posSwitches(j))
end do
! then need to reinit double weight
weights = weights%ptr
! and also with the semi-start
if (csf_i%Occ_int(se) /= 1) then
new_pgen = new_pgen * weight_funcs(se)%ptr(weights, &
csf_i%B_real(se), negSwitches(se), posSwitches(se))
end if
do j = se + 1, i - 1
if (csf_i%Occ_int(j) /= 1) cycle
new_pgen = new_pgen * weight_funcs(j)%ptr(weights, &
csf_i%B_real(j), negSwitches(j), posSwitches(j))
end do
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
pgen = pgen + new_pgen * orb_pgen
end if
end do
! deal with switch specifically:
! figure out orbital pgen
call calc_orbital_pgen_contrib_end(csf_i, [2 * elecInd, 2 * sw], &
holeInd, orb_pgen)
if (.not. near_zero(orb_pgen)) then
step = csf_i%stepvector(sw)
if (step == 1) then
! then a -2 branch arrived!
call getMixedFullStop(2, 1, -2, csf_i%B_real(sw), x1_element=end_mat)
! also reduce negative switches then
! only everything to the left or?
negSwitches(se:sw - 1) = negSwitches(se:sw - 1) - 1.0_dp
else
! +2 branch arrived!
call getMixedFullStop(1, 2, 2, csf_i%B_real(sw), x1_element=end_mat)
! reduce positive switchtes otherwise
posSwitches(se:sw - 1) = posSwitches(se:sw - 1) - 1.0_dp
end if
! loop to get correct pgen
new_pgen = 1.0_dp
weights = init_semiStartWeight(csf_i, se, sw, negSwitches(se), &
posSwitches(se), csf_i%B_real(se))
! deal with the start and semi-start seperately
if (csf_i%Occ_int(st) /= 1) then
new_pgen = new_pgen * weight_funcs(st)%ptr(weights, &
csf_i%B_real(st), negSwitches(st), posSwitches(st))
end if
do j = st + 1, se - 1
if (csf_i%Occ_int(j) /= 1) cycle
new_pgen = new_pgen * weight_funcs(j)%ptr(weights, &
csf_i%B_real(j), negSwitches(j), posSwitches(j))
end do
weights = weights%ptr
! and also with the semi-start
if (csf_i%Occ_int(se) /= 1) then
new_pgen = new_pgen * weight_funcs(se)%ptr(weights, &
csf_i%B_real(se), negSwitches(se), posSwitches(se))
end if
do j = se + 1, sw - 1
if (csf_i%Occ_int(j) /= 1) cycle
new_pgen = new_pgen * weight_funcs(j)%ptr(weights, &
csf_i%B_real(j), negSwitches(j), posSwitches(j))
end do
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
pgen = pgen + new_pgen * orb_pgen
end if
end if
pgen = pgen / real(ElecPairs, dp)
if (csf_i%stepvector(elecInd) == 3) pgen = pgen * 2.0_dp
end subroutine calc_mixed_end_contr_pgen