Deriving the partition function for a harmonic oscillator
I'm confused why you're interpreting the partition function as a count of states. It can't be a count; it's continuous.

The zero point energy doesn't actually matter because you can just shift the energy scale so that it starts at zero. The main point of zero point energy is that the ground state of the harmonic oscillator is such that there is energy, and the system is not stationary. I'm going to use it below anyway because you are.

You can get the answer you want, but you'll want to look at the probability $P_i$ of being in state $i$, where the partition function is used to normalize:

$$P_{i} = \frac{\epsilon_{i}}{q} = \frac{\mathrm{e}^{-\beta\hbar\omega (i + 1)/2}}{\sum_{i}\mathrm{e}^{-\beta\hbar\omega (i+1)/2}}$$

Substitution with your convergent sum:

$$P_{i} = \mathrm{e}^{-\beta\hbar\omega (i+1)/2} \frac{1 - \mathrm{e}^{-\beta\hbar\omega}}{\mathrm{e}^{-\beta\hbar\omega/2}} = \mathrm{e}^{-\beta\hbar\omega i}(1- e^{-\beta\hbar\omega})$$

For $T\rightarrow 0$, $P_{i} = \delta_{i0}$, which is exactly what you are looking for.