I'm trying to figure out how to write the partition function as a function of temperature for just a single hydrogen atom with a bound electron.
I know the Energies will be $$E_n=\frac{-13.6eV}{n^2}$$ and the partition function is $$Z = \sum_{n=1}^\infty e^{\beta(13.6eV)/n^2}$$ but I think I must have gone wrong somewhere because that sum diverges but I don't what else to do...
There are too many states (an infinite number of states) accumulating near the zero energy level, and thus this level (given temperature $T>0$) acts, as time goes on, as strong "attractor" of probability when the initial state is one of negative energy. This implies in equilibrium, if it exists, the states at or above zero should be occupied, and those below the zero energy level should be unoccupied, even though they have lower energy. The zero or positive energy states "win" the probability from the negative energy states, because there is so many more of them in the partition sum.
Also, there are infinitely many states at any interval of positive energies (describing ionized states).
A lone atom interacting with background radiation of temperature $T>0$ is not stable (in the long term, in the sense of time averages), because equilibrium radiation at any temperature will eventually kick the system into one of the ionized states, and there are too many of those, so the atom will stay there with probability one (assuming low temperature, this decay may take a long time though). Thus if single hydrogen atom is in contact with a thermal reservoir at non-zero temperature, after long enough time, it is almost certainly ionized.
The partition function of a single hydrogen atom in infinite empty space thus indeed diverges, you've got the right result, only in the sum, you are missing terms due to positive energies. The "correct" sum should be $$ \sum_{n=1}^\infty g_n e^{\beta \frac{hcR}{n^2}} + \int_0^\infty G(E)e^{-\beta E}dE $$ where $g_n$ is degeneracy of the level $E_n = -\frac{hcR}{n^2}$ and $G(E)dE$ gives some measure of the number of states with energy in the infinitesimal interval $E,E+dE$.
If the hydrogen atom particles (nucleus, electron) are restricted to a finite volume of space (by means of some walls or by presence of other hydrogen atoms nearby, such as when the atom is part of a gas of hydrogen atoms with finite non-zero volume density), the structure of the Hamiltonian eigenstates changes; the energy levels no longer accumulate overwhelmingly around zero, so density of states around zero and above is no longer quasi-continuous, but is discrete instead. The easy way to see this is plausible is to recall that the spectrum of eigenvalues for a particle in a box is discrete, and depends on the quantum numbers as $n_1^2+n_2^2+n_3^2$. Then the partition function may converge, and the most probable state, assuming low enough temperatures, may be one of the negative energy discrete states.
Choosing the energy scale to be $Ry$ (Rydberg energy): $$ E_n=-\frac1{n^2} $$ and the degeneracy is $n^2$, so: $$ Z = \sum_{n=1}^\infty n^2e^{\beta/n^2} $$ The issue is that the density of states scales as: $$ D(E)\sim \frac1{2(-E)^{5/2}} $$ for $E\to0^-$. Indeed, the degeneracy of $E_n$ is $n^2$, so the number of states for $n<N$ scales as $\frac{N^3}3$. The scaling is exact in the classical limit, where it follows from dimensional analysis. The partition function of the boundary states therefore diverges at $0$: $$ Z = \int_{-\infty}^0e^{-\beta E}D(E)dE $$ You'll need to regularise it by imposing a cutoff energy $E_c$. Alternatively, you can physically change the potential say by imposing a finite volume, i.e. a cutoff radius $r_c$. Actually, both are equivalent as $E_c\sim -\frac1{r_c}$. Using the cutoff, you regularise the partition function: $$ Z \sim \frac{e^{\beta E_c}}{3(-E_c)^{3/2}} $$ which diverges as expected when $E_c\to0$.
Another way to get the result is to notice that formally, the second derivative of $Z$ converges: $$ \frac{d^2Z}{d\beta^2} = \sum_{n=1}^\infty \frac{e^{\beta/n^2}}{n^2} := Z_2 $$ You can therefore formally define $Z$ as the double antiderivative of $Z_2$, which makes sense since $Z_2$ is continuous. This defines $Z$ up to an undetermined affine term $A+B\beta$.
You can recover these constants by the previous regularisation approach, and check the consistency of the two methods. For example, you can introduce a cutoff quantum number $N$ (which corresponds to $N\sim(-E_c)^{-1/2}$), you get when $N\to\infty$: $$ A = \frac{N^3}3 \\ B = N $$ But without needing regularisation, you can view $A,B$ as adjustable parameters of your theory. While it may not seem predictive, it still severely restricts how mean energy, entropy, heat capacity or free energy depend on temperature for example. In particular, the low temperature behaviour does not depend on these constants, as you can check that $Z\sim e^\beta$.
Hope this helps.
