The Hamiltonian of a real gas is usually taken of the form
$$
H(p_1,\dots,p_N,q_1,\dots,q_N) = \sum_{i=1}^N \frac{p_i^2}{2m} + \sum_{i<j} V(q_i,q_j).
$$
Since it is the sum of a momenta-dependent term and a position-dependent term, the canonical probability density (at temperature $T$) factorizes:
$$
f^{\rm can}_T(p_1,\dots,p_N,q_1,\dots,q_N) = f^{{\rm can},p}_T(p_1,\dots,p_N)f^{{\rm can},q}_T(q_1,\dots,q_N),
$$
where $f^p_T(p_1,\dots,p_N)$ is the marginal density for the momenta and is simply given by a Gaussian density, independently of the interaction $V$:
$$
f_T^{{\rm can},p}(p_1,\dots,p_N) = (2\pi m/\beta)^{-N/2} e^{-\beta\sum_{i}p_i^2/2m}.
$$
In particular, the distribution of the momenta is the same as for an ideal gas, which implies that the Maxwellian distribution applies regardless of the interaction term.
The distribution of the positions is however quite complicated. This makes an approach along the lines you want very intricate. The usual way to compute the pressure is thus via the virial expansion, which leads to the virial equation
$$
PV=Nk_BT\left( 1 + \frac{N}{V}B_2(T) + \frac{N^2}{V^2}B_3(T) + \frac{N^3}{V^3}B_4(T)+ \cdots \right),
$$
where the virial coefficients $B_i(T)$ are given by explicit expressions and can be (in principle) computed.
Some references:
- You also asked whether the van der Waals equation can be deduced from theoretical considerations or whether it is of a purely phenomenological nature. This equation can indeed be deduced, as an approximation by comparing it to the (exact) virial expansion described above. Alternatively, it is possible to derive it as a type of mean-field limit. A possible reference is our book, where the latter is discussed in Chapter 4 and the former in Chapter 5), although we only discuss the lattice gas.
- A very interesting (although rather old) book where the kind of derivation you are after can be found is this one.
