My question is what are the steps for taking an integral of the following form?
$$\int e^{-\alpha|\mathbf r- \mathbf R_a|^2} {1\over|\mathbf r- \mathbf R_b|} e^{-\beta|\mathbf r- \mathbf R_b|^2} dV$$
This integral is commonly seen when attempting to do Quantum Chemistry calculations with a Gaussian type basis set. I have tried to use wolfram alpha to solve this problem but it fails to give a solution.
Since the integral is over all space, it makes sense to shift the origin of your coordinate system to $\mathbf{R}_b$. Then your integral becomes $$I = \int e^{-\alpha|\mathbf{r}'-\Delta\mathbf{R}|^2}\frac{1}{|\mathbf{r}'|}e^{-\beta|\mathbf{r}'|^2}dV',$$ where $\Delta\mathbf{R}\equiv\mathbf{R}_a-\mathbf{R}_b$.
The next trick to use would be to choose to align the $z'$ axis along the direction of $\Delta\mathbf{R}$. Then you'll have that $$I = 2\pi\int_0^\infty r'^2dr'\int_0^\pi\sin\theta' d\theta' \, e^{-\alpha(r'^2-2r'\Delta R\cos\theta'+\Delta R^2)}\frac{1}{r'}e^{-\beta r'^2},$$ where I've already performed the integration over $d\phi'$. The above integral should be do-able by hand (or by Wolfram Alpha).
Let's call this integral $C = C({\bf R}_a,{\bf R}_b,\alpha,\beta)$, since the arguments are the free variables in the integral. Let's also assume that we're working in $\mathbb R^3$.
Our first attack consists of performing the substitution ${\bf r}\rightarrow{\bf r}+{\bf R}_b$. This won't change the volume element, so the result is now: $$ C = \int d^3{\bf r}\ \ e^{-\alpha|\mathbf r- \Delta \mathbf R|^2} {1\over r} e^{-\beta r^2} $$ where $r = |{\bf r}|$ and we have defined $\Delta{\bf R}\equiv {\bf R}_a-{\bf R}_b$. That looks a lot simpler, doesn't it? We see that the integral only depends on the separation between a and b. $C = C(\Delta{\bf R},\alpha,\beta)$ But there's that pesky magnitude in the first exponential...
We can handle this if we work in spherical coordinates, so that $d^3{\bf r}=r^2\sin\theta \ dr\ d\theta\ d\phi$. Don't forget the Jacobian factor in the volume element or you will get 7 years of bad luck. Now, recall that $$ \begin{align} |{\bf r}-\Delta{\bf R}|^2&=({\bf r}-\Delta{\bf R})\cdot({\bf r}-\Delta{\bf R})\\ &=r^2+\Delta R^2-2r\Delta R\cos\theta \end{align} $$ If we plug this into the integral, along with the volume element, use the substitution $u=\cos\theta$, and integrate over $\phi$, then we get: $$ C = 2\pi e^{-\alpha\Delta R^2}\int_0^\infty dr \ re^{-(\alpha+\beta)r^2}\int_{-1}^1du \ e^{2\alpha r\Delta R u} $$ The $u$ integral is pretty simple, leaving us with $$ C = \frac{\pi}{a}\frac{e^{-\alpha\Delta R^2}}{\Delta R}\int_0^\infty dr\ \left(e^{-(\alpha+\beta)r^2+2\alpha r\Delta R}-e^{-(\alpha+\beta)r^2-2\alpha r\Delta R}\right) $$ Now we have two incomplete Gaussian integrals, so the results will introduce error functions. Note that in both integrals we can complete the square and write them as: $$ \begin{align} e^{\frac{\alpha^2}{\alpha+\beta}\Delta R^2}\int_0^\infty dr e^{-(\alpha+\beta)(r\pm \frac{\alpha\Delta r}{\alpha+\beta})^2}&=e^{\frac{\alpha^2}{\alpha+\beta}\Delta R^2}\int_{\pm \frac{\alpha\Delta r}{\alpha+\beta}}^\infty dr \ e^{-(\alpha+\beta)r^2}\\ &=\frac{1}{2}\sqrt\frac{\pi}{\alpha+\beta}\text{erfc}(\pm\frac{\alpha\Delta R}{\sqrt{\alpha+\beta}}) \end{align} $$ where we have introduced the complimentary error function. Recalling the identity $\text{erfc}(-x)-\text{erfc}(x)=2\text{erf}(x)$, we can plug everything in to find: $$ C(\Delta R,\alpha,\beta)=\frac{\pi^{3/2}e^{-\frac{\alpha\beta}{\alpha+\beta}\Delta R^2}}{\alpha\sqrt{\alpha+\beta}\Delta R}\text{erf}(\frac{\alpha\Delta R}{\sqrt{\alpha+\beta}}) $$ One can quickly check that the units for both answers are $Length^2$, as they should be. The argument of the error function must be dimensionless. All is well in the world.
Hope that helps!
Here is a more formal complete solution going off the approach suggested from the comments.
$I=\int e^{-\alpha|\mathbf r-\mathbf R_a|^2} {1 \over |\mathbf r-\mathbf R_b|} e^{-\beta|\mathbf r-\mathbf R_b|^2}dV $
First changing the arrangement of the equation we can define $\mathbf \Delta \mathbf R$ and $\mathbf r'$as:
$\mathbf \Delta \mathbf R = \mathbf R_a - \mathbf R_b$
$\mathbf r'=\mathbf r - \mathbf R_b$
Substitution leads to:
$I=\int e^{-\alpha|\mathbf r'-\mathbf \Delta \mathbf R|^2} {1 \over |\mathbf r'|} e^{-\beta|\mathbf r'|^2}dV $
Aligning the $z'$ axis along the direction of $\mathbf \Delta \mathbf R$
$\mathbf r'=r'sin(\phi')cos(\theta')\hat x'+r'sin(\phi')sin(\theta')\hat y'+r'cos(\phi')\hat z'$
$\mathbf \Delta \mathbf R=0\hat x' + 0\hat y' + \Delta R\hat z'$
$|\mathbf r' -\mathbf \Delta \mathbf R|=\sqrt {(r'sin(\phi')cos(\theta')\hat x')^2+(r'sin(\phi')sin(\theta')\hat y')^2+(r'cos(\phi')\hat z'-\Delta R\hat z')^2}$
Returning to the integral in spherical coordinates and plugging in:
$I=\int^\infty_0 r' \int^{\pi}_0 sin(\phi')\int^{2\pi}_0 e^{-\alpha((r'sin(\phi')cos(\theta')\hat x')^2+(r'sin(\phi')sin(\theta')\hat y')^2+(r'cos(\phi')\hat z'-\Delta R\hat z')^2)} e^{-\beta r'^2} d\theta' d\phi' dr' $
This leads to
$I=2\pi\int^\infty_0 r' \int^{\pi}_0 sin(\phi') e^{-\alpha(r'^2-2r'\Delta R cos(\phi')+\Delta R^2)} e^{-\beta r'^2} d\phi' dr' $
$I=2\pi\int^\infty_0 e^{-\alpha(\Delta R+r')^2} {{e^{4\alpha \Delta R r'} -1} \over 2\alpha\Delta R} e^{-\beta r'^2} dr' $'
Finally, if I did everything correctly
$I = {\pi^{3 \over 2} e^{{-\alpha \beta\Delta R^2 \over \alpha+\beta}}Erf({\Delta R \alpha \over \sqrt{\alpha+\beta}}) \over {\Delta R \alpha \sqrt{\alpha+\beta}} }$'
