pure subroutine calc_normal_double_ex(csf_i, csf_j, excitInfo, mat_ele, &
t_hamil, rdm_ind, rdm_mat)
! combined routine to calculate the mixed generator excitations with
! 4 different spatial orbitals. here i have to consider if a
! spin change happened in the overlap. if no, i also have to calc. the
! non-overlap contribution!
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_normal_double_ex"
integer :: step1, step2, db, i, delta_b(nSpatOrbs)
real(dp) :: temp_x0, temp_x1, temp_mat0, temp_mat1, bVal, guga_mat
logical :: t_hamil_
def_default(t_hamil_, t_hamil, .true.)
ASSERT(present(rdm_ind) .eqv. present(rdm_mat))
! set defaults for early exits
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, start1 => excitInfo%fullStart, &
start2 => excitInfo%secondStart, ende1 => excitInfo%firstEnd, &
ende2 => excitInfo%fullEnd, gen1 => excitInfo%gen1, &
gen2 => excitInfo%gen2, firstgen => excitInfo%firstgen, &
lastgen => excitInfo%lastgen, order => excitInfo%order, &
order1 => excitInfo%order1, typ => excitInfo%typ)
if (present(rdm_ind) .and. present(rdm_mat)) then
allocate(rdm_ind(2), source=0_int_rdm)
allocate(rdm_mat(2), source=0.0_dp)
! this does get tricky now with the rdm inds and mats..
! the indices which end up here should always be intertwined..
! as we always deal with an overlap range in the excitation
! generation
ASSERT(max(ii, jj) > min(kk, ll))
! so the first two entries correspond to the overlap version
! of the generators (in the case of mixed R and L)
rdm_ind(1) = contract_2_rdm_ind(ii, jj, kk, ll)
! and the second to the non-overlap (again only in the case of
! mixed generator combinations!)
rdm_ind(2) = contract_2_rdm_ind(ii, ll, kk, jj)
end if
! i have to check if the deltaB value is in its correct bounds for
! the specific parts of the excitations
if (any(abs(delta_b(start1:start2 - 1)) > 1) .or. &
any(abs(delta_b(start2:ende1 - 1)) > 2) .or. &
any(abs(delta_b(ende1:ende2)) > 1)) return
! depending on which type of excitation the non-overlap has a specific
! tbd generator! i should deal with that in the starting if statement
! according to my stochastic implementation it is not so hard luckily..
! do first single part
guga_mat = 1.0_dp
! have to do the start specifically as i need the outgoing deltaB
! value
step1 = csf_i%stepvector(start1)
step2 = csf_j%stepvector(start1)
db = delta_b(start1)
bVal = csf_i%b_real(start1)
guga_mat = guga_mat * getSingleMatrixElement(step2, step1, db, firstgen, bVal)
if (near_zero(guga_mat)) return
do i = start1 + 1, start2 - 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, firstgen, bVal)
if (near_zero(guga_mat)) return
end do
temp_x0 = guga_mat
temp_x1 = guga_mat
! do overlap part
do i = start2, ende1
step1 = csf_i%stepvector(i)
step2 = csf_j%stepvector(i)
bVal = csf_i%b_real(i)
db = delta_b(i - 1)
call getDoubleMatrixElement(step2, step1, db, gen1, gen2, 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
! now do second single overlap part
do i = ende1 + 1, ende2
step1 = csf_i%stepvector(i)
step2 = csf_j%stepvector(i)
bVal = csf_i%b_real(i)
db = delta_b(i - 1)
temp_mat0 = getSingleMatrixElement(step2, step1, db, lastgen, bVal)
temp_x0 = temp_x0 * temp_mat0
temp_x1 = temp_x1 * temp_mat0
if (near_zero(temp_x0) .and. near_zero(temp_x1)) return
end do
! now the excitation-type and spin question come into play..
! i actually could combine all "normal" double excitations in one
! routine with a if statement at the end here..
if (present(rdm_mat)) then
select case (typ)
case (excit_type%double_lowering, &
excit_type%double_raising)
! in both the orbital picker and also excitation identifier
! the order quantities are setup to be 1.0..
! i hope I did everything right there.. then this would
! mean here
ASSERT(order.isclose.1.0_dp)
ASSERT(order1.isclose.1.0_dp)
! and now I have to fill in the correct combinations of
! signs here..
! to check that everything is correctly set to the default
! ordering I should assert that here
#ifdef DEBUG_
if (typ == excit_type%double_raising) then
ASSERT(ii < jj)
ASSERT(kk < ll)
ASSERT(ii > kk)
ASSERT(jj > ll)
end if
#endif
! be sure with the rdm_sign function:
! rdm_mat(1) = temp_x0 + generator_sign(ii,jj,kk,ll) * temp_x1
! rdm_mat(2) = temp_x0 + generator_sign(ii,ll,kk,jj) * temp_x1
rdm_mat(1) = temp_x0 + temp_x1
rdm_mat(2) = temp_x0 - temp_x1
case (excit_type%double_L_to_R_to_L, &
excit_type%double_R_to_L_to_R, &
excit_type%double_L_to_R, &
excit_type%double_R_to_L)
if (excitInfo%spin_change) then
! if we have a spin-change the non-overlap contribution
! must be 0! which is already intiated to 0 above
rdm_mat(1) = temp_x1
else
! if there is not spin-change the non-overlap is also
! 0, and in this case is -2 the original x0
rdm_mat(1) = temp_x0 + temp_x1
rdm_mat(2) = -2.0_dp * temp_x0
end if
end select
end if
select case (typ)
case (excit_type%double_lowering)
! double lowering
! wait a minute.. i have to do that at the end apparently..
! since i need to know the x0 and x1 matrix element contributions
mat_ele = (temp_x0 * (get_umat_el(ende1, ende2, start1, start2) + &
get_umat_el(ende2, ende1, start2, start1) + &
get_umat_el(ende2, ende1, start1, start2) + &
get_umat_el(ende1, ende2, start2, start1)) + &
order * order1 * temp_x1 * ( &
get_umat_el(ende1, ende2, start1, start2) + &
get_umat_el(ende2, ende1, start2, start1) - &
get_umat_el(ende2, ende1, start1, start2) - &
get_umat_el(ende1, ende2, start2, start1))) / 2.0_dp
case (excit_type%double_raising)
! double raising
mat_ele = (temp_x0 * (get_umat_el(start1, start2, ende1, ende2) + &
get_umat_el(start2, start1, ende2, ende1) + &
get_umat_el(start1, start2, ende2, ende1) + &
get_umat_el(start2, start1, ende1, ende2)) + &
order * order1 * temp_x1 * ( &
get_umat_el(start1, start2, ende1, ende2) + &
get_umat_el(start2, start1, ende2, ende1) - &
get_umat_el(start1, start2, ende2, ende1) - &
get_umat_el(start2, start1, ende1, ende2))) / 2.0_dp
case (excit_type%double_L_to_R_to_L)
! L -> R -> L
if (excitInfo%spin_change) then
! if a spin-change happenend -> no non-overlap!
mat_ele = temp_x1 * (get_umat_el(ende2, start2, start1, ende1) + &
get_umat_el(start2, ende2, ende1, start1)) / 2.0_dp
else
mat_ele = (-temp_x0 * (get_umat_el(start2, ende2, start1, ende1) + &
get_umat_el(ende2, start2, ende1, start1)) * 2.0_dp + &
(temp_x0 + temp_x1) * (get_umat_el(ende2, start2, start1, ende1) + &
get_umat_el(start2, ende2, ende1, start1))) / 2.0_dp
end if
case (excit_type%double_R_to_L_to_R)
! R -> L -> R
if (excitInfo%spin_change) then
mat_ele = temp_x1 * (get_umat_el(start1, ende1, ende2, start2) + &
get_umat_el(ende1, start1, start2, ende2)) / 2.0_dp
else
mat_ele = (-temp_x0 * (get_umat_el(start1, ende1, start2, ende2) + &
get_umat_el(ende1, start1, ende2, start2)) * 2.0_dp + &
(temp_x0 + temp_x1) * (get_umat_el(start1, ende1, ende2, start2) + &
get_umat_el(ende1, start1, start2, ende2))) / 2.0_dp
end if
case (excit_type%double_L_to_R)
! L -> R
if (excitInfo%spin_change) then
mat_ele = temp_x1 * (get_umat_el(ende1, start2, start1, ende2) + &
get_umat_el(start2, ende1, ende2, start1)) / 2.0_dp
else
mat_ele = (-temp_x0 * (get_umat_el(start2, ende1, start1, ende2) + &
get_umat_el(ende1, start2, ende2, start1)) * 2.0_dp + &
(temp_x0 + temp_x1) * (get_umat_el(ende1, start2, start1, ende2) + &
get_umat_el(start2, ende1, ende2, start1))) / 2.0_dp
end if
case (excit_type%double_R_to_L)
! R -> L
if (excitInfo%spin_change) then
mat_ele = temp_x1 * (get_umat_el(start1, ende2, ende1, start2) + &
get_umat_el(ende2, start1, start2, ende1)) / 2.0_dp
else
mat_ele = (-temp_x0 * (get_umat_el(start1, ende2, start2, ende1) + &
get_umat_el(ende2, start1, ende1, start2)) * 2.0_dp + &
(temp_x0 + temp_x1) * (get_umat_el(start1, ende2, ende1, start2) + &
get_umat_el(ende2, start1, start2, ende1))) / 2.0_dp
end if
! combine the "normal" double RR/LL also in here, since the rest
! of the routine is totally the same!
end select
end associate
end subroutine calc_normal_double_ex