Normally we do NOT calculate the phase space density of a system. In the phase space formulation of classical statistical mechanics, the phase space density $\rho(p,q;t)$ has its specified form for different ensembles. Normally for systems at equilibrium the density $\rho$ has no explicit time dependence and thus we work with $\rho(p,q)$.
(1) For microcanonical ensemble, the energy of the system is fixed and according to the 'equal a priori probability', each microstate is equally likely, therefore, $\rho=constant$ throughout the relevant region in the phase space. The exactly value for this constant is not important since in statistical mechanics all we care about the relative chance or probabilities, you can always normalize it in the end. More precisely, for the average value of some physical quantity $f(p,q)$, we have
$$<f>=\frac{\int dqdp f(p,q)\rho}{\int dqdp \rho}=\frac{\int dqdp f(p,q)}{\int dqdp},$$
since $\rho$ is constant in this case.
(2) For canonical ensemble, which you might be more interested in, the density $\rho(p,q)$ is no longer uniform, instead we have
$$\rho(p,q)\propto exp(-H(p,q)/kT),$$
which simply means the higher the energy, the less likely the level will be occupied. This is just the usual Boltzmann factor to be consistent with the standard thermodynamic results. Then for different systems you just need to put in their respective Hamiltonians. Again, the proportionality constant in front is not important, it will be cancelled out anyway when we do the average calculation.
To summarize, the phase space density simply characterizes the relative probability density for the classical states labelled by $(p,q)$ with Hamiltonian $H(p,q)$ to be occupied. But we already have the probability density from standard thermodynamics, which can be directly used in the phase space approach.