calc_pgen_back_spawn_hubbard Function

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

        integer :: d_elecs(2), d_src(2), d_ispn, src(2), loc, tgt(2), d_orb
        real(dp) :: pgen_elec, paij, mult

        ! the hubbard pgen recalculation is pretty similar to the ueg one
        ! below..
        if (ic /= 2) then
            pgen = 0.0_dp
            return
        end if

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

            ! can i use the already provided routine also for hubbard models?
            ! yes!
            call CalcPGenLattice(ex, pgen)
        else
            ! in the hubbard case it can also be that we pick the electrons
            ! with the old back-spawn method
            src = get_src(ex)

            if (t_back_spawn) then

                call pick_virtual_electrons_double_hubbard(nI, part_type, d_elecs, d_src, &
                                                           d_ispn, pgen_elec, .true.)

                loc = -1
            else

                ! otherwise:
                pgen_elec = 1.0_dp / real(nOccBeta * nOccAlpha, dp)

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

            ! otherwise i have to calculate stuff
            if ((loc == 2) .or. (loc == 1 .and. occ_virt_level /= -1) .or. &
                (loc == 0 .and. occ_virt_level >= 1)) then
                ! i only do that in the back-spawn flex case..

                call pick_occupied_orbital_hubbard(ilutI, part_type, paij, &
                                                   d_orb, .true.)

                tgt = get_tgt(ex)

                ! ok i just realized.. the excitation matrix gets always
                ! sorted.. thats not good for my implementation..
                ! atleast that means that i cant assume that the second
                ! picked orbital is always at position tgt(2)
                ! is it enough to check if both are in the reference?
                ! because when i am here, i know that atleast one is
                ! definetly in the reference.. and if the second also is
                ! i could have picked the orbitals the other way around..
                if (is_in_ref(tgt(1), part_type) .and. is_in_ref(tgt(2), part_type)) then
                    mult = 2.0_dp
                else
                    mult = 1.0_dp
                end if

            else
                ! otherwise we can pick the orbital (a) freely
                paij = 1.0_dp / real(nbasis - nel, dp)
                mult = 2.0_dp

            end if

            pgen = mult * paij * pgen_elec

        end if

    end function calc_pgen_back_spawn_hubbard