DFT simulation of the Helium atom

This is a density functional theory program for the calculation of the ground state energy of the helium atom, using the Slater X-alpha exchange potential and the parametrization of the Ceperley-Alder correlation potential (Parametrisation given by Perdew and Zunger). In the DFT, the electronic orbitals are solution to the Schrödinger equation (SE) which depends on the electron density rather than on the individual electron orbitals. As such the DFT orbitals have no individual meaning but are used to construct the charge density. The main result of DFT is that there exists a form of the exchange correlation potential, depending only on the electron density n(r), that yields the exact ground state energy and density. The form of this functional is however unknown and we should relay on approximation, such as the local density approximation (LDA, in this example we will consider the CA-LDA parameterization).

Walter Kohn and Lu Jeu Scham have proposed to write the total functional E[n] as E[n] = T_s[n] + E_H[n] + E_XC[n],

where the terms are respectively the non interacting system kinetic energy functional, the Hartree functional and the exchange functional.From δE[n]=0 we obtain the total potential

,

where V_N is the culombic potential from the nucleus, the second term is the Hartree potential and the last term is the exchange potential. If E_XC would be well known, then the self-consistent solution of the equation would give the exact ground-state density. Using the local density approximation (LDA) the potential is given by the sum of a exchange term and a correlation term: V_XC^LDA = V_X + V_C. The Slater X potential was used for the Exchange potential while for the correlation interaction the Quantum Monte Carlo Parametrization from Ceperley-Adler was used.

Structure of the program

Considering that the two electrons occupy the 1s-orbital, both the density and the Hartree potential are radially symmetric; thus, we can exploit this radial symmetry and solve the radial Schrödinger equation. The SE is solved directly using Verlet's algorithm. The program is organized in three steps:

1. Solution of the hydrogen-like radial SE using a simple integration algorithm combined with an interpolation routine in order to find the stationary states.
2. Calculation of the Hartree potential from the radial electron density.
3. Incorporation of the exchange-correlation potential (within the local density approximation, LDA).

Solution of the radial SE

The radial symmetry is exploited in order to solve the radial Schrödinger equation on a radial vector r using the Verlet algorithm in order to integrate inward and find the electron energy E

using the boundary conditions . The bisection method is used until the given precision on E is reached, checking the number of nodes in order to know if the correct energy is higher or lower with respect to the current one.

Computing the potentials

The nuclear potential is V_N = n_e/r. If we define U''_H(r) = -u(r)/r². We find U_H(r) using the Verlet algorithm with the boundary conditions .

Within the LDA V_X + V_C depend locally on n(r). The Slater potential was used fo the exchange interaction:

The correlation term is implemented using the Ceperley Alder parametrization:

with r_S = (4πn/3)^1/3.

The potential enegy for each potential is given by E_i = ∫dr V_i(r)u²(r) and the total energy of the system is therefore

E_tot = 2E - E_H - E_X/2 - E_C/2.

Structure of the project

The main files of the program are:

1. configuration.txt : This file contains all the major parameters used in the simulation such as the precision of the SE eignvalues (PREC_HSE) and the ground energy (PREC_DFT). The suggested higher bound for PREC_HSE is 10^-5 in order for the algorithm to work correctly. Also the range (R_MAX) and the number of points (SAMPLES) in the radial space are specified, as well as the energy range (HSE_E_MIN) in which to look for the SE eignvalues. In addition to these parameters, the paths where the output files are saved are specified. If the user wants to modify these parameters, he has to change them in the text file.
2. He.py : This file contains the main class of the programm. The He class has as attributes the radial WF (u) and density (rho), as well as the kinetic nergy of electrons (E_K), total energy (total_energy) the potentials (V_N,V_X,V_C,V_H,V). There are methods in order to compute the potentials as well as solving the SE and normalzing the WF. The main method is the hdft method which solves iteratively the SE updating the potentials as a new density is computed. Some methods in order to plot an save the results to files are present. The parameters used are taken from the configuration file.
3. DFT.py : This is the script that actually runs the simulation. After instancing an object of the class He, the methods to perform the simulation and to save the results are called.
4. testing.py : contains all the tests for the methods of He.py using hypothesis testing.
5. plotting.py : Contains the methods in order to plot the density and potentials. There are also methods to save the plots to .pdf files and the data about the attributes to a .txt file.
6. requirements.txt : Containes the dependencies needed in order for the script to run. The file was obtained using pipreqs. The oldest version of python able to run the script is python 3.7.

To run the script first install the dependencies:

pip install -r requirements.txt

And then launch the main script, with the possibility to specify the configuration file as sys.argv[1], for example as:

python DFT.py configuration.txt

The default configuration file if not specified through sys.argv[1] is configuration.txt .

Main results

The ground energy of the He atom using the Ceperley-Alder LDA approximation and the Slater exchange potential is:

E = -2.817 a.u