pure subroutine calc_fullstart_mixed_ex(ilutI, csf_i, ilutJ, csf_j, excitInfo, mat_ele, &
t_hamil, rdm_ind, rdm_mat)
integer(n_int), intent(in) :: ilutI(0:niftot), ilutJ(0:niftot)
type(CSF_Info_t), intent(in) :: csf_i, csf_j
type(ExcitationInformation_t), intent(inout) :: excitInfo
HElement_t(dp), intent(out) :: mat_ele
logical, intent(in), optional :: t_hamil
integer(int_rdm), intent(out), allocatable, optional :: rdm_ind(:)
real(dp), intent(out), allocatable, optional :: rdm_mat(:)
integer :: step1, step2, db, i
real(dp) :: bVal, temp_mat0, temp_mat1, temp_x0, temp_x1
HElement_t(dp) :: integral
integer(n_int) :: tmp_I(0:nifguga), tmp_J(0:nifguga)
logical :: t_hamil_
integer :: st, en, se, delta_b(nSpatOrbs)
real(dp) :: guga_mat
def_default(t_hamil_, t_hamil, .true.)
! set defaults for early exit
mat_ele = h_cast(0.0_dp)
delta_b = csf_i%B_int - csf_j%B_int
! i can not associate to all stuff of excitInfo, since it will
! get changed later on..
st = excitInfo%fullStart
en = excitInfo%fullEnd
se = excitInfo%firstEnd
associate (ii => excitInfo%i, jj => excitInfo%j, kk => excitInfo%k, &
ll => excitInfo%l, gen => excitInfo%lastGen, typ => excitInfo%typ)
if (any(abs(delta_b(st:se - 1)) > 2) .or. &
any(abs(delta_b(se:en)) > 1)) return
! do the full-start, and i know here that it is singly occupied
step1 = csf_i%stepvector(st)
step2 = csf_j%stepvector(st)
! to indicate the mixed fullstart matrix elements in the routine!
db = -1
bVal = csf_i%b_real(st)
call getDoubleMatrixElement(step2, step1, -1, gen_type%L, gen_type%R, bVal, 1.0_dp, &
temp_x0, temp_x1)
if (near_zero(temp_x0) .and. near_zero(temp_x1)) return
! then do the double overlap range
do i = st + 1, se
step1 = csf_i%stepvector(i)
step2 = csf_j%stepvector(i)
db = delta_b(i - 1)
bVal = csf_i%b_real(i)
call getDoubleMatrixElement(step2, step1, db, gen_type%L, gen_type%R, bVal, 1.0_dp, &
temp_mat0, temp_mat1)
temp_x0 = temp_x0 * temp_mat0
temp_x1 = temp_x1 * temp_mat1
if (near_zero(temp_x1)) return
end do
! i think here i should also check if the x0 matrix element is non-zero..
if (.not. near_zero(temp_x0)) return
guga_mat = temp_x1
! do single range
do i = se + 1, en
step1 = csf_i%stepvector(i)
step2 = csf_j%stepvector(i)
db = delta_b(i - 1)
bVal = csf_i%b_real(i)
guga_mat = guga_mat * getSingleMatrixElement(step2, step1, db, gen, bVal)
if (near_zero(guga_mat)) return
end do
! then also reuse the stochastic routines to get the integral contribs
call convert_ilut_toGUGA(ilutI, tmp_I)
call convert_ilut_toGUGA(ilutJ, tmp_J)
block
real(dp) :: discard
discard = 1.0_dp
if (t_hamil .or. (tFillingStochRDMOnFly .and. present(rdm_mat))) then
! call calc_mixed_start_contr_sym(tmp_I, csf_i, tmp_J, excitInfo, discard, &
! also_discard, integral, rdm_ind, rdm_mat)
call calc_mixed_start_contr_integral(tmp_I, csf_i, tmp_J, excitInfo, &
integral, rdm_ind, rdm_mat)
if (present(rdm_mat)) rdm_mat = rdm_mat * guga_mat
if (typ == excit_type%fullstart_L_to_R) then
mat_ele = guga_mat * ((get_umat_el(st, se, en, st) + &
get_umat_el(se, st, st, en)) / 2.0_dp + integral)
else if (typ == excit_type%fullstart_R_to_L) then
mat_ele = guga_mat * ((get_umat_el(st, en, se, st) + &
get_umat_el(en, st, st, se)) / 2.0_dp + integral)
end if
end if
end block
end associate
end subroutine calc_fullstart_mixed_ex