7/2 versus 9/2 for diatomic heat capacity Question
I calculated the classical heat capacity of a diatomic gas as $C_V = (9/2)Nk_B$, however the accepted value is $C_V = (7/2)Nk_B$.
I assumed the classical Hamiltonian of two identical atoms bound together as
$$
H = \dfrac{1}{2m}( |\bar{p}_2|^2 + |\bar{p}_2|^2)+
\dfrac{\alpha}{2} |\bar{q}_1-\bar{q}_2|^2.
$$
I calculated the partition function of $N$ particles as
$$
Z
=
\left(
\iiint_{-\infty}^{\infty}
\iiint_{-\infty}^{\infty}
\iiint_{-\infty}^{\infty}
\iiint_{-\infty}^{\infty}
e^{-\beta H}
~d^3q_1~d^3p_1~d^3q_2~d^3p_2
\right)^N
\propto
V^N T^{(9/2)N}.
$$
I calcuated the heat capacity as
$$
C_V
=
\dfrac{\partial }{\partial T}
\left(
k_B
T^2
\dfrac{\partial \ln(Z)}{\partial T}
\right)
=
\dfrac{9}{2}k_BN.
$$
Why does the classical argument fail?
Classical Derivation
The partition function is
\begin{align}
Z
&=&
\left(
\frac{1}{h^6} \int \mathrm{e}^{- \beta H(\bar{q}_1,\bar{q}_2,\bar{p}_1,\bar{p}_2)} ~d^{3}q_1 ~d^{3}q_2 ~d^{3}p_1   ~d^{3}p_2
\right)^N
\\&=&
\left(
\frac{1}{h^6} \int \mathrm{e}^{- \beta ((|\bar{p}_1|^2+|\bar{p}_2|^2)/(2m)+\alpha |\bar{q}_1-\bar{q}_2|^2/2)} ~d^{3}q_1 ~d^{3}q_2 ~d^{3}p_1   ~d^{3}p_2
\right)^N
\end{align}
A useful gaussian integral
\begin{align}
\int_{-\infty}^{\infty} e^{-\gamma (x-x_0)^2}dx = \sqrt{\dfrac{\pi}{\gamma}}
\end{align}
The partition function can be evaluated using separated integrals
\begin{align}
\iiint_{-\infty}^{\infty} \mathrm{e}^{- \beta |\bar{p}_1|^2} ~d^{3}p_1 = \iiint_{-\infty}^{\infty} \mathrm{e}^{- \beta |\bar{p}_2|^2} ~d^{3}p_2 = \left(\sqrt{\dfrac{\pi}{\beta}}\right)^3
\end{align}
and
\begin{align}
\iiint_{-\infty}^{\infty} \iiint_{-\infty}^{\infty} \mathrm{e}^{- \beta \alpha |\bar{q}_1-\bar{q}_2|^2/2 } ~d^{3}q_1 ~d^{3}q_2
=
\left(
\sqrt{\dfrac{\pi}{\beta \alpha/2}}
\right)^3
\iiint_{-\infty}^{\infty}  ~d^{3}q_1
=
\left(
\sqrt{\dfrac{\pi}{\beta \alpha/2}}
\right)^3
V
\end{align}
The last set of integrals are improper integrals. One has to take the limit as the space approaches infinite containment. In that limit, integrating one set of variables $d^3q_2$ approaches the limit of a finite Gaussian term, while the other $d^3q_1$ approaches the diverging value of the total volume of the gas.
The partition function is
\begin{align}
Z
&=&
\left(
h^{-6}
\left(\sqrt{\dfrac{\pi}{\beta}}\right)^3
\left(\sqrt{\dfrac{\pi}{\beta}}\right)^3
\left(
\sqrt{\dfrac{\pi}{\beta \alpha/2}}
\right)^3
V
\right)^N
\\&=&
\left(
h^{-6}
\left(k_B T \pi\right)^{9/2}
\left(
\dfrac{2}{\alpha}
\right)^{3/2}
V
\right)^N
\\&=&
\left(
h^{-6}
\left(k_B \pi\right)^{9/2}
\left(
\dfrac{2}{\alpha}
\right)^{3/2}
\right)^N
V^N
T^{9N/2}
\end{align}
 A: The potential energy for a diatomic molecule is not
