gen_double_back_spawn_ueg_new Subroutine

    subroutine gen_double_back_spawn_ueg_new(nI, ilutI, part_type, nJ, ilutJ, tPar, ex, pgen)
        integer, intent(in) :: nI(nel), part_type
        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) :: tPar
        real(dp), intent(out) :: pgen

        integer :: elecs(2), src(2), ispn, sum_ml, loc, orb_a, orb_b
        real(dp) :: p_elec, p_orb, dummy, cum_arr(nBasis), cum_sum

        logical :: t_temp_back_spawn

        t_temp_back_spawn = .false.

        if (t_back_spawn) then
            call pick_virtual_electrons_double(nI, part_type, elecs, src, ispn, &
                                               sum_ml, p_elec)

            loc = -1
        else
            call pick_uniform_elecs(elecs, p_elec)
            src = nI(elecs)

            ispn = get_ispn(src)

            loc = check_electron_location(src, 2, part_type)
        end if

        if (elecs(1) == 0) then
            nJ(1) = 0
            pgen = 0.0_dp
            return
        end if

        if ((loc == 2) .or. (loc == 1 .and. occ_virt_level /= -1) .or. &
            (loc == 0 .and. occ_virt_level >= 1)) then

            t_temp_back_spawn = .true.
            call pick_occupied_orbital_ueg(ilutI, src, iSpn, part_type, p_orb, &
                                           dummy, orb_a)
            ! it can happen that there are no valid orbitals
            if (orb_a == 0) then
                nJ(1) = 0
                pgen = 0.0_dp
                return
            end if

            ! i think in this case i have to multiply if both are in the
            ! reference or?

        else

            ! i could refactor that in a smaller function:
            call create_ab_list_ueg(ilutI, src, cum_arr, cum_sum)

            if (cum_sum < EPS) then
                nJ(1) = 0
                pgen = 0.0_dp
                return
            end if

            orb_a = binary_search_first_ge(cum_arr, genrand_real2_dsfmt() * cum_sum)

            if (orb_a == 1) then
                p_orb = cum_arr(1) / cum_sum
            else
                p_orb = (cum_arr(orb_a) - cum_arr(orb_a - 1)) / cum_sum
            end if

        end if

        orb_b = get_orb_from_kpoints(src(1), src(2), orb_a)

        ! thats it i guess..
        ! but i might have to be careful:
        ! in the back-spawn mechanism any order of the orbitals is possible..
        ! in the new ueg by NB not.. there orb_b > orb_a is strictly
        ! enforced.. hm.. how should we handle that?
        ! wait a minute.. since when is this called with the bare
        ! eleci and elecj?? yes it is apparently..
        call make_double(nI, nJ, elecs(1), elecs(2), orb_a, orb_b, ex, tpar)

        ilutJ = make_ilutJ(ilutI, ex, 2)

        pgen = p_elec * p_orb

        if (t_temp_back_spawn .and. is_in_ref(orb_b, part_type)) then
            pgen = 2.0_dp * pgen
        end if

    end subroutine gen_double_back_spawn_ueg_new