pick_ab_orbitals_par_hubbard Subroutine

    subroutine pick_ab_orbitals_par_hubbard(nI, ilutI, src, orbs, p_orb)
        integer, intent(in) :: nI(nel), src(2)
        integer(n_int), intent(in) :: ilutI(0:niftot)
        integer, intent(out) :: orbs(2)
        real(dp), intent(out) :: p_orb
#ifdef DEBUG_
        character(*), parameter :: this_routine = "pick_ab_orbitals_par_hubbard"
        real(dp) :: test, test_arr(nBasis / 2)
        integer :: ex(2, 2)
#endif
        real(dp) :: cum_arr(nbasis / 2)
        real(dp) :: cum_sum
        integer :: orb_list(nbasis / 2, 2)
        integer :: ind

        ! without transcorrelation factor this is uniform, but with a
        ! transcorrelation factor the matrix element might change and so also
        ! the pgen should change.

        call create_ab_list_par_hubbard(nI, ilutI, src, orb_list, cum_arr, cum_sum)

        if (cum_sum < EPS) then
            p_orb = 0.0_dp
            orbs(1) = ABORT_EXCITATION
            return
        end if

        ! this stuff is also written so often i should finally make a routine
        ! out of that
        call pick_from_cum_list(cum_arr, cum_sum, ind, p_orb)

        orbs = orb_list(ind, :)

#ifdef DEBUG_
        ! check that the other way of picking the orbital has the same
        ! probability..
        call create_ab_list_par_hubbard(nI, ilutI, src, orb_list, test_arr, cum_sum, &
                                        orbs(2), test)

        if (abs(test - p_orb) > 1.e-8) then
            print *, "pgen assumption wrong in ", this_routine
            print *, "nI: ", nI
            print *, "p_orb: ", p_orb
            print *, "test: ", test
            print *, "ij: ", src
            print *, "orbs: ", orbs
            print *, "cum_arr: ", cum_arr
            print *, "test_arr: ", test_arr
        end if

        !todo: also call the calc_pgen_k_space_hubbard here and check
        ! pgens
        ex(1, :) = src
        ex(2, :) = orbs

#endif

        ! do i have to recalc. the pgen the other way around? yes!
        ! effectively reuse the above functionality
        ! i am pretty sure i just have to find the position in the
        ! list.. OR: since in the hubbard it is just twice the
        ! probability or? i am pretty sure yes.. but for all of them..
        ! so in the end it shouldnt matter again..
        p_orb = 2.0_dp * p_orb

    end subroutine pick_ab_orbitals_par_hubbard