First integral of the Kepler problem Consider the motions of a bounded particle which is under the influence of the 
gravitional interaction of a second particle fixed at the origin
$$
\ddot q = -\nabla V(q) 
$$
where $V(q) = - \frac{\mu}{|q|}$.
Usually we define the angular momentum and the eccentricity vector as
$$
C := q \times \dot q \qquad  e:=-\frac{q}{|q|} + 1/\mu\; p \times C
$$
Is it common to define $C$ and $e$ where the potential $V$ may contains for instance higher terms in the potential (like the J2 term for the Earth) ?
Would it be the usual way to define orbital elements for this Kepler problem with
higher order terms ?
 A: Yes, it does make sense to define orbital elements even for potentials with higher order terms of $q$, which I'll rename into $r$ since it's more common in these types of problems.
Example: let the potential be given by $U(r)=-\frac{\mu}{r}+\frac{h}{r^2}$, where $h>0$. For the initial conditions of the particle having mass $m$, with an initial position of $T(a,0)$ and initial velocity of $v_0=\sqrt{\frac{\mu}{2am}}$ in the direction $\phi_0=\pi/4$, and let the relation between $h$ and other parameters be given as $h=\frac{3}{8}ak$ (parameters are chosen to give a nice output formula and graph), the orbit equation can be calculated to be 
$$r(\theta)=\frac{a}{1-0.5\sin(2\theta)}$$
and looks like this: 

On the image you may see the initial position and velocity $\vec{r}_0$ and $\vec{v}_0$ respectively, as well as $r_{min}$ and $r_{max}$ which correspond to perihelions and aphelions.
As for defining angular momentum and the Laplace-Runge-Lenz vector (which you call $\vec{e}$), since the potential has rotational symmetry (in other words, the force is central), angular momentum is conserved (proof of which I may include if you want), so you may calculate it from initial conditions: 
$$\vec{M}=\vec{r}_0\times\vec{v}_0=mav_0\sin(\pi/4)\hat{z}=0.5\sqrt{mak}\hat{z}$$ 
The L-R-L vector is a conserved quantity in the Kepler problem, but in other ones it's usually not conserved, and has a time evolution which can be obtained perturbatively. You may find more information here.
With more higher-order terms the integrals of motion usually become very hard or unsolvable, but even in the case where we resort to numerical calculations, a lot of information can be obtained from those two vectors.
