two_body_contrib_kvec Function

    HElement_t(dp) function two_body_contrib_kvec(J, p, k)
        ! two body contribution for the J optimization!
        real(dp), intent(in) :: J
        integer, intent(in) :: p(N_DIM), k(N_DIM)

        integer :: q(N_DIM)

        q = lat%subtract_k_vec(p, k)

        ! i still have to decide how to loop over the HF.. maybe i dont
        ! double count and then i dont need the /2 here!
        two_body_contrib_kvec = real(uhub, dp) / 2.0_dp - real(bhub, dp) * &
                                ((exp(J) - 1.0_dp) * epsilon_kvec(p) + (exp(-J) - 1.0) * epsilon_kvec(q))

    end function two_body_contrib_kvec