Solution of the atom within LDA

In this excercise we will solve the multielectron atom in LDA approximation.

We will test it on He and oxygen by computing the total energy and charge density.

We will plot charge density and compute the total energy, which will be compared to the reference data at NIST database: https://www.nist.gov/pml/atomic-reference-data-electronic-structure-calculations/atomic-reference-data-electronic-7

We want to solve the Schroedinger equation for an atom with nucleous charge Z. We will approximate electron-electron interaction with an effective potential, which is computed by so-called ”local density approximation” (LDA). In this theory, the classical (Hartree) potential is treated exactly, while the rest of the interaction is ”hidden” into approximate exchange-correlation functional. We will not go into details of this functional, we will just use it here.

The Schroedinger equation we are solving is \begin{eqnarray} (-\frac{\hbar^2}{2m}\nabla^2-\frac{Z e^2}{4\pi\varepsilon_0 r} + V_H(r) + V_{xc}(r))\psi(\vec{r})=E \psi(\vec{r}) \end{eqnarray}

The first two terms are appearing in Hydrogen problem, and we already coded them. The Hartree is the electrostatic potential, and the exchange-correlation potential has an approximate form, which depends only the charge density of the system. We will use the module excor.py, where the function $V_{xc}(\rho)$ is tabulated. We will use it as $V_{xc}(r)== V_{xc}(\rho(r))$, corresponding to the local density approximation.

First we take the code from the Hydrogen project and repeat.

Hartree term

The Hartree term is treated exactly in this approximation. It describes the electrostatic interaction of one electron with the cloud of all electrons (including the electron itself). Mathematically, this term is \begin{eqnarray} && \frac{1}{2}\int d\vec{r} d\vec{r}' \psi^\dagger (\vec{r})\psi^\dagger (\vec{r}') v_c(\vec{r}-\vec{r}') \psi(\vec{r}')\psi(\vec{r}) \rightarrow\\ && \int d\vec{r} \psi^\dagger(\vec{r}) \psi(\vec{r}) \int d\vec{r}' \langle\psi^\dagger(\vec{r}') \psi(\vec{r}')\rangle v_c(\vec{r}-\vec{r}') \equiv \int d\vec{r} \psi^\dagger(\vec{r}) V_{H}(\vec{r}) \psi(\vec{r})\nonumber \end{eqnarray} with \begin{equation} V_H(\vec{r}) = 2 \int d\vec{r}' \frac{\rho(\vec{r}')}{|\vec{r}-\vec{r}'|} \end{equation} where $2$ is due to Rydberg units sinc $v_c = 2/r$.

For any atom, the electron density is spherically symetric and hence $V_{H}$ depends only on radial distance. (In solids, the hartree potential should be expanded in spheric harmonics to sufficiently high $l$, maybe $l=6$).

Step 1: Using $\rho(r)$ computed in previous homework, compute the Hartree potential.

This is usually achieved by solving the Poisson equation. From clasical electrostatic we know \begin{eqnarray} \nabla^2 V_{H}(\vec{r}) = -8\pi \rho(\vec{r}) \end{eqnarray} In Hartree approximation, we have \begin{equation} \frac{1}{r^2}\frac{d}{dr}(r^2 \frac{d V_H}{dr})= -8\pi\rho(r) \end{equation} which simplifies to \begin{equation} U^{''}(r) = -8\pi r \rho(r) \end{equation} where $U(r) = V_{H}(r) r$.

This second order differential equation has the following boundary conditions $U(0)=0$ and $U(\infty)=2 Z$.

The two point boundary problem does not require shooting because we know solution to the homogenous differential equation $U^{''}(r)=0$. The Hartree potential can be obtained from any particular solution by \begin{equation} U(r) = U_p(r) + \alpha r \end{equation} where $\alpha = \lim_{r\rightarrow\infty}(2 Z-U_{p}(r))/r$.

Numerov again

Poisson equation does not have the first order derivative, hence it can also be more efficiently solved by the Numerov algorithm.

Recall the Numerov algorithm as appropriate for the Poisson equation

\begin{eqnarray} x(h)+x(-h) = 2x(0)+h^2 (f(0)x(0)+u(0))+\frac{2}{4!}h^4 x^{(4)}(0)+O(h^6) \end{eqnarray}\begin{equation} x^{(4)}\sim \frac{u_{i+1}-2 u_i+u_{i-1}}{h^2} \end{equation}

Inserting the fourth order derivative into the above recursive equation (forth equation in his chapter), we get

\begin{equation} x_{i+1}-2 x_i+x_{i-1}=h^2 u_i +\frac{h^2}{12}(u_{i+1}-2 u_i+u_{i-1}) \end{equation}

If we switch to a new variable $w_i=x_i-\frac{h^2}{12}u_i$ we are left with the following equation

\begin{equation} w_{i+1} -2 w_i + w_{i-1} = h^2 u_i+O(h^6) \end{equation}

The variable $x$ needs to be recomputed at each step with $x_i=(w_i+\frac{h^2}{12}u_i)$.

Step1 routine HartreeU

For generic density the following routine will work:

Step 2 : Compute the exchange correlation potential.

Note that $V_{xc}(r)=V_{xc}(\rho(r))$ is unquely determined by the electron charge density $\rho(r)$. If we know $\rho$, we can instantly compute $V_{xc}$ by the module provided parametrized function.

Download the module excor.py from http://www.physics.rutgers.edu/~haule/509/src_prog/python/homework5/ and import it in your code.

Instantiate the ExchangeCorrelation object by

exc = ExchangeCorrelation()

and used it, for example, by

exc.Vx(rs(rho[i]))+exc.Vc(rs(rho[i]))

where $r_s = ({4\pi\rho/3})^{-1/3}$.

Be careful: The energy unit in "excor.py" is Hartree and not Rydergs. Hence, you need to multiply energy or potential by 2.

Step 3: Find bound states using Hartree and XC potential.

Add the Hartree potential and the exchange-correlation potential to the Schroedinger equation and find bound states of the atom.

The Schroedinger equation is \begin{equation} u^{''}(r) = \left(\frac{l(l+1)}{r^2}-\frac{2 Z}{r} + V_{H}(r)+V_{XC}(r)-\varepsilon\right)u(r)=0. \end{equation} or \begin{equation} u^{''}(r) = \left(\frac{l(l+1)}{r^2}+\frac{U_H - 2 Z +r V_{XC}}{r}-\varepsilon\right)u(r)=0. \end{equation}

Step 4: Compute the new electron density

by filling the lowest $Z$ eigenstatates.

Step 5: Admix the new and the old density

(50\% of the old and 50\% of the new should work) and use the resulting density to compute the new Hartree and exchange-correlation potential.

Iterate Steps 1 to Step 5 until self-consistency is achieved.

For oxygen, the total energy in this implementation is : -74.47303426133809 Hartree, while NIST shows -74.473077 Hartree. This is in excellent agreement. The small loss of accuracy is mainly due to the mesh used for R.