Cannot solve this integral used in quantum chemistry

Question

I am writing computer code for an implementation of the Hartree Fock algorithm and I am stuck on a certain integral. This is a great walkthrough to get some background : HFTheory

Anyway, the set-up is this. I am using a single Guassian type orbital as my basis function for an approximation to helium. This is defined (in cartesian coordinates) by

$$\displaystyle \phi_2(x,y,z) = ne^{-a(x^2 + y^2 + z^2)}$$

And I want to plug that into the following

$$U_\mathrm{eff}(1)= \int \mathrm{d}(2)\,\frac{\left|\phi_2(x,y,z)\right|^2}{r_{12}}$$ This leads to something like so $$U_\mathrm{eff}(1)= \int \int \int \frac{\left|\phi_2(x,y,z)\right|^2}{r_{12}} dx dydz$$ $$U_\mathrm{eff}(1)= \int \int \int \frac{\left|ne^{-a(x^2 + y^2 + z^2)}\right|^2}{r_{12}} dx dydz$$ Now $r_{12}$ is the distance between the nuclei so this is just $\left|r_1 - \sqrt{x^2+y^2+z^2}\right|$ so plugging this in I get $$U_\mathrm{eff}(1)= \int \int \int \frac{\left|ne^{-a(x^2 + y^2 + z^2)}\right|^2}{\left|r_1 - \sqrt{x^2+y^2+z^2}\right|} dx dydz$$

I've plugged this into many a program and it seems like it is unsolvable, any ideas?

Answer 1

Choose the coordinate axes and introduce polar coordinates such that

• the $1^{st}$ electron is located at $\vec{p}_1 = (0,0,r_1)$.
• the $2^{nd}$ electron is located at $\vec{p}_2 = (x,y,z) = (r \sin\theta\cos\phi, r \sin\theta\sin\phi, r\cos\theta)$

It is known that $\frac{1}{r_{12}} = \frac{1}{|\vec{p}_1 - \vec{p}_2|} = \frac{1}{\sqrt{x^2+y^2+(r_1-z)^2}}$ has following Multipole expansion (see $\S 3.3$ of Classical Electrodynamics by J.D.Jackson ). $$\frac{1}{r_{12}} = \sum_{\ell=0}^\infty \frac{r_<^{\ell}}{r_>^{\ell+1}}P_\ell(\cos\theta)$$ where $r_< = \min(r_1, r)$, $r_> = \max(r_1,r)$ and $P_\ell(t)$ is the Legendre polynomial of degree $\ell$.

The sort of integral you want to calculate can be expressed as

$$\int \frac{e^{-\beta r^2}}{r_{12}} d^3\vec{p}_2 = \sum_{\ell=0}^\infty \left( \int_0^\infty e^{-\beta r^2} \frac{r_<^\ell}{r_>^{\ell+1}} r^2 dr \right)\left(\int_0^{\pi} P_\ell(\cos\theta) \sin\theta d\theta \right)\left(\int_0^{2\pi} d\phi\right)\tag{*1}$$

The key is

$$\int_0^{\pi} P_\ell(\cos\theta) \sin\theta d\theta = \int_{-1}^{1} P_\ell(t) dt = \begin{cases} 2, &\ell = 0,\\ 0, &\text{ otherwise } \end{cases} $$ and we only need to keep the $\ell = 0$ term in $(*1)$. As a result, $$\int \frac{e^{-\beta r^2}}{r_{12}} d^3\vec{p}_2 = 4\pi \int_0^\infty e^{-\beta r^2} \frac{r^2}{r_>} dr = 4\pi \left(\int_0^{r_1} e^{-\beta r^2} \frac{r^2}{r_1} dr + \int_{r_1}^\infty e^{-\beta r^2} r dr \right) $$ I will leave the rest of computation to you.