# Cannot solve this integral used in quantum chemistry

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?

• you are working with nuclei…are they spheric by any chance? (are your bounds in spherical coordinates?) – Kuifje Nov 26 '15 at 20:43
• Yeah good call, maybe you can simplify this in spherical coordinates, but still if an integral can be solved in spherical it can be solved in cartesian right? So this should still be solvable. – user2879934 Nov 26 '15 at 20:45
• no not necessarily, for example $\iint e^{x^2+y^2}dxdy$ is trivial in polar coordinates but very hard if not impossible in cartesian. Since your function has terms in $x^2+y^2+z^2=\rho^2$ I would definitely try it in spherical coordinates. – Kuifje Nov 26 '15 at 20:48
• The formula for $r_{12}$ looks wrong. What is $(x,y,z)$, location of what? – achille hui Nov 26 '15 at 20:49
• @Kuifje could you maybe give some hints on how to do that... sorry haven't looked at integrals since high school :) – user2879934 Nov 26 '15 at 20:50

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.