When you say you "got the first to terms" I guess that means the non-interacting part. So, your trouble seems to be with the interaction terms. I don't want to work the whole thing out since this looks like homework, but maybe this will help:
Probably best to start by working out the commutator:
$$
[\sum_{k_1,k_2,q} V_q a_{k_1+q}^{\dagger} a_{k_2+q}^\dagger a_{k_2} a_{k_1},a_{k-Q}^\dagger a_{k}]
$$,
which it not really so horrible.
You can just repeatedly use the fact that:
$$
[AB,C] = A[B,C] + [A,C]B
$$
I think it might help to rewrite the commutator of interest like:
$$
[a^\dagger_1 a^\dagger_2 a_3 a_4,a_5^\dagger a_6]
$$
since it will make it harder to make transcription errors in the writing.
You can do this reduction many ways. For example, you could start off like:
$$
[a^\dagger_1 a^\dagger_2 a_3 a_4,a_5^\dagger a_6]
=
a^\dagger_1 [a^\dagger_2 a_3 a_4,a_5^\dagger a_6]
+
[a^\dagger_1 ,a_5^\dagger a_6]a^\dagger_2 a_3 a_4
$$
$$
=
a^\dagger_1 [a^\dagger_2 a_3 a_4,a_5^\dagger a_6]
+
a_5^\dagger[a^\dagger_1 , a_6]a^\dagger_2 a_3 a_4
$$
$$
=
a^\dagger_1 (a^\dagger_2[ a_3 a_4,a_5^\dagger a_6]+[a^\dagger_2 ,a_5^\dagger a_6]a_3 a_4)
+
a_5^\dagger[a^\dagger_1 , a_6]a^\dagger_2 a_3 a_4
$$
$$
=
a^\dagger_1 (a^\dagger_2[ a_3 a_4,a_5^\dagger ]a_6+a_5^\dagger[a^\dagger_2 , a_6]a_3 a_4)
+
a_5^\dagger[a^\dagger_1 , a_6]a^\dagger_2 a_3 a_4
$$
$$
=
a^\dagger_1 (a^\dagger_2(a_3[ a_4,a_5^\dagger ]+[ a_3 ,a_5^\dagger ]a_4)a_6+a_5^\dagger[a^\dagger_2 , a_6]a_3 a_4)
+
a_5^\dagger[a^\dagger_1 , a_6]a^\dagger_2 a_3 a_4
$$
$$
=
a^\dagger_1 (a^\dagger_2(a_3\delta_{4,5}+\delta_{35}a_4)a_6-a_5^\dagger\delta_{26}a_3 a_4)
-
a_5^\dagger\delta_{16}a^\dagger_2 a_3 a_4
$$
$$
=
a^\dagger_1 a^\dagger_2 a_3 a_6\delta_{4,5}+a^\dagger_1 a^\dagger_2 a_4 a_6\delta_{35}-a^\dagger_1 a_5^\dagger a_3 a_4\delta_{26}
-
a_5^\dagger a^\dagger_2 a_3 a_4\delta_{16}
$$
Ok, so now what? I guess now you can start using the Hartree Fock approximation to figure out the expectation value for:
$$
< a^\dagger_1 a^\dagger_2 a_3 a_6\delta_{4,5}+a^\dagger_1 a^\dagger_2 a_4 a_6\delta_{35}-a^\dagger_1 a_5^\dagger a_3 a_4\delta_{26}
-
a_5^\dagger a^\dagger_2 a_3 a_4\delta_{16}
>$$
Then substitute back in the correct index names and do all the summations.
In some cases we can approximate the sum by an integral appropriate:
$$\sum_j e^{-\beta\epsilon_j} =\int e^{-\beta E}D(\epsilon)d\epsilon=\int e^{-\beta E(\omega)}D(\omega)d\omega$$
In fact that is how we determine the density function. And I think there is no reason we can write $\sum_i Z_{\omega_i}=\int Z_\omega D(\omega)d\omega$ with the same $D(\omega)$ from the egality above.
And also there is no reason to believe that $Z=\sum_i Z_{\omega_i}$, in fact $Z=\prod_i Z_{\omega_i}$ like you have demonstrated. So that $\ln Z=\sum_i \ln Z_{\omega_i}$, therefore:
$$\ln Z=\int (\ln Z_{\omega_i}) D(\omega)d\omega $$
In fact $D(\omega)d\omega$ counts the number of states having frequency between $\omega$ and $\omega+d\omega$, and similar for $D(\epsilon)d\epsilon$. You can imagine that
$$Z=\prod_{\omega_i}Z_{\omega_i}^{D(\omega)d\omega}$$
which mean we have approximate all the frequency in $[\omega,\omega+d\omega]$ by $\omega$.
Best Answer
Let us first consider the trace of the density matrix: $$ \text{Tr}\; \exp\left\{-\beta H\right\}=\sum_{n_1,n_2,\ldots}\langle n_1,n_2,\ldots|\exp\left\{-\beta\sum_k\left(a_k^\dagger a_k+\frac{1}{2}\right)\hbar\omega_n\right\}| n_1,n_2,\ldots\rangle\\ =\prod_k \sum_{n_1,n_2,\ldots}\langle n_1,n_2,\ldots|\exp\left\{-\beta\left(a_k^\dagger a_k+\frac{1}{2}\right)\hbar\omega_n\right\}| n_1,n_2,\ldots\rangle\\ =\prod_k\left(\sum_{n_k}\langle n_k|\exp\left\{-\beta\left(a_k^\dagger a_k+\frac{1}{2}\right)\hbar\omega_k\right\}|n_k\rangle\right) $$ Crucially, the different harmonic oscillators decouple and if you compute the expectation value of some operator $\mathcal{O}$ acting only on one of the oscillators all the other oscillator contributions will cancel in $\langle \mathcal{O} \rho\rangle/\langle \rho\rangle$. Let us ignore the vacuum energy for now and compute $$ \sum_{n_k}\langle n_k|\exp\left\{-\beta a_k^\dagger a_k\hbar\omega_k\right\}|n_k\rangle\\ =\sum_{n_k} \exp\left\{-\beta n_k \hbar\omega_k\right\}\\ =\sum_{n_k} \left(\exp\left\{-\beta \hbar\omega_k\right\}\right)^{n_k}\\ =\frac{1}{1-e^{-\beta \hbar\omega_k}}. $$ Adding the vacuum energy you get $$ \sum_{n_k}\langle n_k|\exp\left\{-\beta \left(a_k^\dagger a_k+\frac{1}{2}\right) \hbar\omega_k\right\}|n_k\rangle\\ =\frac{e^{-\frac{1}{2}\hbar\omega_k}}{1-e^{-\beta \hbar\omega_k}}. $$ Now, we tackle the actual expectation value you want to find: $$ \langle a_1^\dagger a_1\rangle=\frac{\text{Tr}\;a_1^\dagger a_1 \rho}{\text{Tr}\;\rho}. $$
As stated above, only the contribution from the first oscillator will not cancel. We need $$ \sum_{n_1}\langle n_1|a_1^\dagger a_1\exp\left\{-\beta\left(a_1^\dagger a_1+\frac{1}{2}\right)\hbar\omega_1\right\}|n_1\rangle\\ =\sum_{n_1} n_1 e^{-\beta\left(n_1+\frac{1}{2}\right)\hbar\omega_1}\\ =e^{-\frac{1}{2}\beta\hbar\omega_1} \left(-\frac{1}{\hbar\omega_1}\right)\frac{\partial}{\partial\beta}\sum_{n_1} e^{-\beta n_1 \hbar\omega_1}\\ =e^{-\frac{1}{2}\beta\hbar\omega_1}\left(-\frac{1}{\hbar\omega_1}\right)\frac{\partial}{\partial\beta}\frac{1}{1-e^{-\beta \hbar\omega_1}}\\ =e^{-\frac{1}{2}\beta\hbar\omega_1}\frac{e^{-\beta \hbar\omega_1}}{\left(1-e^{-\beta \hbar\omega_1}\right)^2} $$
Now we just divide by the contribution to the trace of the density matrix from the first oscillator, which is $$ \frac{e^{-\frac{1}{2}\hbar\omega_k}}{1-e^{-\beta \hbar\omega_k}}, $$ and we get $$ \langle a_1^\dagger a_1\rangle = \frac{e^{-\beta \hbar\omega_1}}{1-e^{-\beta \hbar\omega_1}}=\frac{1}{e^{\beta \hbar\omega_1}-1}. $$ This differs from your result by a minus sign but I think mine should be correct since your result would give negative values for the expectation value of the number operator.
The crucial point was: Both the numerator and denominator in $$ \langle\mathcal{O}\rangle=\frac{\text{Tr}\;\mathcal{O}\rho}{\text{Tr}\;\rho} $$ look somewhat like $$ \prod_k\left(\sum_{n_k}\langle n_k|\exp\left\{-\beta\left(a_k^\dagger a_k+\frac{1}{2}\right)\hbar\omega_k\right\}|n_k\rangle\right)=\prod_k \frac{e^{-\frac{1}{2}\hbar\omega_k}}{1-e^{-\beta \hbar\omega_k}}, $$ except the numerator gets changed in the factors where $\mathcal{O}$ is acting. All the other factors cancel.