Numerically integrating the Hartree potential

Question

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!

Answer 1

Based on your clarifications starting here, I see that you are starting with the simple case in which $n(r^\prime) = 1$. In this case we have the integral:

$$\tag{1} \int \frac{1}{|r-r^\prime|}\textrm{d}r^\prime, $$

for which can check to see whether or not our numerical integration procedure works by comparing with the analytic expression for the integral:

To do the integral in Eq. 1 numerically, we need to know how to numerically integrate a function of two variables, $r$ and $r^\prime$, with respect to one of them $r^\prime$. This search for "numerical integration of two variable function with respect to one variable" results in many immediate results, including (in order) on the following high-quality Q/A venues:

• Matlab Answers: "Integration of function with two variables with respect to one of them"
  • How to implement the trapezoid rule in MATLAB is the same as what you would want to do in NumPy.
• Mathematica StackExchange: "Numerical integration of two variable function"
  • Although I admit that Mathematica's NIntegrate function mentioned in the accepted answer, shelters you from getting to see the details, but the details for how to do it in MATLAB are indeed given in the question).
• StackOverflow: "numerical integration of two variable function over one variable only in C++ (using Numerical Recipes libraries)"
• ResearchGate: "What are the methods for the numerical integration of a function of more than one variable?

Basically what you are hoping to do is perfectly reasonable and has been implemented in many places, many times. You are on the right track, but as for your question:

"but won't the denominator just be zero?"

you may just have allowed $r = r^\prime$ at some point, which would have lead you to receive an error suggesting that your denominator is zero.

If you look at the analytic expression for the integral, in my above screenshot, you'll see that even the analytic expression becomes ugly when $r=r^\prime$, particularly because you'd be evaluating $\log(0)$.

Notice also that the analytic expression (in the above screenshot) for the integral, looks remarkably similar to what you would get if you calculate the integral of a related 1-variable function:

\begin{align}\tag{2} u &= r - r^\prime \\ \textrm{d}u &= -\textrm{d}r^\prime\\ -\int &\frac{1}{|u|}\textrm{d}u. \end{align}

However, while this last piece of information may be helpful in the early stages of your code development, in which you want to test your code for simple functions like $n(r^\prime)$ = 1, the more general numerical techniques for calculating the integral of a 2-variable function with respect to one of the variables, may be the best way forward for the general case.

Answer 2

The $r^{-1}$ operator is not bounded in 1D, as Nike's answer above illustrates: the integral diverges logarithmically. Indeed, it is well-known in the literature that even if you solve this by truncating 1/r to a constant beyind some small r value, your answer depends pathologically on the truncation threshold and does not converge to an answer in the limit of no truncation i.e. the constant going to infinity.

Thus, your physical model is wrong. You need to use a softened Coulomb operator: $$ \frac 1 r \to \frac 1 {\sqrt{r^2 + d^2}}$$ with $d$ a parameter.