subroutine gen_ab_cum_list_3_3(csf_i, occ_orbs, cum_arr, excit_arr, orb_arr)
! specific routine when the occupaton of the already picked orbitals
! is 2 -> still think if i really want to outpout a cum_arr of lenght
! nBasis -> nSpatOrbs would be better and doable... !! todo
type(CSF_Info_t), intent(in) :: csf_i
integer, intent(in) :: occ_orbs(2)
real(dp), intent(out) :: cum_arr(nSpatOrbs)
integer, intent(out) :: orb_arr(nSpatOrbs)
type(ExcitationInformation_t), intent(out) :: excit_arr(nSpatOrbs)
character(*), parameter :: this_routine = "gen_ab_cum_list_3_3"
integer :: ki(3), kj(3), n_id(2), orb_a, ka(3), kb(3), orb_b, st, en
integer :: z_ind
real(dp) :: cum_sum, contrib
type(ExcitationInformation_t) :: excitInfo
! for legacy compatibility
z_ind = ubound(kPointToBasisFn, 3)
ki = G1(occ_orbs(1))%k
kj = G1(occ_orbs(2))%k
cum_sum = 0.0_dp
orb_arr = 0
n_id = gtID(occ_orbs)
! due to k-point symmetry alot of excitation types:
! eg. if i and j are seperate a and b also have to be!
! not true! ki + kj = 2*ka can defo be!
do orb_a = 1, n_id(1) - 1
contrib = 0.0_dp
excitInfo%valid = .false.
! avoid doubly occupied orbitals
if (csf_i%stepvector(orb_a) /= 3) then
ka = G1(2 * orb_a)%k
kb = ki + kj - ka
call MomPbcSym(kb, nBasisMax)
! is b allowed by the size of space? -> TODO: what does that mean?
orb_b = gtID(KPointToBasisFn(kb(1), kb(2), z_ind, 1))
! just check if k-point restriction works as i thought
! i think this below still holds..
ASSERT(orb_b /= n_id(1) .and. orb_b /= n_id(2))
! just to be safe also check if (a = b)
if (orb_a == orb_b) then
if (csf_i%stepvector(orb_a) == 0) then
contrib = 2.0_dp
! _RR_(ab) > ^RR(i) > ^R(j)
excitInfo = assign_excitInfo_values_double( &
excit_type%fullstart_raising, &
gen_type%R, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
orb_a, n_id(1), orb_a, n_id(2), orb_a, orb_a, n_id(1), n_id(2), &
0, 2, 1.0_dp, 1.0_dp)
end if
else
! check guga restrictions todo
if (csf_i%stepvector(orb_b) /= 3) then
! no additional restrictions or?...
contrib = 1.0_dp
if (orb_b > n_id(2)) then
! _R(a) > _LR(i) > ^RL(j) > ^L(b)
excitInfo = assign_excitInfo_values_double( &
excit_type%double_R_to_L, &
gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%L, &
orb_a, n_id(2), orb_b, n_id(1), orb_a, n_id(1), n_id(2), orb_b, &
0, 4, 1.0_dp, 1.0_dp)
else if (orb_b < n_id(1)) then
! check if a == b
! check maxima
st = min(orb_a, orb_b)
en = max(orb_a, orb_b)
! _R(min) > _RR(max) > ^RR(i) > ^R(j)
excitInfo = assign_excitInfo_values_double( &
excit_type%double_raising, &
gen_type%R, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
en, n_id(2), st, n_id(1), st, en, n_id(1), n_id(2), &
0, 4, 1.0_dp, 1.0_dp)
else
! b is between i and j
! _R(a) > _LR(i) > ^LR(b) > ^R(j)
excitInfo = assign_excitInfo_values_double( &
excit_type%double_R_to_L_to_R, &
gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
orb_a, n_id(2), orb_b, n_id(1), orb_a, n_id(1), orb_b, n_id(2), &
0, 4, 1.0_dp, 1.0_dp)
end if
end if
end if
end if
cum_sum = cum_sum + contrib
cum_arr(orb_a) = cum_sum
excit_arr(orb_a) = excitInfo
orb_arr(orb_a) = orb_b
end do
if (n_id(1) == 1) then
cum_arr(1) = 0.0_dp
else
cum_arr(n_id(1)) = cum_arr(n_id(1) - 1)
end if
do orb_a = n_id(1) + 1, n_id(2) - 1
contrib = 0.0_dp
excitInfo%valid = .false.
! avoid doubly occupied orbitals
if (csf_i%stepvector(orb_a) /= 3) then
ka = G1(2 * orb_a)%k
kb = ki + kj - ka
call MomPbcSym(kb, nBasisMax)
! is b allowed by the size of space? -> TODO: what does that mean?
orb_b = gtID(KPointToBasisFn(kb(1), kb(2), z_ind, 1))
ASSERT(orb_b /= n_id(1) .and. orb_b /= n_id(2))
if (orb_a == orb_b) then
if (csf_i%stepvector(orb_a) == 0) then
contrib = 2.0_dp
! _L(i) > ^LR_(ab) > ^R(j)
excitInfo = assign_excitInfo_values_double( &
excit_type%single_overlap_L_to_R, &
gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%R, &
orb_a, n_id(1), orb_a, n_id(2), n_id(1), orb_a, orb_a, n_id(2), &
0, 2, 1.0_dp, 1.0_dp, 1)
end if
else
! check guga restrictions todo
if (csf_i%stepvector(orb_b) /= 3) then
! no additional restrictions or?...
contrib = 1.0_dp
if (orb_b < n_id(1)) then
! _R(b) > _LR(i) > ^LR(a) > ^R(j)
excitInfo = assign_excitInfo_values_double( &
excit_type%double_R_to_L_to_R, &
gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
orb_b, n_id(2), orb_a, n_id(1), orb_b, n_id(1), orb_a, n_id(2), &
0, 4, 1.0_dp, 1.0_dp)
else if (orb_b > n_id(2)) then
! _L(i) > _RL(a) > ^RL(j) > ^L(b)
excitInfo = assign_excitInfo_values_double( &
excit_type%double_L_to_R_to_L, &
gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%L, &
orb_b, n_id(1), orb_a, n_id(2), n_id(1), orb_a, n_id(2), orb_b, &
0, 4, 1.0_dp, 1.0_dp)
else
st = min(orb_a, orb_b)
en = max(orb_a, orb_b)
! _L(i) > _RL(min) > ^LR(max) > ^R(j)
excitInfo = assign_excitInfo_values_double( &
excit_type%double_L_to_R, &
gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%R, &
en, n_id(1), st, n_id(2), n_id(1), st, en, n_id(2), &
0, 4, 1.0_dp, 1.0_dp)
end if
end if
end if
end if
cum_sum = cum_sum + contrib
cum_arr(orb_a) = cum_sum
excit_arr(orb_a) = excitInfo
orb_arr(orb_a) = orb_b
end do
! fill in orbital j
! not so sure about that anymore.. i think i got that wrong and
! this case is possible! damn
cum_arr(n_id(2)) = cum_arr(n_id(2) - 1)
do orb_a = n_id(2) + 1, nSpatOrbs
contrib = 0.0_dp
excitInfo%valid = .false.
if (csf_i%stepvector(orb_a) /= 3) then
ka = G1(2 * orb_a)%k
kb = ki + kj - ka
call MomPbcSym(kb, nBasisMax)
! there is no "spin" restriction in the guga case
! so both possible b orbs have to be checked
! its actually NOT possible that ka = ki !! so do not have
! to check that case!
orb_b = gtID(KPointToBasisFn(kb(1), kb(2), z_ind, 1))
! this should work as we pick beta orbital first
ASSERT(orb_b /= n_id(1) .and. orb_b /= n_id(2))
if (orb_a == orb_b) then
if (csf_i%stepvector(orb_a) == 0) then
contrib = 2.0_dp
! _L(i) > _LL(j) > ^LL^(ab)
excitInfo = assign_excitInfo_values_double( &
excit_type%fullstop_lowering, &
gen_type%L, gen_type%L, gen_type%L, gen_type%L, gen_type%L, &
orb_a, n_id(1), orb_b, n_id(2), n_id(1), n_id(2), orb_a, orb_a, &
0, 2, 1.0_dp, 1.0_dp)
end if
else
if (csf_i%stepvector(orb_b) /= 3) then
! only then its a possible excitation
contrib = 1.0_dp
! check type of excitation
if (orb_b < n_id(1)) then
! _R(b) > _LR(i) > ^RL(j) > ^L(a)
excitInfo = assign_excitInfo_values_double( &
excit_type%double_R_to_L, &
gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%L, &
orb_b, n_id(2), orb_a, n_id(1), orb_b, n_id(1), n_id(2), orb_a, &
0, 4, 1.0_dp, 1.0_dp)
else if (orb_b > n_id(2)) then
st = min(orb_a, orb_b)
en = max(orb_a, orb_b)
! _L(i) > _LL(j) > ^LL(min) > ^L(max)
excitInfo = assign_excitInfo_values_double( &
excit_type%double_lowering, &
gen_type%L, gen_type%L, gen_type%L, gen_type%L, gen_type%L, &
st, n_id(1), en, n_id(2), n_id(1), n_id(2), st, en, &
0, 4, 1.0_dp, 1.0_dp)
else
! _L(i) > _RL(b) > ^RL(j) > ^L(a)
excitInfo = assign_excitInfo_values_double( &
excit_type%double_L_to_R_to_L, &
gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%L, &
orb_a, n_id(1), orb_b, n_id(2), n_id(1), orb_b, n_id(2), orb_a, &
0, 4, 1.0_dp, 1.0_dp)
end if
end if
end if
end if
cum_sum = cum_sum + contrib
cum_arr(orb_a) = cum_sum
excit_arr(orb_a) = excitInfo
orb_arr(orb_a) = orb_b
end do
! make orbital (j) invalid for excitation also... not actually
! necessary, as i never pick that orbital anyway,,
excitInfo%valid = .false.
excit_arr(n_id(2)) = excitInfo
end subroutine gen_ab_cum_list_3_3