# Numerically integrating the Hartree potential

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!

• +1 and welcome to our new community! Thank you for contributing your question here, and we hope to see much more of you in the future !!! How did you make matrix versions o the other terms? Mar 19, 2022 at 17:56
• Thanks for your answer. I've moved it to chat and answered there :) Mar 19, 2022 at 18:10
• You already have a very nice answer to your question but I'm unsure whether it actually helps you. I guess with 1D you mean that you have a periodic situation in this dimension. In that case the integral would go from -infinity to +infinity and it is obvious that the result of the integral would also be inifinite. The point here is that in such a situation you cannot treat the Hartree potential separately from the other contributions to the potential. You will also see diverging integrals in those terms and you can only calculate the potential if you treat cancelling terms in a smart way. Mar 20, 2022 at 17:04

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.

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.

• +1. Perhaps (if you get enough time later) you can show the user some literature that shows those results using truncations? No rush though, take your time! Mar 20, 2022 at 23:08
• Looking up the literature is left as an exercise for the reader ;) Mar 21, 2022 at 0:46