subroutine gen_double_back_spawn(nI, ilutI, part_type, nJ, ilutJ, ex, tPar, pgen)
! the double excitation routine for the back-spawn method
integer, intent(in) :: nI(nel)
integer(n_int), intent(in) :: ilutI(0:NIfTot)
integer, intent(in) :: part_type
integer, intent(out) :: nJ(nel), ex(2, maxExcit)
integer(n_int), intent(out) :: ilutJ(0:NIfTot)
logical, intent(out) :: tpar
real(dp), intent(out) :: pgen
character(*), parameter :: this_routine = "gen_double_back_spawn"
integer :: elecs(2), sym_product, ispn, sum_ml, src(2), loc, cc_a, cc_b, &
orbs(2)
real(dp) :: int_cpt(2), cum_sum(2), sum_pair(2), cpt_pair(2), &
cum_arr(nbasis)
if (t_back_spawn_flex) then
! Pick the electrons in a weighted fashion
call pick_weighted_elecs(nI, elecs, src, sym_product, ispn, sum_ml, &
pgen)
loc = check_electron_location(src, 2, part_type)
else
call pick_virtual_electrons_double(nI, part_type, elecs, src, ispn, &
sum_ml, pgen)
! here it could be that there are no valid electrons..
if (elecs(1) == 0) then
nJ(1) = 0
pgen = 0.0_dp
return
end if
! TODO!! thats the problem here! there are circular dependencies
! when I call that from here.. so i guess i need to setup the
! back_spawn excitation routines in a seperate file
sym_product = RandExcitSymLabelProd(SpinOrbSymLabel(src(1)), &
SpinOrbSymLabel(src(2)))
! indicate no flex option is used
loc = -1
end if
! at some point i have to make this whole logic easier here..
! in the end i think i have to remove some possibilities and stick to
! one back-spawn method.
if ((t_back_spawn_occ_virt) .or. (t_back_spawn_flex .and. ( &
(loc == 1 .and. occ_virt_level /= -1) .or. (loc == 2) .or. &
(loc == 0 .and. occ_virt_level >= 1)))) then
call pick_occupied_orbital(ilutI, src, ispn, part_type, int_cpt(1), cum_sum(1), &
orbs(1))
else
orbs(1) = pick_a_orb(ilutI, src, iSpn, int_cpt(1), cum_sum(1), cum_arr)
end if
if (orbs(1) /= 0) then
cc_a = ClasSCountInd(orbs(1))
cc_b = get_paired_cc_ind(cc_a, sym_product, sum_ml, iSpn)
if (t_back_spawn_flex .and. ((loc == 2 .and. occ_virt_level /= -1) .or. &
(occ_virt_level == 2))) then
call pick_second_occupied_orbital(ilutI, cc_b, orbs(1), ispn, &
part_type, int_cpt(2), cum_sum(2), orbs(2))
else
orbs(2) = select_orb(ilutI, src, cc_b, orbs(1), int_cpt(2), &
cum_sum(2))
end if
ASSERT((.not. (is_beta(orbs(2)) .and. .not. is_beta(orbs(1)))))
end if
if (any(orbs == 0)) then
nJ(1) = 0
pgen = 0.0_dp
return
end if
! can i exit right away if this happens??
! i am pretty sure this means
if (any(cum_sum < EPS)) then
cum_sum = 1.0_dp
int_cpt = 0.0_dp
if (.not. same_spin(orbs(1), orbs(2))) then
nJ(1) = 0
pgen = 0.0_dp
return
end if
end if
! only on parallel excitations.. and symmetric exciation generator is
! turned off for now in the back-spawning
if (same_spin(orbs(1), orbs(2))) then
if (t_back_spawn_occ_virt .or. (t_back_spawn_flex .and. ( &
(loc == 1 .and. (occ_virt_level == 0 .or. occ_virt_level == 1)) &
.or. (loc == 2 .and. occ_virt_level == -1) .or. &
(loc == 0 .and. occ_virt_level == 1)))) then
if (is_in_ref(orbs(2), part_type)) then
! if (b) is also in the occupied manifold i could have
! picked the other way around..
! with the same uniform probability:
cpt_pair(1) = int_cpt(1)
sum_pair(1) = cum_sum(1)
! and then (a) would have been picked according to the
! "normal" procedure
call pgen_select_orb(ilutI, src, orbs(2), orbs(1), &
cpt_pair(2), sum_pair(2))
else
! if (b) is not in the occupied this does not work or?
! since i am forcing (a) to be in the occupied..
! so remove this pgen:
cpt_pair = 0.0_dp
sum_pair = 1.0_dp
end if
else if (t_back_spawn_flex .and. ((loc == 0 .and. occ_virt_level == 2) &
.or. (loc == 1 .and. occ_virt_level == 2) .or. &
(loc == 2 .and. occ_virt_level /= -1))) then
! we have picked both in the occupied manifold
cpt_pair = int_cpt
sum_pair = cum_sum
else
! otherwise "normal"
call pgen_select_a_orb(ilutI, src, orbs(2), iSpn, cpt_pair(1), &
sum_pair(1), cum_arr, .false.)
call pgen_select_orb(ilutI, src, orbs(2), orbs(1), &
cpt_pair(2), sum_pair(2))
end if
else
cpt_pair = 0.0_dp
sum_Pair = 1.0_dp
end if
if (any(sum_pair < EPS)) then
cpt_pair = 0.0_dp
sum_pair = 1.0_dp
! can something like that slip through:
if (any(cum_sum < EPS)) then
pgen = 0.0_dp
nJ(1) = 0
return
end if
end if
pgen = pgen * (product(int_cpt) / product(cum_sum) + &
product(cpt_pair) / product(sum_pair))
! And generate the actual excitation
call make_double(nI, nJ, elecs(1), elecs(2), orbs(1), orbs(2), &
ex, tpar)
ilutJ = make_ilutJ(ilutI, ex, 2)
end subroutine gen_double_back_spawn