Calculating quantum partition functions 
...By quantizing we the get the following Hamiltonian operator
$$\hat{H}=\sum_{\mathbf{k}}\hbar \omega(\mathbf{k})\left(\hat{n}(\mathbf{k})+\frac{1}{2} \right)$$
  where $\hat{n}(\mathbf{k})=\hat{a}^{\dagger}(\mathbf{k})\hat{a}(\mathbf{k})$ is the number operator of oscillator mode $\mathbf{k}$ with eigenvalues $n_{\mathbf{k}}=0,1,2,\dots$.
Using the quantum canonical ensemble show that the internal energy $E(T)$ is given by>
$$E(T)=\langle H \rangle = E_0 + \sum_{\mathbf{k}}\frac{\hbar \omega(\mathbf{k})}{e^{\beta\hbar \omega(\mathbf{k})}-1}$$
where $E_0$ is the sum of ground state energies of all the oscillators.

I started this by calculating the partition function 
$$\begin{align} Z &= \sum_{\Gamma}e^{-\beta \mathcal{H}(\Gamma)} \\
&= \sum_{\Gamma}e^{-\beta (\sum_{\mathbf{k}}\hbar \omega(\hat{n}(\mathbf{k})+\frac{1}{2}))}
\end{align}$$
($\Gamma$ is a microstate of the system)
but I cannot see the thought process behind evaluating these, particularly with respect to the summations. This is a common problem I have found.
I would then go on to use $E=-\frac{\partial \ln Z}{\partial \beta}$
 A: Quantum mechanically the general expression you want for the partition function is 
$$ Z = \mathrm{Tr} \left( \mathrm{e}^{-\beta H} \right),$$
where $\mathrm{Tr}$ means the trace (i.e. sum over micro-states). Now you can use the fact that the modes are independent, so that quantum Boltzmann operator $\mathrm{e}^{-\beta H}$ factorises into a product. This means that you can evaluate the trace over each oscillator mode separately:
$$ Z =  \mathrm{Tr} \left( \mathrm{e}^{-\beta H} \right) = \mathrm{Tr} \left( \prod_\mathbf{k}\mathrm{e}^{-\beta H_\mathbf{k}}\right) =  \prod_\mathbf{k}  \mathrm{Tr}_\mathbf{k} \left( \mathrm{e}^{-\beta H_\mathbf{k}}\right) = \prod_\mathbf{k} Z_\mathbf{k} $$
where $\mathrm{Tr}_\mathbf{k}$ means the trace over only the Hilbert space of mode $\mathbf{k}$, and
$$H_\mathbf{k} = \hbar\omega(\mathbf{k})\left(\hat{n}(\mathbf{k}) + \frac{1}{2}\right).$$
Now $\mathrm{Tr}_\mathbf{k}$ means simply averaging over all the possible states in the Hilbert space, which you might as well choose to be the eigenstates of the number operator $\hat{n}(\mathbf{k})\lvert m_\mathbf{k}\rangle= m_\mathbf{k} \lvert m_\mathbf{k}\rangle$, with $m_\mathbf{k} = 0,1,2,\ldots$. So you have to evaluate
$$ Z_\mathbf{k} = \mathrm{Tr}_\mathbf{k} \left( \mathrm{e}^{-\beta H_\mathbf{k}}\right)  = \sum_{m_\mathbf{k}=0}^\infty \langle m_\mathbf{k} \rvert \mathrm{e}^{-\beta H_\mathbf{k}} \lvert m_\mathbf{k} \rangle. $$
A: My original answer was incorrect and Mark gave a good answer above.
Here is a slightly different approach by considering the total energy of the system. Say we have $\alpha$ total harmonic oscillators. We can write the total energy of the whole system as
$$ E = \sum_\alpha \frac{1}{2}\hbar\omega_\alpha + \sum_\alpha n_\alpha \hbar \omega_\alpha .$$
For each oscillator $\alpha$, $n_\alpha = 0,1,2,3...$
The first term in the expression for the total energy is the sum of ground state energies of all the oscillators, $E_0$.
$$ \Rightarrow E_0 = \sum_\alpha \frac{1}{2}\hbar\omega_\alpha $$
The we can write the partition function as the sum of all possible subsytems that have this total energy.
$$ Z = \sum_{\{n_\alpha\}} e^{-\sum_\alpha \frac{1}{2}\hbar\omega_\alpha /kT}
e^{-\sum_\alpha n_\alpha \hbar \omega_\alpha /kT}$$
This simplifies to...
$$ Z = e^{-E_0/kT} \prod_\alpha \left( 1+e^{-\hbar\omega_\alpha /kT}+ e^{-2\hbar\omega_\alpha /kT}+e^{-3\hbar\omega_\alpha /kT}+... \right)
$$
Rewriting the series of exponential terms gives us...
$$ Z = e^{-E_0/kT} \prod_\alpha \left( \frac{1}{1-e^{-\hbar\omega_\alpha /kT}}\right)$$
From here you can use the formula $U = -\frac{\partial \ln Z}{\partial\beta}$ to find the internal energy.
