The integral of $1/r_{12}$ in the cartesian space. So I have two electrically charged point particles $p$ and $p'$, with positions ${\vec x} = (x,y,z)$ and ${\vec x}'=(x',y',z')$ respectively, each with a charge $-1$. Assume we are working on Hartree atomic units such that Coulomb's constant equals $1$. That means that the Coulombic potential $v(p,p')$ between the two particles is given by: 
$$
v(p,p') = \frac{1}{|{\vec x} - {\vec x}'|}
$$
Where the notation $|\vec v|$ denotes the norm of vector $\vec v$. I tried finding closed-form for the following definite integrals.
$$
\int_{x_0}^{x_f} \int_{y_0}^{y_f} \int_{z_0}^{z_f} \frac{1}{|{\vec x}-{\vec x}'|} dx'dy'dz'
$$
$$
\int_{x_0}^{x_f} \int_{y_0}^{y_f} \int_{z_0}^{z_f} \int_{x_0}^{x_f} \int_{y_0}^{y_f} \int_{z_0}^{z_f} \frac{1}{|{\vec x}-{\vec x}'|} dx'dy'dz'  dxdydz
$$
Where the integration limits $x_0,x_f,y_0,y_f,z_0,z_f$ are the same for both particles, meaning that they are constrained to the same region of the 3D space. All of the variables and integration limits are elements of the set $\mathbb{R}$.
This seems like a pair of very complicated integrals, and I was ready for the possibility of not finding any closed-forms for them. However, the symmetry of the potential and the fact that these integrals probably appear so often in physics the the mathematics around them is probably very well developed made me look a bit harder. 
I fired the Mathematica and tried integrating the first one but the symbolic computation came out empty-handed. I tried integrating only for $z'$, but again the result was a conditional answer, with very weird convergence criteria involving various things about the real part of some algebraic function of the terms. I tried looking it up in physics and mathematical physics textbooks but to no avail. Finally I come here to settle the matter once and for all.
So the question is: Do the definite integrals above have closed forms? If so, what are they?
Schwinger Parametrization trial [edit]
So using the Schwinger Parametrization I wrote my solution $S$ like this:
$$
S=\frac{1}{\Gamma[1/2]} \int_0^\infty \lambda^{-1/2} \int_{x_0}^{x_f} \int_{y_0}^{y_f} \int_{z_0}^{z_f} \exp(-\lambda |\vec x - \vec x'|^2) \rm d \vec x ' \rm d \lambda
$$
But Mathematica did not manage to integrate it, again. I started trying it out by simplifying, by hand, as much as I could to see if the system managed to work its way to a solution. It turns out that the integral 
$$
\int_{x_0}^{x_f} \int_{y_0}^{y_f} \int_{z_0}^{z_f} \exp(-\lambda |\vec x - \vec x'|^2) \rm d \vec x '
$$
Is separable into three integrals of the type:
$$
\mathcal{I}_{x'} =\int_{x_0}^{x_f} \exp(-\lambda |x-x'|^2) \rm d x '
$$
With similar ones for y' and z'. So the final integration solution $S$ will be like this:
$$
S=\frac{1}{\Gamma[1/2]} \int_0^\infty \lambda^{-1/2} \mathcal{I}_{x'} \mathcal{I}_{y'} \mathcal{I}_{z'} \rm d \lambda
$$
The Integral $\mathcal{I}_{x'}$ has the form:
$$
\mathcal{I}_{x'}=-\frac{1}{2\sqrt{\lambda}} \sqrt{\pi} ({\rm Erf}[(x_0-x]\sqrt{\lambda})-{\rm Erf}[(x_f-x)\sqrt{\lambda}])
$$
When expanding the triple product $\mathcal{I}_{x'} \mathcal{I}_{y'} \mathcal{I}_{z'}$ we get a sum of a bunch of $S_i$ terms similar to term $S_0$ below.
$$
S_0=-\frac{\pi^{3/2}}{8 \lambda^{3/2}} (\rm Erf[(x_0 - x)\sqrt{\lambda}]\rm Erf[(y_0 - y)\sqrt{\lambda}]\rm Erf[(z_0 - z)\sqrt{\lambda}])
$$
For all combinations of integration limits $x_0,x_f,y_0,y_f,z_0,z_f$ and variables $x,y,z$ in their proper places, giving 8 terms in the total summation representation of $\mathcal{I}_{x'} \mathcal{I}_{y'} \mathcal{I}_{z'}$. Since 
$$
\mathcal{I}_{x'} \mathcal{I}_{y'} \mathcal{I}_{z'}=\sum_i S_i
$$
Then
$$
S=\frac{1}{\Gamma(1/2) }\sum_i \int_0^\infty \lambda^{-1/2} S_i \rm d \lambda
$$
I tried to integrate the term $S_0$ on Mathematica, as shown below, hoping that the result could help me solve the remaining $S_i$ terms and, finally, the complete integral $S$.
$$
\int_0^\infty \lambda^{-1/2} S_0 \rm d \lambda =\int_0^\infty -\lambda^{-1/2}\frac{\pi^{3/2}}{8 \lambda^{3/2}} (\rm Erf[(x_0 - x)\sqrt{\lambda}]\rm Erf[(y_0 - y)\sqrt{\lambda}]\rm Erf[(z_0 - z)\sqrt{\lambda}]) \rm d \lambda
$$
But Mathematica did not find out a way to resolve the above integral for $S_0$, returning to me empty-handed. Does it make sense that the system didn't manage to integrate that?
 A: Too long for a comment:
Have you considered any integral transformation first, for instance, the Fourier transform?
Why do I think this could help? Well, the dipole-dipole energy of a Bose-Einstein condensate can be calculated as
$$ E=\int_{-\infty}^{\infty}\mathrm{d}^3x\int_{-\infty}^{\infty}\mathrm{d}^3 x_1\, n(\vec{x})V_{dd}(\vec{x}-\vec{x}_1)n(\vec{x}_1) \tag{1}$$
Aplying Fourier transfom this reduces to 
$$E= \int_{-\infty}^{\infty}\frac{\mathrm{d}^3k}{(2\pi)^3}\, \tilde{n}(\vec{k})\tilde{V}_{dd}(\vec{k})\tilde{n}(-\vec{k}),$$
This last integral is simpler than the first.
Edit:
Here is how I did the Fourier Transform of Eq. $ (1) $. Maybe this can be useful.
(I did this a while ago, okay?)
Using the definition of Fourier transform:
\begin{align}
E&=\int_{-\infty}^{\infty}d^3x\int_{-\infty}^{\infty} d^3x_1\int_{-\infty}^{\infty}\frac{d^3k}{(2\pi)^3}\tilde{n}(\vec{k})e^{i\vec{k} \cdot \vec{x}}\int_{-\infty}^{\infty}\frac{d^3k_1}{(2\pi)^3}V_{dd}(\vec{k}_1)e^{i\vec{k}_1 \cdot (\vec{x}-\vec{x}_1)}\int_{-\infty}^{\infty}\frac{d^3k_2}{(2\pi)^3}\tilde{n}(\vec{k}_2)e^{i\vec{k}_2 \cdot \vec{x}_1}\\ &=\underbrace{\int_{-\infty}^{\infty}d^3x\,e^{i\vec{x}\cdot(\vec{k}+\vec{k}_1)}}_{{\equiv (2\pi)^3\delta(\vec{k}+\vec{k}_1)}}\underbrace{\int_{-\infty}^{\infty} d^3x_1\,e^{i\vec{x}_1\cdot(\vec{k}_2-\vec{k}_1)}}_{{\equiv (2\pi)^3\delta(\vec{k}_2-\vec{k}_1)}}\int_{-\infty}^{\infty}\frac{d^3k}{(2\pi)^3}\tilde{n}(\vec{k})\int_{-\infty}^{\infty}\frac{d^3k_1}{(2\pi)^3}V_{dd}(\vec{k}_1)\int_{-\infty}^{\infty}\frac{d^3k_2}{(2\pi)^3}\tilde{n}(\vec{k}_2) \\ &= \int_{-\infty}^{\infty}\frac{d^3k}{(2\pi)^3}\tilde{n}(\vec{k})\int_{-\infty}^{\infty}\frac{d^3k_1}{(2\pi)^3}V_{dd}(\vec{k}_1)(2\pi)^3\delta(\vec{k}+\vec{k}_1)\int_{-\infty}^{\infty}\frac{d^3k_2}{(2\pi)^3}\tilde{n}(\vec{k}_2)(2\pi)^3\delta(\vec{k}_2-\vec{k}_1) \\ &= \int_{-\infty}^{\infty}\frac{\mathrm{d}^3k}{(2\pi)^3}\, \tilde{n}(\vec{k})\tilde{V}_{dd}(\vec{k})\tilde{n}(-\vec{k}).
\end{align}
