pure subroutine calc_mixed_contr_integral(ilut, csf_i, t, start, ende, integral, &
rdm_ind, rdm_mat)
integer(n_int), intent(in) :: ilut(0:nifguga), t(0:nifguga)
type(CSF_Info_t), intent(in) :: csf_i
integer, intent(in) :: start, ende
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_contr_integral"
real(dp) :: inter, tempWeight, tempWeight_1
HElement_t(dp) :: temp_int
integer :: i, j, k, step1, step2, bVector(nSpatOrbs), first, last
logical :: rdm_flag
integer :: max_num_rdm, rdm_count
integer(int_rdm), allocatable :: tmp_rdm_ind(:)
real(dp), allocatable :: tmp_rdm_mat(:)
real(dp) :: tmp_mat
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 it differently... since its always a mixed double overlap region
! and only 2 indices are involved.. just recalc all possible
! matrix elements in this case..
! so first determine the first and last switches and calculate
! the overlap matrix element in this region, since atleast thats
! always the same
first = findFirstSwitch(ilut, t, start, ende)
last = findLastSwitch(ilut, t, first, ende)
if (rdm_flag) then
max_num_rdm = first * (nSpatOrbs - last + 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
! calc. the intermediate matrix element..
! but what is the deltaB value inbetween? calculate it on the fly..
bVector = int(calcB_vector_ilut(ilut(0:nifd)) - calcB_vector_ilut(t(0:nifd)))
inter = 1.0_dp
integral = h_cast(0.0_dp)
do i = first + 1, last - 1
if (csf_i%Occ_int(i) /= 1) cycle
step1 = csf_i%stepvector(i)
step2 = getStepvalue(t, i)
call getDoubleMatrixElement(step2, step1, bVector(i - 1), gen_type%L, gen_type%R, &
csf_i%B_real(i), 1.0_dp, x1_element=tempWeight)
inter = inter * tempWeight
end do
! and then add up all the possible integral contribs
do i = 1, first
if (csf_i%Occ_int(i) /= 1) cycle
do j = last, nSpatOrbs
if (csf_i%Occ_int(j) /= 1) cycle
temp_int = (get_umat_el(i, j, j, i) + get_umat_el(j, i, i, j)) / 2.0_dp
if ((.not. rdm_flag) .and. near_zero(temp_int)) cycle
! get the starting matrix element
step1 = csf_i%stepvector(i)
step2 = getStepvalue(t, i)
call getDoubleMatrixElement(step2, step1, -1, gen_type%L, gen_type%R, &
csf_i%B_real(i), 1.0_dp, x1_element=tempWeight)
! then calc. the product:
do k = i + 1, first
if (csf_i%Occ_int(k) /= 1) cycle
step1 = csf_i%stepvector(k)
step2 = getStepvalue(t, k)
! only 0 branch here
call getDoubleMatrixElement(step2, step1, 0, gen_type%L, gen_type%R, &
csf_i%B_real(k), 1.0_dp, x1_element=tempWeight_1)
tempWeight = tempWeight * tempWeight_1
end do
do k = last, j - 1
if (csf_i%Occ_int(k) /= 1) cycle
step1 = csf_i%stepvector(k)
step2 = getStepvalue(t, k)
call getDoubleMatrixElement(step2, step1, bVector(k - 1), gen_type%L, gen_type%R, &
csf_i%B_real(k), 1.0_dp, x1_element=tempWeight_1)
tempWeight = tempWeight * tempWeight_1
end do
! and then do the the end value at j
step1 = csf_i%stepvector(j)
step2 = getStepvalue(t, j)
call getMixedFullStop(step2, step1, bVector(j - 1), csf_i%B_real(j), &
x1_element=tempWeight_1)
! and multiply and add up all contribution elements
integral = integral + tempWeight * tempWeight_1 * inter * temp_int
if (rdm_flag) then
tmp_mat = tempWeight * tempWeight_1 * inter
if (.not. near_zero(tmp_mat)) then
rdm_count = rdm_count + 1
tmp_rdm_ind(rdm_count) = contract_2_rdm_ind(i, j, j, i)
tmp_rdm_mat(rdm_count) = tmp_mat
end if
end if
end do
end do
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))
end if
end subroutine calc_mixed_contr_integral