Electrostatic energy integral for point charges The electric energy stored in a system of two point charges $Q_1$ and $Q_2$ is simply $$W = \frac{1}{4\pi\epsilon_0}\frac{Q_1Q_2}{a}$$ where $a$ is the distance between them.
However, the total energy can also be calculated through the volume integral of magnitude squared of the electric over all space: $$W = \frac{\epsilon_0}{2}\int_{\mathbb{R}^3} E^2\,dV$$
Suppose that $Q_1$ sits on the origin and $Q_2$ a distance $a$ away on the $z$-axis.  Then the electric field is $$\vec{E}(x,y,z) = \frac{1}{4\pi\epsilon_0}\left[\left(\frac{Q_1}{r^3} + \frac{Q_2}{d^3}\right)x\hat{x} +\left(\frac{Q_1}{r^3} + \frac{Q_2}{d^3}\right)y\hat{y} + \left(\frac{Q_1z}{r^3} + \frac{Q_2\left(z-a\right)}{d^3}\right)\hat{z}\right]$$
where
\begin{align}
r &= \sqrt{x^2+y^2+z^2} \\
d&= \sqrt{x^2+y^2+\left(z-a\right)^2} 
\end{align}
Calculating $W$ through $\int_{\mathbb{R}^3} E^2\,dV$ seems extremely difficult; Mathematica, for example, appears stumped.  Yet its result should simply be $\frac{1}{4\pi\epsilon_0}\frac{Q_1Q_2}{a}$, correct?  The integral formula still applies to point charges, correct?
 A: Let us call the two electric charges $Q$ and $q$ with electric fields $|{\bf E}_Q|=\frac{k_e |Q|}{r^2_Q}$ and $|{\bf E}_q|=\frac{k_e |q|}{r^2_q}$, respectively. Here $k_e=\frac{1}{4\pi\varepsilon_0}$. The energy 
$$\tag{0} U~=~\frac{\varepsilon_0}{2}\iiint_{\mathbb{R}^3}\! d^3r ~|{\bf E}_Q+{\bf E}_q|^2~=~U_Q+U_q+U_{Qq} $$
of the total electric field contains three contributions: 


*

*The energy $U_Q$ of the electric field of a charge $Q$ 
$$\tag{1} U_Q~=~\frac{\varepsilon_0}{2}\iiint_{r\geq\delta}\! d^3r ~|{\bf E}_Q|^2~=~ \frac{k_eQ^2}{2} \int_{\delta}^{\infty}\frac{dr}{r^2}~=~\frac{k_eQ^2}{2\delta}.$$

*The energy $U_q$ of the electric field of a charge $q$ 
$$\tag{2} U_q~=~\frac{\varepsilon_0}{2}\iiint_{r\geq\delta}\! d^3r ~|{\bf E}_q|^2~=~ \frac{k_eq^2}{2} \int_{\delta}^{\infty}\frac{dr}{r^2}~=~\frac{k_eq^2}{2\delta}.$$
In equation (1) and (2), we have inserted a regulator $\delta$. If the regulator $\delta\to 0$ is removed the energy becomes infinite.

*The energy from the exchange term
$$
\tag{3} U_{Qq}~=~\varepsilon_0\iiint_{\mathbb{R}^3}\! d^3r ~{\bf E}_Q\cdot {\bf E}_q~=~\frac{k_eQq}{R},$$
where $R$ is the distance between the two charges. The triple integral (3) can be analytically calculated by (among other things) using the azimuthal symmetry of the integrand.
If we slowly vary the separation $R$, we will only detect the variation of the exchange term (3) (the Coulomb energy) as the two other contributions (1) and (2) remain fixed.
A: Point charges pose a problem because their self energy $\rightarrow \infty$ as their size$\rightarrow 0$
It represents the work done in assembling a point charge, whereas you're only interested in the work done in moving charges in the fields of other charges. You therefore need to subtract the energy stored in the electric fields of the isolated points charges from the total. If you do it right, you should end up with an expression for the energy that depends only on their separation from one another which won't blow up when evaluated.
A: Yes, the integral formula applies to point charges and yes, your second integral does describe the energy. That version, however it not the easiest way to calculate the energy. It is a good way to calculate the energy when E = 0 in a large region.
E For point charges:
$$\mathbf E (\mathbf r) = \sum_i \frac{q_i}{4 \pi \varepsilon_0 R^2} \mathbf {\hat R}$$
$\mathbf{\hat R}= \frac{\mathbf r - \mathbf r'}{|\mathbf r - \mathbf r'|}$ where $\mathbf r'$ is where the charge is and $\mathbf r $ is the point you are looking at.
