QFI must always be computed with respect to a parameter "$\theta$". Perhaps it is the temperature that you want to use here, or $\beta$?
Regardless, we can put the state into your desired form by expanding the exponential. Taking the complete set of orthonormal basis states $|\mathbf{n}\rangle=|n_1\rangle\otimes |n_2\rangle\otimes\cdots$ that satisfy $a_i^\dagger a_i|\mathbf{n}\rangle=n_i|\mathbf{n}\rangle$, the state can be rewritten as
$$\rho=\sum_{\mathbf{n}}\frac{\exp[-\beta\hbar\omega(\sum_i n_i+\frac{1}{2})]}{Z}|\mathbf{n}\rangle\langle\mathbf{n}|.$$ All of the probabilities are thus given by $p_{\mathbf{n}}=\frac{\exp[-\beta\hbar\omega(\sum_i n_i+\frac{1}{2})]}{Z}$ and the eigenstates do not depend on any parameter, so you can use classical Fisher information formulas to proceed:
$$I_Q(\rho,\theta)=\sum_{\mathbf{n}}p_{\mathbf{n}}\left(\frac{\partial \ln p_{\mathbf{n}}}{\partial \theta}\right)^2.$$
Now, if you would like a general formula for the QFI $I_Q$ for mixed states that doesn't require this canonical form, there are many options, e.g. collected in this review. If you have an eigen-decomposition of your state $\rho=\sum_k p_k|k\rangle\langle k|$, you can use
$$I_Q(\rho,\theta)=2\sum_{jk}\frac{|\langle j|\partial_\theta\rho|k\rangle|^2}{p_j+p_k}.$$
If you do not, you can use
$$I_Q(\rho,\theta)=\min_\Psi(\langle \partial_\theta\Psi|\partial_\theta\Psi\rangle-|\langle \Psi|\partial_\theta\Psi|^2\rangle)=\min_\Psi(\langle \partial_\theta\Psi|\partial_\theta\Psi\rangle\rangle)$$ where $|\Psi\rangle$ is a purification of $\rho$ in a larger Hilbert space, or a Lyapunov representation
$$I_Q(\rho,\theta)=2\int_0^\infty ds \mathrm{Tr}[(\partial_\theta \rho)e^{-\rho s}(\partial_\theta \rho)e^{-\rho s}].$$
If $\rho$ is invertible (i.e., full-rank) then we can solve the Lyapunov equation to find
$$I_Q(\rho,\theta)=2\mathrm{vec}(\partial_\theta \rho)^\dagger (\rho^*\otimes \mathbb{I}+\mathbb{I}\otimes \rho)^{-1}\mathrm{vec}(\partial_\theta \rho)$$ and if it is not full-rank then we can add a small component of a maximally mixed state and take the limit of that component going to zero (https://doi.org/10.48550/arXiv.1801.00945).