$$
U(\vec{q}_1, \vec{q}_2) = \frac{\alpha}{2} |\vec{q}_1 - \vec{q}_2|^2
$$
but is instead
$$
U(\vec{q}_1, \vec{q}_2) = \frac{\alpha}{2} (|\vec{q}_1 - \vec{q}_2| - r_0)^2,
$$
where $r_0$ is the equilibrium bond distance.  The important difference here is that in your version, any displacement of the vector $\vec{q}_1 - \vec{q}_2$ will result in a quadratic change in the potential energy;  whereas in the correct version, there will be two directions in "configuration space" that correspond to no change in the potential energy.  Remember that the equipartition theorem basically says that every degree of freedom that contributes quadratically to the energy will then contribute $\frac{1}{2} k$ to $C_V$.  These two spurious energetic degrees of freedom are what are giving you $C_V = \frac{9}{2} k N$ instead of $C_V = \frac{7}{2} k N$.
Just to show that I'm not making this up, let's do the integral.  Define $\vec{Q} = \frac{1}{2}(\vec{q}_1 + \vec{q}_2)$ and $\vec{r} = \vec{q}_1 - \vec{q}_2$.
$$
\begin{align}
I = \iiint_{-\infty}^{\infty} \iiint_{-\infty}^{\infty} \mathrm{e}^{- \beta \alpha (|\bar{q}_1-\bar{q}_2|-r_0)^2/2 } ~d^{3}q_1 ~d^{3}q_2
&= \iiint_{-\infty}^{\infty} \iiint_{-\infty}^{\infty} \mathrm{e}^{- \beta \alpha (r-r_0)^2/2 } ~d^{3}Q ~d^{3}r \\
&= \left[ \iiint_{-\infty}^{\infty}~d^{3}Q \right] \left[ \iiint_{-\infty}^{\infty} \mathrm{e}^{- \beta \alpha (r-r_0)^2/2 } ~d^{3}r \right]
\end{align}
$$
The first integral gives a factor of $V$ as before.  The second one is a little more complicated.  The angular contribution is obviously $4\pi$, leaving 
$$
I = 4 \pi V \int_0^\infty r^2 \mathrm{e}^{- \beta \alpha (r-r_0)^2/2 } ~dr
$$
This last integral isn't of the standard "useful Gaussian integral" form, and will not give a result that is exactly proportional to $\beta^{-1/2}$.  However, in the limit of low temperature, it does approach this limit.  Define $\tilde{r} = \sqrt{\beta \alpha} (r - r_0)$;  then the integral becomes
$$
I = \frac{4 \pi V}{\sqrt{\beta \alpha}} \int_{-\sqrt{\beta \alpha} r_0}^\infty \left( \frac{\tilde{r}}{\sqrt{\beta \alpha}} + r_0 \right)^2 e^{-\tilde{r}^2/2} \, d \tilde{r}.
$$
In the low-temperature limit, we have $\beta \to \infty$, meaning that the lower limit of integration becomes $- \infty$ and the first term in the parentheses vanishes;  thus, in this limit,
$$
I \approx \frac{4 \sqrt{2} \pi^{3/2} V r_0^2 }{\sqrt{\beta \alpha}} \propto V T^{1/2} 
$$
as desired.  
EDIT: The exact integral above can't actually be evaluated in closed form, but it can be expressed in terms of the normalized error function erf(x):
$$
I = \frac{4 \pi^{3/2} V}{\sqrt{2}} \left[ \left(\frac{r_0^2}{\sqrt{\alpha \beta}} + \frac{1}{(\alpha \beta)^{3/2}} \right)\left( 1 + \text{erf} \left( \frac{r_0 \sqrt{\alpha \beta}}{\sqrt{2}} \right) \right) + \sqrt{\frac{2}{\pi}} \frac{r_0}{\alpha \beta} e^{-\alpha \beta r_0^2/2}\right].
$$
Note that if we set $r_0 \to 0$, we recover your result above (with $I \propto T^{3/2}$.)  However, for non-zero $r_0$, we get a leading-order result proportional to $\sqrt{T}$, and a leading-order correction proportional to $T^{3/2}$ (as well as even smaller corrections proportional to $e^{-\alpha \beta r_0^2/2}$ times various powers of $T$, arising from the exponential term and the asymptotic expansion of the erf function.)
