pure subroutine calc_mixed_start_contr_integral(ilut, csf_i, t, excitInfo, &
integral, rdm_ind, rdm_mat)
integer(n_int), intent(in) :: ilut(0:nifguga), t(0:nifguga)
type(CSF_Info_t), intent(in) :: csf_i
type(ExcitationInformation_t), intent(in), value :: excitInfo
HElement_t(dp), intent(out) :: integral
integer(int_rdm), intent(out), allocatable, optional :: rdm_ind(:)
real(dp), intent(out), allocatable, optional :: rdm_mat(:)
character(*), parameter :: this_routine = "calc_mixed_start_contr_integral"
integer :: st, se, en, elecInd, holeInd, sw, step, i, max_num_rdm
integer :: rdm_count
real(dp) :: bot_cont
logical :: below_flag, rdm_flag
real(dp) :: mat_ele, start_mat, stay_mat
integer(int_rdm), allocatable :: tmp_rdm_ind(:)
real(dp), allocatable :: tmp_rdm_mat(:)
if (present(rdm_ind) .or. present(rdm_mat)) then
ASSERT(present(rdm_ind) .and. present(rdm_mat))
rdm_flag = .true.
else
rdm_flag = .false.
end if
! whats different here?? what do i have to consider? and how to optimize?
! to make it most similar to the full-start into full-stop calc.
! i could loop from the first switch downwards and stop at
! a d = 1, b = 1 stepvalue and definetly unify pgen and integral
! calculation!
! to similary reuse the already calculated quantities loop from
! switch to start to 1
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)
if (rdm_flag) then
max_num_rdm = sw
allocate(tmp_rdm_ind(max_num_rdm), source=0_int_rdm)
allocate(tmp_rdm_mat(max_num_rdm), source=0.0_dp)
rdm_count = 0
end if
! what can i precalculate beforehand?
step = csf_i%stepvector(st)
integral = h_cast(0.0_dp)
if (step == 1) then
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))
else
bot_cont = -sqrt(2.0_dp / ((csf_i%B_real(st) - 1.0_dp) * &
(csf_i%B_real(st) + 1.0_dp)))
end if
else
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)))
else
bot_cont = -Root2 * sqrt((csf_i%B_real(st) + 3.0_dp) / &
(csf_i%B_real(st) + 1.0_dp))
end if
end if
! 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
mat_ele = 1.0_dp
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 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.
! 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
integral = integral + start_mat * mat_ele * &
(get_umat_el(i, holeInd, elecInd, i) &
+ get_umat_el(holeInd, i, i, elecInd)) / 2.0_dp
if (rdm_flag) then
rdm_count = rdm_count + 1
tmp_rdm_ind(rdm_count) = contract_2_rdm_ind(i, elecInd, holeInd, i)
tmp_rdm_mat(rdm_count) = start_mat * mat_ele * bot_cont
end if
end if
if (below_flag) exit
! also update matrix element on the fly
mat_ele = stay_mat * mat_ele
end do
! and update matrix element finally with bottom contribution
integral = integral * bot_cont
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=mat_ele)
! and calc. x1^-1
! keep tempWweight as the running matrix element which gets updated
! every iteration
! for rdms (in this current setup) I need to make a dummy
! output if sw == st)
if (rdm_flag .and. sw == st) then
rdm_count = rdm_count + 1
tmp_rdm_ind(rdm_count) = contract_2_rdm_ind(sw, elecInd, holeInd, sw)
tmp_rdm_mat(rdm_count) = 1.0_dp
end if
if (.not. near_zero(abs(mat_ele))) then
mat_ele = 1.0_dp / mat_ele
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
step = csf_i%stepvector(i)
! update inverse product
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))
mat_ele = mat_ele / stay_mat
! and also get starting contribution
call getDoubleMatrixElement(step, step, -1, gen_type%L, gen_type%R, &
csf_i%B_real(i), 1.0_dp, x1_element=start_mat)
! because the rest of the matrix element is still the same in
! both cases...
if (.not. near_zero(start_mat)) then
integral = integral + mat_ele * start_mat * &
(get_umat_el(holeInd, i, i, elecInd) + &
get_umat_el(i, holeInd, elecInd, i)) / 2.0_dp
if (rdm_flag) then
rdm_count = rdm_count + 1
tmp_rdm_ind(rdm_count) = contract_2_rdm_ind(i, elecInd, holeInd, i)
tmp_rdm_mat(rdm_count) = start_mat * mat_ele
end if
end if
end do
! handle switch seperately (but only if switch > start)
if (sw > st) then
step = csf_i%stepvector(sw)
! on the switch the original probability is:
if (step == 1) then
call getDoubleMatrixElement(2, 1, 0, gen_type%L, gen_type%R, &
csf_i%B_real(sw), 1.0_dp, x1_element=stay_mat)
call getDoubleMatrixElement(2, 1, -1, gen_type%L, gen_type%R, &
csf_i%B_real(sw), 1.0_dp, x1_element=start_mat)
else
call getDoubleMatrixElement(1, 2, 0, gen_type%L, gen_type%R, &
csf_i%B_real(sw), 1.0_dp, x1_element=stay_mat)
call getDoubleMatrixElement(1, 2, -1, gen_type%L, gen_type%R, &
csf_i%B_real(sw), 1.0_dp, x1_element=start_mat)
end if
! update inverse product
! and also get starting contribution
ASSERT(.not. near_zero(stay_mat))
mat_ele = mat_ele * start_mat / stay_mat
! because the rest of the matrix element is still the same in
! both cases...
integral = integral + mat_ele * (get_umat_el(holeInd, sw, sw, elecInd) + &
get_umat_el(sw, holeInd, elecInd, sw)) / 2.0_dp
if (rdm_flag) then
rdm_count = rdm_count + 1
tmp_rdm_ind(rdm_count) = contract_2_rdm_ind(sw, elecInd, holeInd, sw)
tmp_rdm_mat(rdm_count) = mat_ele
end if
end if
end if
if (present(rdm_mat)) then
allocate(rdm_ind(rdm_count), source=tmp_rdm_ind(1:rdm_count))
allocate(rdm_mat(rdm_count), source=tmp_rdm_mat(1:rdm_count))
deallocate(tmp_rdm_ind)
deallocate(tmp_rdm_mat)
end if
end subroutine calc_mixed_start_contr_integral