Recently I learned that electrostatic potential energy of a system of charges can be calculated like so:
$$E = \int \frac{1}{2} \epsilon_0 \mathbf E^2 dV.$$
This is correct for regular distributions of charge (finite everywhere) and for some singular distributions as well (charge on a surface). It is not valid for point particles, because (as is easily seen) the integral diverges, and more importantly, because the derivation breaks down for point particles. Instead, one should use the original expression due to Coulomb
$$
E = \sum_{i=1}^N \sum_{k=1}^N{}^{'}\frac{1}{2} \frac{Kq_iq_k}{|\mathbf r_i - \mathbf r_k|}
$$
which has no such problems (it does not diverge and can be derived using the Coulomb's law based on empiric observations). The prime next to sum sign means that the case $k=i$ is to be omitted (this also follows from the Coulomb law).
However, for a system of two point charges Q1 and Q2, we have
$$E = \int \frac{1}{2} \epsilon_0 (\mathbf E_1 + \mathbf E_2)^2 dV = \int_{}^{}\frac{1}{2}\epsilon_0 \mathbf E_1^2 dV + \int_{}^{}\frac{1}{2}\epsilon_0 \mathbf E_2^2 dV + \int_{}^{}\frac{1}{2}\epsilon_0 * 2 \mathbf E_1 \cdot \mathbf E_2 dV$$
from which we subtract the first two terms because they are already properties of the charges themselves.
This is a misguided reasoning for point charges, because energy of such system is not given by integral of total field squared. But it does result in the correct value: the integral
$$
\int \epsilon_0 \mathbf E_1 \cdot \mathbf E_2 dV
$$
has the same value as the Coulomb energy of two point charges.
Why do we not need to do something like this for charge distributions? Aren't charge distributions made up of many many little charges of magnitude $e$? Even if we treat the charge distribution as perfect continuous, won't there be infinitely infinitesimal charges?
For regular distributions of charge, energy is defined to be integral of electric field squared and usually, for finite integration volume, comes out finite. So no subtraction is needed.
For point-like distributions, energy cannot be validly defined to be integral of electric field squared, but it can be validly defined by the multilinear Coulomb expression above. So no subtraction is needed either.
Therefore, energy of point particles (which have singular spatial distribution of electric charge) is not, in general, a limit of energy of continuous distribution that converges to that singular distribution (the limit is infinite).
Another way to see this is to consider work of electric forces on a set of closely spaced point charges:
$$
W_{point} = \sum_k q_k \mathbf E_{-k} \cdot \mathbf v_k \Delta t
$$
where $\mathbf E_{-k}$ is electric field experienced by the particle $k$ and $\mathbf v_{k}$ is velocity of the $k-$th particle. In the regular distribution model, one would write the work happening on that system as (we denote the volume where the particles are $\Delta V$):
$$
W_{regular} = \mathbf E \cdot \mathbf j~ \Delta V\Delta t
$$
where $\mathbf E$ is average of $\mathbf E_{-k}$ and $\mathbf j \Delta V$ is average of $q_k\mathbf v_k$. But this will have a very different value from $W_{point}$, because product of averages is not in general equal to average of products. As an example, if the electric field experienced is the same for all particles, and all particles have the same velocity except one, which will have much greater magnitude, the value of work obtained will be greatly biased by the outlier particle and so will be much higher than what the point model would give.
The work of electric field manifests as decrease of energy of electric field. Consequently, although the point charge distribution is a limit of regular distribution, electric field energy of point distribution is not a limit of electric field energy of regular distributions of charge. These are two different models with different valuations of energy.
For more on the point charge model and how energy is defined there, see for example
J. Frenkel, Zur Elektrodynamik punktfoermiger Elektronen, Zeits. f. Phys., 32, (1925), p. 518-534. http://dx.doi.org/10.1007/BF01331692
In English, this article also explains it concisely:
R. C. Stabler, A Possible Modification of Classical Electrodynamics, Physics Letters, 8, 3, (1964), p. 185-187. http://dx.doi.org/10.1016/S0031-9163(64)91989-4