pure subroutine calc_fullstop_alike_ex(csf_i, csf_j, excitInfo, mat_ele, &
t_hamil, rdm_ind, rdm_mat)
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_hamil
integer(int_rdm), intent(out), allocatable, optional :: rdm_ind(:)
real(dp), intent(out), allocatable, optional :: rdm_mat(:)
character(*), parameter :: this_routine = "calc_fullstop_alike_ex"
real(dp) :: bVal, temp_mat, nOpen, guga_mat
HElement_t(dp) :: umat
integer :: i, step1, step2, db, delta_b(nSpatOrbs)
logical :: t_hamil_
def_default(t_hamil_, t_hamil, .true.)
ASSERT(present(rdm_ind) .eqv. present(rdm_mat))
! set some defaults in case of early exit
mat_ele = h_cast(0.0_dp)
delta_b = csf_i%B_int - csf_j%B_int
associate (ii => excitInfo%i, jj => excitInfo%j, kk => excitInfo%k, &
ll => excitInfo%l, typ => excitInfo%typ, st => excitInfo%fullStart, &
se => excitInfo%secondStart, gen => excitInfo%gen1, &
en => excitInfo%fullEnd)
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(1) = contract_2_rdm_ind(ii, jj, kk, ll)
end if
! can i exclude every deltaB > 1, since only db = 0 allowed in
! double overlap region? i think so..
if (any(abs(delta_b) > 1)) return
if (t_hamil_) then
if (typ == excit_type%fullstop_lowering) then
! LL
umat = (get_umat_el(en, en, st, se) + &
get_umat_el(en, en, se, st)) / 2.0_dp
else if (typ == excit_type%fullstop_raising) then
! RR
umat = (get_umat_el(st, se, en, en) + &
get_umat_el(se, st, en, en)) / 2.0_dp
end if
else
umat = h_cast(1.0_dp)
end if
if (t_hamil_ .and. near_zero(umat)) return
guga_mat = 1.0_dp
! have to deal with start specifically
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, gen, 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, gen, bval)
if (near_zero(guga_mat)) return
end do
! do semi-start:
step1 = csf_i%stepvector(se)
step2 = csf_j%stepvector(se)
db = delta_b(se - 1)
bVal = csf_i%b_real(se)
call getDoubleMatrixElement(step2, step1, db, gen, gen, bVal, 1.0_dp, temp_mat)
guga_mat = guga_mat * temp_mat
if (near_zero(guga_mat)) return
nOpen = (-1.0_dp)**real(count_open_orbs_ij(csf_i, se + 1, en - 1), dp)
! is this the same for both type of gens?
mat_ele = guga_mat * nOpen * Root2 * umat
! since the x1 element is 0 there is no sign influence from the
! order of generators!
if (present(rdm_mat)) rdm_mat = guga_mat * nOpen * Root2
end associate
end subroutine calc_fullstop_alike_ex