How to compute Coulomb repulsion from electron density numerically The energy due to mutual repulsion of electrons (Hartree potential) is usually formulated similar to
$$U_{ee}=\frac{e^2}{2}\int\frac{n(r)n(r')}{|r-r'|}\mathrm{d^3}r~ \mathrm{d^3}r'$$
I'm not a physicist and I don't know how to interpret this integral. How can I compute it numerically?
I've tried the following:
tmp=0:
for j=1:length(n)
    for k=1:length(n)
           if j==i
               continue
           end                  
           tmp=tmp+((n(j)*n(k))/abs(r(j)-r(k)));
    end
end
U_ee=tmp/2;

However, this results in a huge electrostatic energy which is more than $10^6$ times the kinetic energy at the same temperature and pressure.
I've also heard, that $U_{ee}$ can be calculated in a different way using the Poisson equation. 
Any hints?
 A: We have Maxwell's equation (I'll be using CGS, sorry for that)
$
\nabla\cdot\mathbf{E}=4\pi\rho,
$
where $E=-\nabla\phi$ (where $\phi$ is the scalar potential). Which gives us: $$\nabla^2\phi=-4\pi\rho.$$ 
This is the Poisson's equation and its general solution (given by the Green's function) is:
$$
\phi(\mathbf{r})=\int\frac{\rho(\mathbf{r}')}{\left|\mathbf{r}-\mathbf{r}'\right|}\mathrm{d}^3\mathbf{r}'.
$$
And hence the overall energy can be found as:
\begin{equation}
W=\frac{1}{2}\int\rho\phi\mathrm{d}^3\mathbf{r}.
\end{equation}
If you have a system of quasi-discrete charges, then $\rho(\mathbf{r}_i)=e\cdot n(\mathbf{r}_i)=e\cdot n_i$. If there are point charges (as far as I understood from your code), then you'll simply have a sum of $\delta$-functions, which will give you just a discrete sum (only over $\mathbf{r}_i$) like:
$$
\phi_i=\phi(\mathbf{r}_i)=\sum_{j\ne i}\frac{e\cdot n_j}{\left|\mathbf{r}_i-\mathbf{r}_j\right|},
$$
and for the energy we'll get:
$$
W=\frac{e}{2}\sum_i n_i\cdot\phi_i.
$$
So your code will look like this (in Python):
def phi(ind):
   phi0 = 0
   for j in range(NP):
      if ind != j:
         phi0 += e * n(j) / DIST(ind, j)
   return phi0
W = 0
for i in range(NP):
   W += 0.5 * e * n(i) * phi(i)

