pure subroutine calc_mixed_end_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(inout) :: 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_end_contr_integral"
integer :: st, se, en, step, sw, elecInd, holeInd, i
real(dp) :: top_cont, mat_ele, stay_mat, end_mat
logical :: above_flag
integer(int_rdm), allocatable :: tmp_rdm_ind(:)
real(dp), allocatable :: tmp_rdm_mat(:)
logical :: rdm_flag
integer :: rdm_count, max_num_rdm
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)
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
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
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
! 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)
! 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
end if
if (above_flag) exit
! otherwise update your running pgen and matrix element vars
mat_ele = mat_ele * stay_mat
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
do i = en - 1, sw + 1, -1
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)
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
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
end if
end do
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)
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)
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
end if
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_integral