gen_triple_hubbard Subroutine

    subroutine gen_triple_hubbard(nI, ilutI, nJ, ilutJ, ex, tParity, pgen)
        ! i think i should calculat the matrix element in here already!
        ! in this case.. otherwise i have to carry the tParity and ex
        ! all the way through the rest of the code and this makes problems
        ! i guess, since triples are actually never considered .. todo
        integer, intent(in) :: nI(nel)
        integer(n_int), intent(in) :: ilutI(0:niftot)
        integer, intent(out) :: nJ(nel), ex(2, 3)
        integer(n_int), intent(out) :: ilutJ(0:niftot)
        logical, intent(out) :: tParity
        real(dp), intent(out) :: pgen
        integer :: elecs(3), orbs(3), src(3), sum_ms
        real(dp) :: p_elec, p_orb(2)

        ! first we pick 3 electrons in this case ofc.
        ! with the restriction, that they must not be all the same spin!
        call pick_three_opp_elecs(nI, elecs, p_elec, sum_ms)

        src = nI(elecs)
        ! then i pick 1 orbital? maybe..
        ! if this orbital is of the minority spin type, the other 2 orbitals
        ! MUST be of parallel spin..
        ! if the orbital is of the majority spin type the remaining
        ! orbitals must be of opposite spin type..
        ! can i take that into account with a tailored get_orb_from_kpoints()
        ! function for 3 electrons or should i decide here if we always
        ! want to do a specific picking order (restricting pgens, but making
        ! it easier to handle algorithmically) or if we want full flexibility
        ! (increasing pgens, but kind of making it a bit more difficult..)
        ! NO: we decide to always pick the minority spin first!
        call pick_a_orbital_hubbard(ilutI, orbs(1), p_orb(1), sum_ms)

        ! and pick the remaining two orbitals (essentially it is only
        ! one, since the last is restricted due to momentum conservation!)
        call pick_bc_orbitals_hubbard(nI, ilutI, src, orbs(1), orbs(2:3), p_orb(2))

        if (orbs(2) == 0) then
            nJ(1) = ABORT_EXCITATION
            pgen = 0.0_dp
            return
        end if

        ! so now.. did Robert make a routine like:
        call make_triple(nI, nJ, elecs, orbs, ex, tParity)

        ilutJ = make_ilutJ(ilutI, ex, 3)

        ! i am not sure about the factor of 2 here.. maybe this is already
        ! taken into account in the orbital creation
        pgen = p_elec * product(p_orb)

    end subroutine gen_triple_hubbard