Integrating a periodic, Gaussian-esque distribution I came across the following integral in this paper as part of my research and would like to understand how its value is found.
Let $\mathbf{k},\mathbf{k}' \in \mathbb{R}^3$, $\mathbf{G} \in \mathbb{Z}^3$, and $\alpha \in \mathbb{R}$, with the squared magnitudes of both $\mathbf{k}$ and $\mathbf{k}'$ strictly smaller than 1 (we say that $\mathbf{k}$ and $\mathbf{k}'$ are in the irreducible Brillouin zone, or BZ for short). Then
$$
\displaystyle\int\limits_{BZ} \displaystyle\sum_{\mathbf{G}} \frac{e^{-\alpha |\mathbf{k} - \mathbf{k}' - \mathbf{G}|^2}}{|\mathbf{k} - \mathbf{k}' - \mathbf{G}|^2}\, d\mathbf{k}' = 2\pi \sqrt{\frac{\pi}{\alpha}}.
$$
Why is this so? I can see the connection to the Gaussian integral (whence the factor of $\sqrt{\pi/\alpha}$), and the paper states that this function is periodic, which is why we can essentially ignore $\mathbf{k}$. But the infinite sum over $\mathbf{G}$-vectors is beyond my ken.
The context is condensed-phase theoretical chemistry; I'm working with Hartree–Fock exchange-style integrals in the plane-wave basis. The paper also references an earlier work, but that doesn't derive the integral either.
 A: The answer lies in seeing how the integration and summation combine. Since the integration connects the G terms they combine to be an integral over all space. This is more obvious if you switch the sum and integral order to
$$ \sum_{\boldsymbol{G}} \int_{BZ} d\boldsymbol{k}' \,f(\boldsymbol{k}-\boldsymbol{k}'-\boldsymbol{G});\ f(\boldsymbol{x})=\frac{e^{-\alpha \vert x\vert^2}}{\vert x\vert^2} $$
Switching the order shouldn't be a problem due to Fubini/Tonelli theorems since your function $f(\boldsymbol{x}) = \vert f(\boldsymbol{x}) \vert $. Then recognizing that for a given $\boldsymbol{k}$, the BZ integration connects a $\boldsymbol{G}=\langle n_x,n_y,n_z\rangle$ to the $ \langle n_x+1,n_y,n_z\rangle $, $ \langle n_x,n_y+1,n_z\rangle $, and $ \langle n_x,n_y,n_z+1\rangle $ vectors. This means the combination of the sum and integral becomes
$$ \sum_{\boldsymbol{G}} \int_{BZ} d\boldsymbol{k}' \,f(\boldsymbol{k}-\boldsymbol{k}'-\boldsymbol{G})= \int_{\boldsymbol{G}\times BZ} d(\boldsymbol{k}',\boldsymbol{G}) \,f(\boldsymbol{k}-\boldsymbol{k}'-\boldsymbol{G}) = \int_{\mathbb{R}^3} d\boldsymbol{x} f(\boldsymbol{x}-\boldsymbol{k})$$
Then since the integral is over all space the offset $\boldsymbol{k}$ doesn't matter and we get 
$$\int_{\mathbb{R}^3} d\boldsymbol{x} f(\boldsymbol{x}-\boldsymbol{k}) = \int_{\mathbb{R}^3} d\boldsymbol{x} f(\boldsymbol{x}) = \int_{\mathbb{R}^3} d\boldsymbol{x} \frac{e^{-\alpha \vert x\vert^2}}{\vert x\vert^2}, $$
which can be integrated in spherical coordinates.
$$ \int_0^{2\pi} d\phi  \int_0^\pi d\theta\ \text{sin}(\theta) \int_{0}^{\infty} dr  \, r^2 \frac{e^{-\alpha r^2}}{r^2} = 4\pi \int_{0}^{\infty} dr \ e^{-\alpha r^2} = 2\pi \sqrt{\frac{\pi}{\alpha}}$$
