subroutine calc_mixed_end_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_end_contr_sym"
integer :: st, se, en, step, sw, elecInd, holeInd, i, j
real(dp) :: top_cont, mat_ele, 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
integer(int_rdm), allocatable :: tmp_rdm_ind(:)
real(dp), allocatable :: tmp_rdm_mat(:)
logical :: rdm_flag, test_skip
integer :: rdm_count, max_num_rdm
test_skip = .false.
if (present(rdm_ind) .or. present(rdm_mat)) then
ASSERT(present(rdm_ind))
ASSERT(present(rdm_mat))
rdm_flag = .true.
else
rdm_flag = .false.
end if
! 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
integral = h_cast(0.0_dp)
! 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)
if (rdm_flag) then
max_num_rdm = (nSpatOrbs - sw + 1)
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
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.
mat_ele = 1.0_dp
! 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)
if (.not. test_skip) then
if (near_zero(orb_pgen) .and. (.not. rdm_flag)) then
! still have to update the switches before cycling
! 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
! also have to update the matrix element, even if
! the orb pgen is 0
mat_ele = mat_ele * stay_mat
cycle
end if
end if
! 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
integral = integral + end_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) = top_cont * end_mat * mat_ele
end if
! 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
! otherwise update your running pgen and matrix element vars
mat_ele = mat_ele * stay_mat
! 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
integral = integral * top_cont
end if
end if
if (rdm_flag .and. sw == en) 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 (sw < en) then
step = csf_i%stepvector(en)
! inverse fullstop matrix element
call getMixedFullStop(step, step, 0, csf_i%B_real(en), x1_element=mat_ele)
ASSERT(.not. near_zero(mat_ele))
mat_ele = 1.0_dp / mat_ele
! 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)
! 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)
call getMixedFullStop(step, step, 0, csf_i%B_real(i), x1_element=end_mat)
! update matrix element
ASSERT(.not. near_zero(stay_mat))
mat_ele = mat_ele / stay_mat
! dont i still have to atleast update the matrix element
! even if the orbital pgen is 0??
if (.not. test_skip) then
if (near_zero(orb_pgen) .and. (.not. rdm_flag)) cycle
end if
if (.not. near_zero(end_mat)) then
integral = integral + end_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) = end_mat * mat_ele
end if
! 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) .or. rdm_flag .or. test_skip) then
step = csf_i%stepvector(sw)
if (step == 1) then
! then a -2 branch arrived!
call getDoubleMatrixElement(2, 1, -2, gen_type%L, gen_type%R, csf_i%B_real(sw), &
1.0_dp, x1_element=stay_mat)
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 getDoubleMatrixElement(1, 2, 2, gen_type%L, gen_type%R, csf_i%B_real(sw), &
1.0_dp, x1_element=stay_mat)
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
ASSERT(.not. near_zero(stay_mat))
mat_ele = mat_ele * end_mat / stay_mat
integral = integral + mat_ele * (get_umat_el(sw, holeInd, elecInd, sw) + &
get_umat_el(holeInd, sw, sw, elecInd)) / 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
! 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
if (rdm_flag) 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_end_contr_sym