Using Legendre polynomials for the Gaussian quadrature is not very efficient. There are several reasons for this:

• You can easily end up in situations where the integrand diverges
• The limits $\pm \infty$ have to be approximated with a finite number

It is very useful here to change to spherical coordinates $$d{\bf r}_1d{\bf r}_2 = r_1^2dr_1 r_2^2dr_2 dcos(\theta_1)dcos(\theta_2)d\phi_1d\phi_2,$$ and $$\frac{1}{r_{12}}= \frac{1}{\sqrt{r_1^2+r_2^2-2r_1r_2cos(\beta)}}$$ with $$\cos(\beta) = \cos(\theta_1)\cos(\theta_2)+\sin(\theta_1)\sin(\theta_2)\cos(\phi_1-\phi_2))$$