two_body_contrib_ksym Function

    HElement_t(dp) function two_body_contrib_ksym(J, p, a)
        ! same as above just with symmetry symbols
        ! AND: we have  to do the p - k before calling this function!
        real(dp), intent(in) :: J
        type(symmetry), intent(in) :: p, a

        if (.not. t_symmetric) then
            two_body_contrib_ksym = real(uhub, dp) / 2.0_dp + real(bhub, dp) * &
                                    ((exp(J) - 1.0_dp) * epsilon_kvec(a) + (exp(-J) - 1.0) * epsilon_kvec(p))
        else
            two_body_contrib_ksym = real(uhub, dp) / 2.0_dp + real(bhub, dp) * &
                                    ((exp(J) - 1.0_dp) * epsilon_kvec(a) + (exp(-J) - 1.0) * epsilon_kvec(p))
        end if

    end function two_body_contrib_ksym