Given an ensemble of $N$ diatomic molecules, we know that the rotational partition function is given by $$Z_r = z_r^N$$ where $$z_r = \sum_{l} (2l+1)e^{-Kl(l+1)}$$ where $K = \beta \hbar ^2 /2I$.
I want to apply this to calculate the $C_v$ for an ensemble of $H_2$ (hydrogen) molecules in the low $T$ limit and the high $T$ limit for a certain quantum mechanical simulation. Before I run the simulation, I want to have an idea of how the results should look like.
For a general canonical ensemble, $C_v = \left( \partial E / \partial T \right)_V$, where $E = -(\partial Z / \partial \beta$). How would I evaluate the $C_v$ for a system of one mole of para-$H$ (singlet state) and one mole of ortho-$H$ (triplet state) in the high temperature limit and low temperature limit? How would I find the partition function for an HD molecule?
After doing some research, I see that $$Z_{para}= \sum_{l \,\text{even}} (2l+1)e^{-\beta \hbar ^2 l(l+1)/2I}$$ and $$Z_{ortho}= \sum_{l \,\text{even}} (2l+1)e^{-\beta \hbar ^2 l(l+1)/2I}$$
I know that because of the different multiplicities of these two states at equilibrium, the rotational partition function at low $T$ looks like
$$z_r = \frac{1}{4} \sum_{l \,\text{even}} (2l+1)e^{-\beta \hbar ^2 l(l+1)/2I} + \frac{3}{4} \sum_{l \, \text{odd}} (2l+1)e^{-\beta \hbar ^2 l(l+1)/2I} $$
After seeing this document, I can see how $C_v$ looks like for ortho, para-$H_2$, and their equilibrium. How do you find the partition function - and then $C_v$ for $HD$?
