Another approach worth trying may be using the convolution property of the Fourier Transform.
The 3-D cartesian Fourier Transform:
$$F(x, y, z) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x,y,z) e^{-2\pi i (xu + yv + zw)} \space dx \space dy \space dz$$
when there is 3-D spherical symmetry, can be expressed as:
$$F(q) = 4\pi \int_0^\infty f(r)r^2\mathrm{sinc}(2qr)dr$$
with the inverse transform expressed as:
$$f(r) = 4\pi \int_0^\infty F(q)q^2\mathrm{sinc}(2rq)dq$$
with
$$\mathrm{sinc}(x) \equiv \dfrac{\sin(\pi x)}{\pi x}$$
So starting with the case $n = 0$
$$f(r) = e^{-\frac{3}{4\pi a^2}\pi r^2} \quad \text{so} \quad F(q) = \sqrt{\dfrac{4\pi a^2}{3}}e^{-\frac{4\pi a^2}{3}\pi q^2}$$
$$g(r) = Ae^{-br} \quad \text{so} \quad G(q) = A\dfrac{8\pi b}{\left[b^2+(2\pi q)^2\right]^2}$$
(For $n > 0$, a derivation of a derivative property for the 3-D Fourier Transform will probably be needed to get $G(q)$.)
So now to attempt the convolution for the $n= 0$ case:
$$\begin{align*}(f *g)(r) &= \mathscr{F}^{-1}\left\{F(q)G(q)\right\} \\
\\
&= 4\pi \int_0^\infty \sqrt{\dfrac{4\pi a^2}{3}}e^{-\frac{4\pi a^2}{3}\pi q^2} A\dfrac{8\pi b}{\left[b^2+(2\pi q)^2\right]^2} q^2 \mathrm{sinc}(2rq) dq\\
\\
&= \dfrac{32\pi Aab}{r}\sqrt{\dfrac{\pi}{3}}\int_0^\infty e^{-\frac{4\pi a^2}{3}\pi q^2} \dfrac{1}{\left[b^2+(2\pi q)^2\right]^2} q \sin(2\pi rq) dq
\end{align*}$$
And that integral looks a little more tractable than the one in the original question. However it is still not obvious if one can get a closed form solution using substitutions and integration by parts.