subroutine gen_double_4ind_ex2(nI, ilutI, nJ, ilutJ, ex, par, pgen)
use GenRandSymExcitNUMod, only: RandExcitSymLabelProd
use SymExcitDataMod, only: SpinOrbSymLabel
use lattice_models_utils, only: pick_spin_opp_elecs
integer, intent(in) :: nI(nel)
integer(n_int), intent(in) :: ilutI(0:NIfTot)
integer, intent(out) :: nJ(nel), ex(2, maxExcit)
integer(n_int), intent(out) :: ilutJ(0:NIfTot)
logical, intent(out) :: par
real(dp), intent(out) :: pgen
character(*), parameter :: this_routine = 'gen_double_4ind_ex2'
integer :: elecs(2), src(2), orbs(2), ispn, sum_ml, cc_a, cc_b
integer :: sym_product
real(dp) :: int_cpt(2), cum_sum(2), cpt_pair(2), sum_pair(2)
! This array stores the cumulative list used in selecting the first
! orbital so that it can be used for calculating the reverse
! probability if required.
real(dp) :: cum_arr(nbasis)
real(dp) :: scratch_cpt, scratch_sm
integer :: scratch_orb
integer :: loc
! Pick the electrons in a weighted fashion
if (t_mixed_hubbard .or. t_olle_hubbard) then
call pick_biased_elecs(nI, elecs, src, sym_product, ispn, sum_ml, pgen)
else
call pick_weighted_elecs(nI, elecs, src, sym_product, ispn, sum_ml, &
pgen)
end if
! then first pick (a) orbital:
! for opposite spin excitations (a) is restricted to be a beta orbital!
! and the probability is split p(a|ij) = p(j)*p(a|i)
! except in the symmetric excitation generator, which isn't used
! ever anyway..
orbs(1) = pick_a_orb(ilutI, src, iSpn, int_cpt(1), cum_sum(1), cum_arr)
! Select the B orbital, in the same way as before!!
! The symmetry of this second orbital depends on that of the first.
if (orbs(1) /= 0) then
cc_a = ClassCountInd(orbs(1))
cc_b = get_paired_cc_ind(cc_a, sym_product, sum_ml, iSpn)
! pick the last orbitals weighted with the exact matrix
! element
orbs(2) = select_orb(ilutI, src, cc_b, orbs(1), int_cpt(2), &
cum_sum(2))
end if
! what does this assert do? do i have to pick the electrons in a
! certain order??
ASSERT((.not. (is_beta(orbs(2)) .and. .not. is_beta(orbs(1)))) .or. tGen_4ind_2_symmetric)
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
end if
! Calculate the pgens. Note that all of these excitations can be
! selected as both A--B or B--A. So these need to be calculated
! explicitly.
if (.not. tGen_guga_crude) then
ASSERT(tGen_4ind_part_exact)
end if
if ((is_beta(orbs(1)) .eqv. is_beta(orbs(2))) .or. tGen_4ind_2_symmetric) then
! in the case of parallel spin excitations or symmetrice excitation
! generation(but does actually someone use that?) we have to
! calculate the probability of picking the holes the other
! way around
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))
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
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, par)
ilutJ = ilutI
clr_orb(ilutJ, src(1))
clr_orb(ilutJ, src(2))
set_orb(ilutJ, orbs(1))
set_orb(ilutJ, orbs(2))
end subroutine gen_double_4ind_ex2