pick_bc_orbitals_hubbard Subroutine

    subroutine pick_bc_orbitals_hubbard(nI, ilutI, src, orb_a, orbs, p_orb)
        ! this is the main routine, which picks the final (b,c) orbital for
        ! the 3-body excitation
        integer, intent(in) :: nI(nel), src(3), orb_a
        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_bc_orbitals_hubbard"
        real(dp) :: test
#endif
        real(dp) :: cum_arr(nbasis / 2), cum_sum
        integer :: orb_list(nbasis / 2, 2), ind

        ! do it similar to the other routines..
        call create_bc_list_hubbard(nI, ilutI, src, orb_a, orb_list, cum_arr, cum_sum)

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

        call pick_from_cum_list(cum_arr, cum_sum, ind, p_orb)

        orbs = orb_list(ind, :)

#ifdef DEBUG_
        ! the influence of orb_a is important in the pgen recalc!!
        call create_bc_list_hubbard(nI, ilutI, src, orb_a, orb_list, cum_arr, cum_sum, &
                                    orbs(2), test)

        if (abs(test - p_orb) > 1.e-8) then
            print *, "pgen assumption wrong: in ", this_routine
            print *, "p_orb: ", p_orb
            print *, "test: ", test
            print *, "ijk: ", src
            print *, "a: ", orb_a
            print *, "orbs: ", orbs
            print *, "cum_arr: ", cum_arr
            print *, "orb_list(:,1): ", orb_list(:, 1)
            print *, "orb_list(:,2): ", orb_list(:, 2)
        end if
#endif

        p_orb = 2.0_dp * p_orb

    end subroutine pick_bc_orbitals_hubbard