pure subroutine calc_single_overlap_mixed_ex(csf_i, csf_j, excitInfo, mat_ele, &
t_calc_full, rdm_ind, rdm_mat)
! routine to exactly calculate the matrix element between 2 CSFs
! connected by a single overlap excitation with mixed generators
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_overlap_mixed_ex"
real(dp) :: bVal, temp_mat, guga_mat
HElement_t(dp) :: umat
integer :: i, db, gen2, step1, step2, delta_b(nSpatOrbs)
logical :: t_calc_full_
! in the case of rdm calculation, i know that this type of exitation
! only has one (or two, with switches 2-body integrals..??)
! rdm-contribution..
ASSERT(present(rdm_ind) .eqv. present(rdm_mat))
def_default(t_calc_full_, t_calc_full, .true.)
! 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, st => excitInfo%fullStart, &
ss => excitInfo%secondStart, en => excitInfo%fullEnd, &
gen => excitInfo%firstGen, fe => excitInfo%firstEnd, &
typ => excitInfo%typ)
if (present(rdm_ind) .and. present(rdm_mat)) then
! i am not sure yet if I will use symmetries in the RDM
! calculation (some are also left out in the SD based implo..
! so for now sample both combinations
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
if (any(abs(delta_b) > 1)) return
if (t_calc_full_) then
if (typ == excit_type%single_overlap_L_to_R) then
umat = (get_umat_el(fe, ss, st, en) + &
get_umat_el(ss, fe, en, st)) / 2.0_dp
else if (typ == excit_type%single_overlap_R_to_L) then
umat = (get_umat_el(st, en, fe, ss) + &
get_umat_el(en, st, ss, fe)) / 2.0_dp
else
call stop_all(this_routine, "shouldnt be here!")
end if
else
umat = h_cast(1.0_dp)
end if
! for the hamiltonian matrix element i can exit here if umat is 0
! but for the rdm-contribution i need to calc. the GUGA element
if (t_calc_full .and. near_zero(umat)) return
guga_mat = 1.0_dp
! i have to do the start specifically, due to the access of the
! single matrix elements
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, ss - 1
step1 = csf_i%stepvector(i)
step2 = csf_j%stepvector(i)
bVal = csf_i%b_real(i)
db = delta_b(i - 1)
guga_mat = guga_mat * getSingleMatrixElement(step2, step1, db, gen, bVal)
if (near_zero(guga_mat)) return
end do
step1 = csf_i%stepvector(ss)
step2 = csf_j%stepvector(ss)
bVal = csf_i%b_real(ss)
db = delta_b(ss - 1)
gen2 = -gen
call getDoubleMatrixElement(step2, step1, db, gen, gen2, bVal, 1.0_dp, temp_mat)
guga_mat = guga_mat * temp_mat
if (near_zero(guga_mat)) return
do i = ss + 1, en
step1 = csf_i%stepvector(i)
step2 = csf_j%stepvector(i)
bVal = csf_i%b_real(i)
db = delta_b(i - 1)
guga_mat = guga_mat * getSingleMatrixElement(step2, step1, db, gen2, bVal)
if (near_zero(guga_mat)) return
end do
end associate
mat_ele = guga_mat * umat
! both coupling coeffs are the same
if (present(rdm_mat)) rdm_mat = guga_mat
end subroutine calc_single_overlap_mixed_ex