calc_pgen_back_spawn_ueg Function

    function calc_pgen_back_spawn_ueg(ilutI, ex, ic, part_type) result(pgen)
        integer, intent(in) :: ex(2, 2), ic, part_type
        integer(n_int), intent(in) :: ilutI(0:niftot)
        real(dp) :: pgen

        real(dp) :: cum_sum, pAIJ
        integer :: dummy_orb, ispn, loc, src(2), tgt(2)

        ! i have to write a generation probability calculator for the hphf
        ! implementation this should only be called if it is a non-init
        ! and back-spawn is on.. otherwise just use the already provided
        ! ones
        if (ic /= 2) then
            pgen = 0.0_dp
            return
        end if

        if (test_flag(ilutI, get_initiator_flag(part_type))) then

            ! i have to do some symmetry setup beforehand..
            ! or i do it by hand to avoid the unnecessary overhead..
            ! nope.. i actually only need:
            ! although i should not land here i guess..
            ! this functionality i could actually unit-test.. damn..
            call CalcPGenLattice(ex, pgen)
        else
            ! do the back-spawn pgen..
            ! elec were picked randomly(but in both orders!)
            ! and then we have to check the electron location, to know if
            ! there was a restriction on the first orbitals
            ! the electrons are in the first column or?
            src = get_src(ex)
            loc = check_electron_location(src, 2, part_type)

            ! and with the implementation right now it is only restricted
            ! if both were in the occup
            ! we extend this also in the ueg for occ_virt_level
            if ((loc == 2) .or. (loc == 1 .and. occ_virt_level /= -1) .or. &
                (loc == 0 .and. occ_virt_level >= 1)) then

                ! i could just call  the orbital picking routine..
                ! although this is inefficient..
                ! and i have to see if it would have been possible the
                ! other way around too..
                ! i need ispn here..
                ispn = get_ispn(src)
                call pick_occupied_orbital_ueg(ilutI, src, ispn, part_type, pAIJ, cum_sum, &
                                               dummy_orb, .true.)

                ! i think i can't just use 2*p(a|ij) since this assumption
                ! comes from p(b|ij) = p(a|ij) which is not true by
                ! default if we exclude a to the occupied manifold..
                ! i actually have to check if b is also in the
                ! occupied
                ! only mult by 2 if b is also in the occupied manifolf and
                ! thus it would have been possible to pick it the other
                ! way around
                pgen = 2.0_dp * pAIJ / real(nel * (nel - 1), dp)
                ! i should check the hole location too..
                tgt = get_tgt(ex)
                ! is tgt(1) always the first picked? yes it is!
                if (is_in_ref(tgt(1), part_type) .and. is_in_ref(tgt(2), part_type)) then
                    ! in this case we can recalc p(b|ij)
                    pgen = 2.0_dp * pgen
                end if

            else
                ! otherwise the orbital (a) is free
                ! so it is the number of available orbitals, depending on
                ! the spin
                ! here i can just use the formulas for the standard ueg..
                ! or just use the above provided routine
                call CalcPGenLattice(ex, pgen)
            end if
        end if

    end function calc_pgen_back_spawn_ueg