The trace can be calculated in any basis. Therefore, we will calculate it in the oscillator eigenbasis $\{|n\rangle\}$ as follows. First,
$$\text{Tr}\left( e^{-H/ k_B T} \right)
= \sum_{n=0}^{\infty}\langle n|e^{-H/ k_B T}|n\rangle
=\sum_{n=0}^{\infty}\langle n|e^{-\hbar\omega(n+1/2)/ k_B T}|n\rangle
=\sum_{n=0}^{\infty}e^{-\hbar\omega(n+1/2)/ k_B T}.
$$
This last sum can be evaluated because it can be written as a geometric series:
$$\text{Tr}\left( e^{-H/ k_B T} \right)
=e^{-\hbar\omega/2k_B T}\sum_{n=0}^{\infty}\left(e^{-\hbar\omega / k_B T}\right)^n
=e^{-\hbar\omega/2k_B T}\frac{1}{1-e^{-\hbar\omega / k_B T}}
=\frac{e^{\hbar\omega / 2k_B T}}{e^{\hbar\omega / k_B T}-1}
$$
Then, note that
$$\text{Tr}(\rho H)
= \text{Tr}\left(\frac{e^{-H / k_B T}}{\text{Tr}\left( e^{-H/ k_B T} \right) } H\right)
=\frac{\text{Tr}\left(e^{-H / k_B T}H\right)}{\text{Tr}\left( e^{-H/ k_B T} \right) }.
$$
So, we have to do the top:
\begin{align*}
\text{Tr}\left(e^{-H / k_B T}H\right)
&=\sum_{n=0}^{\infty}\langle n|e^{-H / k_B T}H|n\rangle
=\sum_{n=0}^{\infty}e^{-\hbar\omega(n+1/2)/ k_B T}\hbar\omega\left(n+\frac{1}{2}\right)\\
&=\frac{\hbar\omega}{2}\sum_{n=0}^{\infty}e^{-\hbar\omega(n+1/2)/ k_B T} + \hbar\omega\sum_{n=0}^{\infty}e^{-\hbar\omega(n+1/2)/ k_B T}n\\
&=\frac{\hbar\omega}{2}\frac{e^{\hbar\omega / 2k_B T}}{e^{\hbar\omega / k_B T}-1}
+ \hbar\omega e^{-\hbar\omega / 2k_B T}
\frac{e^{\hbar\omega / 2k_B T}}{e^{\hbar\omega / k_B T}-1}
\frac{1}{e^{\hbar\omega / k_B T}-1}
\end{align*}
Putting these together yields
$$\text{Tr} (\rho H) = \frac{\hbar \omega}{2} + \frac{\hbar \omega}{e^{\hbar \omega / k_B T} - 1} \, .$$