get_orb_from_kpoints_three Function

    function get_orb_from_kpoints_three(src, orb_a, orb_b) result(orb_c)
        ! momentum conservation for three-body terms
        integer, intent(in) :: src(3), orb_a, orb_b
        integer :: orb_c
#ifdef DEBUG_
        character(*), parameter :: this_routine = "get_orb_from_kpoints_three"
#endif
        integer :: sum_ms, kc(3), spin_c, spin_ab
        type(symmetry) :: kc_sym

        ! implement that generally for also all-spin parallel excitation, which
        ! might be necessary in the future..
        sum_ms = sum(get_spin_pn(src))

        ASSERT(sum_ms == -3 .or. sum_ms == -1 .or. sum_ms == 1 .or. sum_ms == 3)

        ! momentum conservation: ka + kb + kc = ki + kj + kk
        ! to this better with the symmetry table or a new addition
        ! functionality for k-vectors, to never leave the first BZ
        ! do something like:
        ! and write two routines, interfaced, one take a vector of
        ! k-values the other takes the encoded symmetry symbol
        ! and setup a symmetry table and inverse symmetry list like in the
        ! sym mod
        kc_sym = SymTable(G1(src(1))%sym%s, SymTable(G1(src(2))%sym%s, G1(src(3))%sym%s)%s)
        kc_sym = SymTable(kc_sym%s, SymConjTab(SymTable(G1(orb_a)%sym%s, G1(orb_b)%sym%s)%s))

        kc = lat%sym_to_k(kc_sym%s, :)

        !do it over the mult table now!!

        ! perdiodic BC:
        if (tHub) then
            if (t_k_space_hubbard) then
                ! use the new lattice implementation
                ! with the above correct addition it should still be in the
                ! first BZ..
            else
                call mompbcsym(kc, nBasisMax)
            end if
        end if

        ! and now we want the spin:
        if (sum_ms == -3) then
            ! assert we can still reach the desired spin:
            ASSERT(get_spin_pn(orb_a) + get_spin_pn(orb_b) == -2)

            spin_c = 1
        else if (sum_ms == 3) then
            ASSERT(get_spin_pn(orb_a) + get_spin_pn(orb_b) == 2)

            spin_c = 2

        else
            spin_ab = get_ispn([orb_a, orb_b])

            ! we know we are not totally parallel so if:
            if (spin_ab == 1) then
                ! if we have already two beta, we need alpha
                spin_c = 2
            else if (spin_ab == 3) then
                ! if we have two alpha we need a beta
                spin_c = 1

            else
                ! in this case we need a spin to reach the final ms
                if (sum_ms == -1) then
                    spin_c = 1
                else
                    spin_c = 2
                end if
            end if
        end if

        if (t_k_space_hubbard) then
            orb_c = lat%get_orb_from_k_vec(kc, spin_c)
        else
            orb_c = KPointToBasisFn(kc(1), kc(2), kc(3), spin_c)
        end if

    end function get_orb_from_kpoints_three