subroutine calc_mixed_start_contr_sym(ilut, csf_i, t, excitInfo, branch_pgen, &
pgen, 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(inout) :: excitInfo
real(dp), intent(inout) :: branch_pgen
real(dp), intent(out) :: pgen
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_sym"
integer :: sw, i, st, se, step, en, elecInd, holeInd
real(dp) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs), mat_ele, &
new_pgen, zero_weight, switch_weight, stay_mat, start_mat, bot_cont, &
orb_pgen, start_weight, stay_weight
type(WeightObj_t) :: weights
logical :: below_flag
integer(int_rdm), allocatable :: tmp_rdm_ind(:)
real(dp), allocatable :: tmp_rdm_mat(:)
HElement_t(dp), allocatable :: temp_int
integer :: rdm_count, max_num_rdm
logical :: rdm_flag, test_skip
test_skip = .false.
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)
! 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
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 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.
! check if orb_pgen is non-zero
if (.not. test_skip) then
if (near_zero(orb_pgen) .and. (.not. rdm_flag)) then
temp_int = (get_umat_el(i, holeInd, elecInd, i) &
+ get_umat_el(holeInd, i, i, elecInd)) / 2.0_dp
if (.not. near_zero(temp_int) .and. near_zero(orb_pgen)) then
print *, "start 1 orb_pgen 0, integral: ", temp_int
end if
! still have to update matrix element, even if 0 pgen
mat_ele = mat_ele * stay_mat
cycle
end if
end if
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
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
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 + 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
if (test_skip) then
if (near_zero(orb_pgen)) then
print *, "stay_weight(1): ", stay_weight
print *, "i, st: ", i, st
end if
end if
! update new_pgen for next cycle
new_pgen = stay_weight * new_pgen
! 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
! 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
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
! check if orb_pgen is non-zero
! still have to update matrix element in this case..
! so do the cycle only afterwards..
if (.not. test_skip) then
if (near_zero(orb_pgen) .and. (.not. rdm_flag)) then
temp_int = (get_umat_el(i, holeInd, elecInd, i) &
+ get_umat_el(holeInd, i, i, elecInd)) / 2.0_dp
if (.not. near_zero(temp_int) .and. near_zero(orb_pgen)) then
print *, "start 2 orb_pgen 0, integral: ", temp_int
end if
cycle
end if
end if
! 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)
! 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
! and update probWeight
! i think i should not put in the intermediate start_weight
! permanently..
if (test_skip) then
if (near_zero(orb_pgen)) then
print *, "stay_weight(2): ", stay_weight
print *, "i, st, sw:", i, st, sw
end if
end if
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)
temp_int = (get_umat_el(sw, holeInd, elecInd, sw) &
+ get_umat_el(holeInd, sw, sw, elecInd)) / 2.0_dp
if (near_zero(orb_pgen) .and. .not. near_zero(temp_int)) then
print *, "start 3 orb_pgen 0, integral: ", temp_int
end if
if (.not. near_zero(orb_pgen) .or. rdm_flag .or. test_skip) 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)
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
switch_weight = weights%proc%minus(negSwitches(sw), &
csf_i%B_real(sw), weights%dat)
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
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
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_sym