How to evaluate this monsterous integral? How to prove this?
$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}  \frac {2(x-x^2)e^{-(x^2+y^2+z^2)}}{((x-1)^2+(y-1)^2+(z-1)^2)^{3/2}}dxdydz=\frac {-4\pi e^{-3}}{3}$$
I tried converting to polar coordinates but still I couldn't evaluate.
 A: Another answer:
I noticed it also follows from this identity see identity
$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}  \frac {\partial f(x',y',z')}{\partial x'}\frac {x'-x}{4\pi r'^3}dx'dy'dz'+\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}  \frac {\partial f(x',y',z')}{\partial y'}\frac {y'-y}{4\pi r'^3}dx'dy'dz'+
\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}  \frac {\partial f(x',y',z')}{\partial z'}\frac {z'-z}{4\pi r'^3}dx'dy'dz'=-f(x,y,z)$$
with $$f(x,y,z)=e^{-(x^2+y^2+z^2)}$$
it is easily seen when $x=y=z=1$ $\quad$ $f(1,1,1)=e^{-3}$
$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}  \frac {\partial f(x',y',z')}{\partial x'}\frac {x'-1}{4\pi r'^3}dx'dy'dz'=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}\frac {2(x'-x'^2)e^{-(x'^2+y'^2+z'^2)}}{4\pi((x'-1)^2+(y'-1)^2+(z'-1)^2)^{3/2}}dx'dy'dz'=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}\frac {2(y'-y'^2)e^{-(x'^2+y'^2+z'^2)}}{4\pi((x'-1)^2+(y'-1)^2+(z'-1)^2)^{3/2}}dx'dy'dz'=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}\frac {2(z'-z'^2)e^{-(x'^2+y'^2+z'^2)}}{4\pi((x'-1)^2+(y'-1)^2+(z'-1)^2)^{3/2}}dx'dy'dz'$$
Therefore: $$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}\frac {2(x'-x'^2)e^{-(x'^2+y'^2+z'^2)}}{4\pi((x'-1)^2+(y'-1)^2+(z'-1)^2)^{3/2}}dx'dy'dz'=\frac {-f(1,1,1)}{3}$$
$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}\frac {2(x'-x'^2)e^{-(x'^2+y'^2+z'^2)}}{((x'-1)^2+(y'-1)^2+(z'-1)^2)^{3/2}}dx'dy'dz'=\frac {-4\pi e^{-3}}{3}$$
A: Let $r = \sqrt{x^2+y^2+z^2}$ and $r_1 = \sqrt{(x-1)^2+(y-1)^2 + (z-1)^2}$. 
For any $R, \epsilon > 0$, let


*

*$B_0(R) = \{ (x,y,z) : r < R \}$ be the ball of radius $R$ centered at origin.

*$B_1(\epsilon) = \{ (x,y,z) : r_1 < \epsilon \}$ be the ball of radius $\epsilon$ centered at $(1,1,1)$.


Notice the integrand blows up as $\frac{1}{r_1^2}$ as $(x,y,z)$ approaches $(1,1,1)$ which get suppressed by a factor $r_1^2$ from the volume element. 
The triple integral should be interpreted as a limit of following form:
$$\iiint ( \cdots ) = \lim_{R\to\infty}\lim_{\epsilon\to 0}\int_{B_0(R)\setminus B_1(\epsilon)} (\cdots )
$$
By symmetry, the integral can be rewritten as
$$\begin{align}
  & -\frac23\iiint \frac{x(x-1)+y(y-1)+z(z-1)}{r_1^3} e^{-r^2} dV\\
= & -\frac23 \iiint \left(\vec{x}\cdot\frac{\vec{\nabla}r_1}{r_1^2}\right) e^{-r^2} dV\\
= & -\frac13 \iiint \vec{\nabla}\left(\frac{1}{r_1}\right)\cdot\vec{\nabla}e^{-r^2} dV\\
= & -\frac13 \iiint 
\left[
\vec{\nabla}\cdot\left(e^{-r^2}\vec{\nabla}\left(\frac{1}{r_1}\right)\right)
- e^{-r^2}\nabla^2\left(\frac{1}{r_1}\right)
\right]
dV\\
\end{align}
$$
Notice $\nabla^2\left(\frac{1}{r_1}\right) = 0$ outside $B_1(\epsilon)$, we can use divergence theorem to rewrite the integral as
$$-\frac13\lim_{R \to \infty}\lim_{\epsilon \to 0}
\int_{\partial (B_0(R)\setminus B_1(\epsilon))} e^{-r^2}\vec{\nabla}\left(\frac{1}{r_1}\right)\cdot d\vec{S}
= -\frac13\left[
\lim_{R \to \infty} \int_{\partial B_0(R)}
-
\lim_{\epsilon \to 0} \int_{\partial B_1(\epsilon)}
\right]
\left( e^{-r^2}\vec{\nabla}\left(\frac{1}{r_1}\right)\cdot d\vec{S} \right)
$$
Because of the $e^{-r^2}$ factor, the integral over $\partial B_0(R)$ vanishes as
$R \to \infty$. 
For the remaining integral, since $e^{-r^2}$ is continuous near $(1,1,1)$, we can pull it outside the integral sign. As a result, the integral
we seek is equal to
$$\frac13
\left(\lim_{(x,y,z)\to (1,1,1)} e^{-r^2}\right)
\left(\lim_{\epsilon \to 0}\int_{\partial B_1(\epsilon)}\vec{\nabla}\left(\frac{1}{r_1}\right)\cdot d\vec{S}\right)
= \frac13 e^{-3} (-4\pi) = -\frac{4\pi}{3e^3}
$$
