function get_3_body_helement_ks_hub(ex, tpar) result(hel)
! the 3-body matrix element.. here i have to be careful about
! the sign and stuff.. and also if momentum conservation is
! fullfilled ..
integer, intent(in) :: ex(2, 3)
logical, intent(in) :: tpar
HElement_t(dp) :: hel
integer :: ms_elec, ms_orbs, opp_elec, opp_orb, par_elecs(2), par_orbs(2)
logical :: sgn
type(symmetry) :: p_sym, hole_sym, k_sym, k1_sym, k2_sym
type(symmetry) :: ka_sym, kb_sym, kc_sym, kd_sym
hel = h_cast(0.0_dp)
ms_elec = sum(get_spin_pn(ex(1, :)))
ms_orbs = sum(get_spin_pn(ex(2, :)))
! check spin:
if (.not. (ms_elec == ms_orbs)) return
if (.not. (ms_elec == 1 .or. ms_elec == -1)) return
! check momentum conservation:
if (.not. check_momentum_sym(ex(1, :), ex(2, :))) return
! i have to get the correct momenta for the epsilon contribution
! see the sheets for the k-vec relations.
! we need the momentum p + (s-b) or variations thereof..
! p, we know, since it is the momentum of the minority spin-electron
! and (s-b) is the difference of an majority spin electron and the
! fitting hole, so the total momentum conservation a + b + c = s + p + q
! is fulfilled
! the same ofc is a + (c-q)
! is the minority hole always in ex(2,1)? otherwise we have to find it
! it is not!
! i think i have figured it out with the help of Manu
! the k-vector of the minority spin is always involved
! but of the electron.. or can we transform it?
! any way we have to calculate
! W(k_p + k_s - k_b) - W(k_p + k_q - k_b)
! for this we have to figure out what the minority and majority
! electrons are!
opp_elec = find_minority_spin(ex(1, :))
! although i really can't be sure about the minority whole always
! being at the first position in ex(2,:)..
opp_orb = find_minority_spin(ex(2, :))
p_sym = G1(opp_elec)%sym
hole_sym = G1(opp_orb)%sym
par_elecs = pack(ex(1, :), ex(1, :) /= opp_elec)
par_orbs = pack(ex(2, :), ex(2, :) /= opp_orb)
k_sym = SymTable(p_sym%s, hole_sym%s)
! we have to define an order here too
par_elecs = [minval(par_elecs), maxval(par_elecs)]
par_orbs = [minval(par_orbs), maxval(par_orbs)]
! BZ conserving addition:
k1_sym = SymTable(G1(par_orbs(1))%sym%s, SymConjTab(G1(par_elecs(1))%sym%s))
k2_sym = SymTable(G1(par_orbs(1))%sym%s, SymConjTab(G1(par_elecs(2))%sym%s))
! need to do the correct k additions!
! for some reason the compiler does not recognize the output of
! add_k_vec and subtract_k_vec as a vector...
! so do it intermediately
ka_sym = SymTable(hole_sym%s, k1_sym%s)
kb_sym = SymTable(hole_sym%s, k2_sym%s)
kc_sym = SymTable(k1_sym%s, SymConjTab(p_sym%s))
kd_sym = SymTable(k2_sym%s, SymConjTab(p_sym%s))
hel = three_body_prefac * ( &
epsilon_kvec(ka_sym) - epsilon_kvec(kb_sym) + &
epsilon_kvec(kc_sym) - epsilon_kvec(kd_sym))
! i have to decide on a sign here depending on the order of the
! operators.. todo!
sgn = get_3body_sign(ex)
if (.not. sgn) hel = -hel
if (tpar) hel = -hel
end function get_3_body_helement_ks_hub