I'm trying to write a 1D density functional theory code on Python to predict the electronic structure of certain ions using the local density approximation.
I'm trying to construct the hamiltonian for this system:
$$E \left[ n(r) \right] = -\frac 12 \frac{d^2}{dr^2} + V\_H(r) + V\_N(r) + V\_x(r)\tag{1},$$

where: $$V\_H(r) = \int \frac{ n( r')}{| r - r'|} dr'\tag{2}.$$

I've discretised my space and created matrix versions for each of the components of the energy functional but I'm not able to grasp how the $V\_H$ would look like as a matrix. To my understanding, it won't be diagonal but hermitian. I'm trying to use the trapezium method for numerical integration:

```
sum( n * delta_r_prime /np.sqrt(r - r_prime) )

```

but won't the denominator just be zero? I'm trying to understand what the matrix elements of $V\_H$ are and how to create this in Python but I'm getting really confused. I'd really appreciate it if someone could enlighten me on this!