pure subroutine calc_single_excitation_ex(csf_i, csf_j, excitInfo, mat_ele, &
t_calc_full, rdm_ind, rdm_mat)
! routine to exactly calculate the matrix element between so singly
! connected CSFs, with the option to output also all the indices and
! overlap matrix elements necessary for the rdm calculation
type(CSF_Info_t), intent(in) :: csf_i, csf_j
type(ExcitationInformation_t), intent(in) :: excitInfo
HElement_t(dp), intent(out) :: mat_ele
logical, intent(in), optional :: t_calc_full
integer(int_rdm), intent(out), allocatable, optional :: rdm_ind(:)
real(dp), intent(out), allocatable, optional :: rdm_mat(:)
character(*), parameter :: this_routine = "calc_single_excitation_ex"
integer :: iOrb, db, step1, step2
real(dp) :: bVal
HElement_t(dp) :: integral
logical :: t_calc_full_
! just mimick the stochastic single excitation!
real(dp) :: tmp_mat
integer :: delta_b(nSpatOrbs)
def_default(t_calc_full_, t_calc_full, .true.)
ASSERT(present(rdm_ind) .eqv. present(rdm_mat))
! set defaults for the output if we excit early..
mat_ele = h_cast(0.0_dp)
delta_b = csf_i%B_int - csf_j%B_int
associate (i => excitInfo%i, j => excitInfo%j, st => excitInfo%fullstart, &
en => excitInfo%fullEnd, gen => excitInfo%currentGen)
if (present(rdm_ind) .and. present(rdm_mat)) then
allocate(rdm_ind(1), source=0_int_rdm)
allocate(rdm_mat(1), source=0.0_dp)
rdm_ind = contract_1_rdm_ind(i, j)
rdm_mat = 0.0_dp
end if
! this deltaB info can slip through the excitation identifier..
if (any(abs(delta_b) > 1)) return
if (t_calc_full_) then
integral = getTmatEl(2 * i, 2 * j)
else
integral = h_cast(1.0_dp)
end if
tmp_mat = 1.0_dp
! what do i need from the 2 CSFs to calculate all the matrix elements?
! the 2 stepvalues: d, d' -> write a routine which only calculates those
! the generator type: gen
! the currentB value of I, b -> also calc. that in the routine to get d
! and the deltaB value: db
! since i mostly use this function to calculate the overlap with
! the reference determinant it is probably a good idea to
! calculate the necessary information for the reference determinant
! once and store it.. and maybe use a flag here to check what
! kind of usage of this function is.. or outside in the calling
! function..
! then i can just loop over the whole excitation range without
! distinction between start.. (mabye end i have to deal with specifically!)
! here i have to calculate the starting element
!
step1 = csf_i%stepvector(st)
step2 = csf_j%stepvector(st)
bVal = csf_i%b_real(st)
! i think i have to take the deltaB for the next orb or?
! because i need the outgoing deltaB value..
! nah.. i need the incoming!
! damn i think i need the deltaB value from the previous orbital..
! no.. for the single start i need the outgoing deltab, that the
! convention how to access the matrix elements, but then i need the
! incoming deltaB value.. or lets check that out how this works
! exactly
db = delta_b(st)
tmp_mat = tmp_mat * getSingleMatrixElement(step2, step1, db, gen, bVal)
! i think it can still always happen that the matrix element is 0..
! but maybe for the rdm case i have to do something more involved..
! to set all the corresponding indices to 0..
if (near_zero(tmp_mat)) return
do iOrb = st + 1, en - 1
step1 = csf_i%stepvector(iOrb)
step2 = csf_j%stepvector(iOrb)
bVal = csf_i%b_real(iOrb)
db = delta_b(iOrb - 1)
tmp_mat = tmp_mat * getSingleMatrixElement(step2, step1, db, gen, bVal)
if (near_zero(tmp_mat)) return
if (t_calc_full_) then
if (.not. (t_new_real_space_hubbard .or. t_heisenberg_model &
.or. t_tJ_model .or. t_mixed_hubbard)) then
integral = integral + get_umat_el(i, iOrb, j, iOrb) * csf_i%Occ_real(iOrb)
integral = integral + get_umat_el(i, iOrb, iOrb, j) * &
getDoubleContribution(step2, step1, db, gen, bVal)
end if
end if
end do
step1 = csf_i%stepvector(en)
step2 = csf_j%stepvector(en)
bVal = csf_i%b_real(en)
db = delta_b(en - 1)
tmp_mat = tmp_mat * getSingleMatrixElement(step2, step1, db, gen, bVal)
if (near_zero(tmp_mat)) return
! i think i could also exclude the treal case here.. try!
if (t_calc_full_) then
if (.not. (treal .or. t_new_real_space_hubbard .or. &
t_heisenberg_model .or. t_tJ_model .or. t_mixed_hubbard)) then
call calc_integral_contribution_single(csf_i, csf_j, i, j, st, en, integral)
end if
end if
mat_ele = tmp_mat * integral
if (present(rdm_mat)) rdm_mat = tmp_mat
end associate
end subroutine calc_single_excitation_ex