subroutine calc_mixed_contr_pgen(ilut, csf_i, t, excitInfo, 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), intent(out) :: pgen
integer :: first, last, deltaB(nSpatOrbs), i, j, k, step1
real(dp) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs), &
zeroWeight, minusWeight, plusWeight, branch_weight, &
above_cpt, below_cpt
type(WeightObj_t) :: weights
logical :: above_flag, below_flag
! also need the pgen contributions from all other index combinations
! shich could lead to this excitation
first = findFirstSwitch(ilut, t, excitInfo%fullStart, excitInfo%fullEnd)
last = findLastSwitch(ilut, t, first, excitInfo%fullEnd)
below_flag = .false.
above_flag = .false.
pgen = 0.0_dp
deltaB = int(csf_i%B_real - calcB_vector_ilut(t(0:nifd)))
do j = last, nSpatOrbs
if (csf_i%Occ_int(j) /= 1) cycle
! calculate the remaining switches once for each (j) but do it
! for the worst case until i = 1
! check if this is the last end needed to consider
if (csf_i%stepvector(j) == 2 .and. csf_i%B_int(j) == 0) then
above_flag = .true.
end if
excitInfo%fullStart = 1
excitInfo%secondStart = 1
excitInfo%fullEnd = j
excitInfo%firstEnd = j
! reinit remainings switches and weights
! i think i could also do a on-the fly switch recalculation..
! so only the weights have to be reinited
call calcRemainingSwitches_excitInfo_double(csf_i, excitInfo, &
posSwitches, negSwitches)
weights = init_doubleWeight(csf_i, j)
! i have to reset the below flag each iteration of j..
below_flag = .false.
do i = first, 1, -1
if (csf_i%Occ_int(i) /= 1) cycle
if (below_flag) exit
! this is the only difference for molecular/hubbard/ueg
! calculations
call calc_orbital_pgen_contr(csf_i, [2 * i, 2 * j], above_cpt, &
below_cpt)
! yes they can, and then this orbital does not contribute to the
! obtained excitaiton -> cycle..
! only from the names.. shouldnt here be below_cpt?
! ah ok below is the below_flag!
! if the bottom stepvector d = 1 and b = 1 there is no
! additional contribution from below, since the x1 matrix
! element is 0
! same if d = 2 and b = 0 for fullstop stepvector
if (csf_i%stepvector(i) == 1 .and. csf_i%B_int(i) == 1) then
below_flag = .true.
end if
zeroWeight = weights%proc%zero(negSwitches(i), posSwitches(i), &
csf_i%B_real(i), weights%dat)
! deal with the start seperately:
if (csf_i%stepvector(i) == 1) then
plusWeight = weights%proc%plus(posSwitches(i), &
csf_i%B_real(i), weights%dat)
if (isOne(t, i)) then
branch_weight = zeroWeight / (zeroWeight + plusWeight)
else
branch_weight = plusWeight / (zeroWeight + plusWeight)
end if
else
minusWeight = weights%proc%minus(negSwitches(i), &
csf_i%B_real(i), weights%dat)
if (isTwo(t, i)) then
branch_weight = zeroWeight / (zeroWeight + minusWeight)
else
branch_weight = minusWeight / (zeroWeight + minusWeight)
end if
end if
! loop over excitation range
! distinguish between different regimes
! if i do it until switch - 1 -> i know that dB = 0 and
! the 2 stepvalues are always the same..
do k = i + 1, first - 1
if (csf_i%Occ_int(k) /= 1) cycle
step1 = csf_i%stepvector(k)
zeroWeight = weights%proc%zero(negSwitches(k), &
posSwitches(k), csf_i%B_real(k), weights%dat)
if (step1 == 1) then
plusWeight = weights%proc%plus(posSwitches(k), &
csf_i%B_real(k), weights%dat)
branch_weight = branch_weight * calcStayingProb( &
zeroWeight, plusWeight, csf_i%B_real(k))
else
minusWeight = weights%proc%minus(negSwitches(k), &
csf_i%B_real(k), weights%dat)
branch_weight = branch_weight * calcStayingProb( &
zeroWeight, minusWeight, csf_i%B_real(k))
end if
end do
! then do first switch site seperately, if (i) is not first
! and what if (i) is first??
if (i /= first) then
step1 = csf_i%stepvector(first)
zeroWeight = weights%proc%zero(negSwitches(first), &
posSwitches(first), csf_i%B_real(first), weights%dat)
if (step1 == 1) then
! i know that step2 = 2
plusWeight = weights%proc%plus(posSwitches(first), &
csf_i%B_real(first), weights%dat)
branch_weight = branch_weight * (1.0_dp - calcStayingProb( &
zeroWeight, plusWeight, csf_i%B_real(first)))
else
! i know that step2 = 1
minusWeight = weights%proc%minus(negSwitches(first), &
csf_i%B_real(first), weights%dat)
branch_weight = branch_weight * (1.0_dp - calcStayingProb( &
zeroWeight, minusWeight, csf_i%B_real(first)))
end if
end if
! loop over the range where switch happened
do k = first + 1, last - 1
! in this region i know, that the matrix element is
! definetly not 0, since otherwise the excitation would
! have been aborted before
! combine stepvalue and deltaB info in select statement
if (csf_i%Occ_int(k) /= 1) cycle
zeroWeight = weights%proc%zero(negSwitches(k), &
posSwitches(k), csf_i%B_real(k), weights%dat)
select case (deltaB(k - 1) + csf_i%stepvector(k))
case (1)
! d=1 + b=0 : 1
plusWeight = weights%proc%plus(posSwitches(k), &
csf_i%B_real(k), weights%dat)
if (isOne(t, k)) then
branch_weight = branch_weight * calcStayingProb( &
zeroWeight, plusWeight, csf_i%B_real(k))
else
branch_weight = branch_weight * (1.0_dp - calcStayingProb( &
zeroWeight, plusWeight, csf_i%B_real(k)))
end if
case (2)
! d=2 + b=0 : 2
minusWeight = weights%proc%minus(negSwitches(k), &
csf_i%B_real(k), weights%dat)
if (isTwo(t, k)) then
branch_weight = branch_weight * calcStayingProb( &
zeroWeight, minusWeight, csf_i%B_real(k))
else
branch_weight = branch_weight * (1.0_dp - calcStayingProb( &
zeroWeight, minusWeight, csf_i%B_real(k)))
end if
case (-1)
! d=1 + b=-2 : -1
minusWeight = weights%proc%minus(negSwitches(k), &
csf_i%B_real(k), weights%dat)
if (isOne(t, k)) then
branch_weight = branch_weight * calcStayingProb(minusWeight, &
zeroWeight, csf_i%B_real(k))
else
branch_weight = branch_weight * (1.0_dp - calcStayingProb( &
minusWeight, zeroWeight, csf_i%B_real(k)))
end if
case (4)
! d=2 + b=2 : 4
zeroWeight = weights%proc%zero(negSwitches(k), &
posSwitches(k), csf_i%B_real(k), weights%dat)
plusWeight = weights%proc%plus(posSwitches(k), &
csf_i%B_real(k), weights%dat)
if (isTwo(t, k)) then
branch_weight = branch_weight * calcStayingProb(plusWeight, &
zeroWeight, csf_i%B_real(k))
else
branch_weight = branch_weight * (1.0_dp - calcStayingProb( &
plusWeight, zeroWeight, csf_i%B_real(k)))
end if
end select
end do
! more efficient to do "last" step seperately, since i have to
! check deltaB value and also have to consider matrix element
! but only of (j) is not last or otherwise already dealt with
if (j /= last) then
if (csf_i%stepvector(last) == 1) then
! then i know step2 = 2 & dB = -2!
zeroWeight = weights%proc%zero(negSwitches(last), &
posSwitches(last), csf_i%B_real(last), weights%dat)
minusWeight = weights%proc%minus(negSwitches(last), &
csf_i%B_real(last), weights%dat)
branch_weight = branch_weight * (1.0_dp - calcStayingProb( &
minusWeight, zeroWeight, csf_i%B_real(last)))
else
! i know step2 == 1 and dB = +2
zeroWeight = weights%proc%zero(negSwitches(last), &
posSwitches(last), csf_i%B_real(last), weights%dat)
plusWeight = weights%proc%plus(posSwitches(last), &
csf_i%B_real(last), weights%dat)
branch_weight = branch_weight * (1.0_dp - calcStayingProb( &
plusWeight, zeroWeight, csf_i%B_real(last)))
end if
end if
! then do remaining top range, where i know stepvalues are
! the same again and dB = 0 always!
do k = last + 1, j - 1
if (csf_i%Occ_int(k) /= 1) cycle
step1 = csf_i%stepvector(k)
zeroWeight = weights%proc%zero(negSwitches(k), &
posSwitches(k), csf_i%B_real(k), weights%dat)
if (step1 == 1) then
! i know step2 = 1 als
plusWeight = weights%proc%plus(posSwitches(k), &
csf_i%B_real(k), weights%dat)
branch_weight = branch_weight * calcStayingProb( &
zeroWeight, plusWeight, csf_i%B_real(k))
else
minusWeight = weights%proc%minus(negSwitches(k), &
csf_i%B_real(k), weights%dat)
branch_weight = branch_weight * calcStayingProb( &
zeroWeight, minusWeight, csf_i%B_real(k))
end if
end do
if (t_trunc_guga_pgen .or. &
(t_trunc_guga_pgen_noninits .and. .not. is_init_guga)) then
if (branch_weight < trunc_guga_pgen) then
branch_weight = 0.0_dp
end if
end if
! add up pgen contributions..
pgen = pgen + (below_cpt + above_cpt) * branch_weight
! check if i deal with that correctly...
if (below_flag) exit
end do
! todo: i cant use tthat like that.. or else some combinations
! of i and j get left out! i have to reinit it somehow..
! not yet sure how..
if (above_flag) exit
end do
! multiply by always same probability to pick the 2 electrons
if (.not. (t_heisenberg_model .or. t_tJ_model)) then
pgen = pgen / real(ElecPairs, dp)
end if
end subroutine calc_mixed_contr_pgen