pure subroutine calc_fullstop_mixed_ex(ilutI, csf_i, ilutJ, csf_j, excitInfo, mat_ele, &
t_hamil, rdm_ind, rdm_mat)
! from the excitInfo i know the first switch position.
! this makes things a bit easier for the exact calculation
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 :: i, step1, step2, db
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)
integer :: st, se, en, delta_b(nSpatOrbs)
real(dp) :: guga_mat
logical :: t_hamil_
def_default(t_hamil_, t_hamil, .true.)
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
se = excitInfo%secondStart
en = excitInfo%fullEnd
associate (ii => excitInfo%i, jj => excitInfo%j, kk => excitInfo%k, &
ll => excitInfo%l, firstGen => excitInfo%firstgen, &
typ => excitInfo%typ)
! set defaults in case of early exit
mat_ele = h_cast(0.0_dp)
if (any(abs(delta_b(st:se - 1)) > 1) .or. &
any(abs(delta_b(se:en)) > 2)) return
! first do single overlap region
guga_mat = 1.0_dp
step1 = csf_i%stepvector(st)
step2 = csf_j%stepvector(st)
bVal = csf_i%b_real(st)
db = delta_b(st)
guga_mat = guga_mat * getSingleMatrixElement(step2, step1, db, firstgen, bval)
if (near_zero(guga_mat)) return
do i = st + 1, se - 1
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, firstgen, bVal)
if (near_zero(guga_mat)) return
end do
! i also do not have to consider if there is a d=3 at the end since
! i determine the last spin-coupling change and thus know it is
! definetly a singly occupied orbital at the end
! and actually the x0 matrix element has to be 0, otherwise it is
! not the excitation i thought it was.. do this as a check to
! abort if anything else happens
temp_x0 = guga_mat
temp_x1 = guga_mat
do i = se, en - 1
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%R, gen_type%L, bVal, &
1.0_dp, temp_mat0, temp_mat1)
temp_x0 = temp_x0 * temp_mat0
temp_x1 = temp_x1 * temp_mat1
if (near_zero(temp_x0) .and. near_zero(temp_x1)) return
end do
! to the fullstop
step1 = csf_i%stepvector(en)
step2 = csf_j%stepvector(en)
db = delta_b(en - 1)
bVal = csf_i%b_real(en)
call getMixedFullStop(step2, step1, db, bVal, temp_mat0, temp_mat1)
temp_x0 = temp_x0 * temp_mat0
temp_x1 = temp_x1 * temp_mat1
! so if x0 > 0 abort!
! and here misuse the stochastic routine to calculate the influence
! of all the other singly-occupied orbitals..
! i probably should write a function which does that only for the
! integral and stores the specific indices and matrix elements for
! the rdm calculation, but to that later! todo
! but to use this function i have to transform the 2 iluts to
! GUGA iluts and store the matrix element in them..
call convert_ilut_toGUGA(ilutI, tmp_I)
call convert_ilut_toGUGA(ilutJ, tmp_J)
! is the matrix element transformed within these functions?
! no apparently not and not event touched.. so no need to encode the
! matrix element
block
real(dp) :: discard
discard = 1.0_dp
if (t_hamil_ .or. (tFillingStochRDMOnFly .and. present(rdm_mat))) then
if (typ == excit_type%fullstop_L_to_R) then
! L -> R
! what do i have to put in as the branch pgen?? does it have
! an influence on the integral and matrix element calculation?
! call calc_mixed_end_contr_sym(tmp_I, csf_i, tmp_J, excitInfo, discard, &
! also_discard, integral, rdm_ind, rdm_mat)
call calc_mixed_end_contr_integral(tmp_I, csf_i, tmp_J, &
excitInfo, integral, rdm_ind, rdm_mat)
! need to multiply by x1
if (present(rdm_mat)) rdm_mat = rdm_mat * temp_x1
mat_ele = temp_x1 * ((get_umat_el(en, se, st, en) + &
get_umat_el(se, en, en, st)) / 2.0_dp + integral)
else if (typ == excit_type%fullstop_R_to_L) then
! R -> L
! call calc_mixed_end_contr_sym(tmp_I, csf_i, tmp_J, excitInfo, discard, &
! also_discard, integral, rdm_ind, rdm_mat)
call calc_mixed_end_contr_integral(tmp_I, csf_i, tmp_J, &
excitInfo, integral, rdm_ind, rdm_mat)
if (present(rdm_mat)) rdm_mat = rdm_mat * temp_x1
mat_ele = temp_x1 * ((get_umat_el(en, st, se, en) + &
get_umat_el(st, en, en, se)) / 2.0_dp + integral)
end if
else
mat_ele = h_cast(temp_x1)
end if
end block
end associate
end subroutine calc_fullstop_mixed_ex