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}
$$