pure subroutine calc_integral_contribution_single(csf_i, csf_j, i, j, st, en, integral)
! calculates the double-excitaiton contribution to a single excitation
type(CSF_Info_t), intent(in) :: csf_i, csf_j
integer, intent(in) :: i, j, st, en
HElement_t(dp), intent(inout) :: integral
real(dp) :: botCont, topCont, tempWeight, prod
integer :: iO, jO, step
! calculate the bottom contribution depending on the excited stepvalue
select case (csf_i%stepvector(st))
case (0)
! this implicates a raising st:
if (csf_j%stepvector(st) == 1) then
call getDoubleMatrixElement(1, 0, 0, gen_type%L, gen_type%R, csf_i%B_real(st), &
1.0_dp, x1_element=botCont)
else
call getDoubleMatrixElement(2, 0, 0, gen_type%L, gen_type%R, csf_i%B_real(st), &
1.0_dp, x1_element=botCont)
end if
case (3)
! implies lowering st
if (csf_j%stepvector(st) == 1) then
! need tA(0,2)
botCont = funA_0_2overR2(csf_i%B_real(st))
else
! need -tA(2,0)
botCont = minFunA_2_0_overR2(csf_i%B_real(st))
end if
case (1)
botCont = funA_m1_1_overR2(csf_i%B_real(st))
! check which generator
if (csf_j%stepvector(st) == 0) botCont = -botCont
case (2)
botCont = funA_3_1_overR2(csf_i%B_real(st))
if (csf_j%stepvector(st) == 3) botCont = -botCont
end select
! do top contribution also already
select case (csf_i%stepvector(en))
case (0)
if (csf_j%stepvector(en) == 1) then
topCont = funA_2_0_overR2(csf_i%B_real(en))
else
topCont = minFunA_0_2_overR2(csf_i%B_real(en))
end if
case (3)
if (csf_j%stepvector(en) == 1) then
topCont = minFunA_2_0_overR2(csf_i%B_real(en))
else
topCont = funA_0_2overR2(csf_i%B_real(en))
end if
case (1)
topCont = funA_2_0_overR2(csf_i%B_real(en))
if (csf_j%stepvector(en) == 3) topCont = -topCont
case (2)
topCont = funA_0_2overR2(csf_i%B_real(en))
if (csf_j%stepvector(en) == 0) topCont = -topCont
end select
! depending on i and j calulate the corresponding single and double
! integral weights and check if they are non-zero...
! gets quite involved... :( need to loop over all orbitals
! have to reset prod inside the loop each time!
do iO = 1, st - 1
! no contribution if not occupied.
if (csf_i%stepvector(iO) == 0) cycle
! else it gets a contrbution weighted with orbital occupation
! first easy part:
integral = integral + get_umat_el(i, iO, j, iO) * csf_i%Occ_real(iO)
! also easy is the non-product involving part...
integral = integral - get_umat_el(i, iO, iO, j) * &
csf_i%Occ_real(iO) / 2.0_dp
! the product part is annoying actually... but doesnt help... todo
! have to do a second loop for the product
! for the first loop iteration i have to access the mixed fullstart
! elements, with a deltaB = -1 value!!
! think about how to implement that !
! do i have to do anything for a d = 3 ? since x1-element is 0
! in this case anyway.. there should not be an influence.
! also if its a d = 2 with b = 0 the matrix element is also 0
if (csf_i%stepvector(iO) == 3 .or. csf_i%B_int(iO) == 0) cycle
step = csf_i%stepvector(iO)
call getDoubleMatrixElement(step, step, -1, gen_type%L, gen_type%R, csf_i%B_real(iO), &
1.0_dp, x1_element=prod)
! and then do the remaining:
do jO = iO + 1, st - 1
! need the stepvalue entries to correctly access the mixed
! generator matrix elements
step = csf_i%stepvector(jO)
call getDoubleMatrixElement(step, step, 0, gen_type%L, gen_type%R, &
csf_i%B_real(jO), 1.0_dp, x1_element=tempWeight)
prod = prod * tempWeight
end do
prod = prod * botCont
integral = integral + get_umat_el(i, iO, iO, j) * prod
end do
! also have to calc the top and bottom contribution to the product terms
! BUT they also depend on the excitation -> so i should only calculate
! these terms after i started the excitation! todo
! at the start of excittaion also certain values.
! need for second sum term the occupation of the excited state... ->
! need generator type considerations.. todo
! at i the occupation of the excited state is necessary ->
! which. independently means
! did some stupid double counting down there...
! still something wrong down there...
! i should be able to formulate that in terms of st and en..
integral = integral + get_umat_el(i, i, j, i) * csf_i%Occ_real(i)
integral = integral + get_umat_el(i, j, j, j) * (csf_i%Occ_real(j) - 1.0_dp)
! have to reset prod here!!!
! also do the same loop on the orbitals above the fullEnd. to get
! the double contribution
do iO = en + 1, nSpatOrbs
! do stuff
if (csf_i%stepvector(iO) == 0) cycle
integral = integral + get_umat_el(i, iO, j, iO) * csf_i%Occ_real(iO)
integral = integral - get_umat_el(i, iO, iO, j) * csf_i%Occ_real(iO) / 2.0_dp
! not necessary to do it for d = 3 or b = 1, d=1 end value! since
! top matrix element 0 in this case
if (csf_i%stepvector(iO) == 3 .or. (csf_i%B_int(iO) == 1 &
.and. csf_i%stepvector(iO) == 1)) cycle
! have to reset prod every loop
prod = 1.0_dp
do jO = en + 1, iO - 1
! do stuff
step = csf_i%stepvector(jO)
call getDoubleMatrixElement(step, step, 0, gen_type%L, gen_type%R, csf_i%B_real(jO), &
1.0_dp, x1_element=tempWeight)
prod = prod * tempWeight
end do
! have to seperately access the top most mixed full-stop
step = csf_i%stepvector(iO)
call getMixedFullStop(step, step, 0, csf_i%B_real(iO), &
x1_element=tempWeight)
prod = prod * tempWeight
prod = prod * topCont
integral = integral + get_umat_el(i, iO, iO, j) * prod
end do
end subroutine calc_integral_contribution_single