# gen_triple_hubbard Subroutine

## public subroutine gen_triple_hubbard(nI, ilutI, nJ, ilutJ, ex, tParity, pgen)

### Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nI(nel)
integer(kind=n_int), intent(in) :: ilutI(0:niftot)
integer, intent(out) :: nJ(nel)
integer(kind=n_int), intent(out) :: ilutJ(0:niftot)
integer, intent(out) :: ex(2,3)
logical, intent(out) :: tParity
real(kind=dp), intent(out) :: pgen

## Source Code

    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