DFT for an atom

Introduction

I have a new project on GitHub¹. The project is using Density Functional Theory to do calculations for an atom. The project is actually not so new, I’ve put it on GitHub more than three months ago, but it had some issues I had to solve and also I did not have patience until now to write a new blog entry. I mentioned the post on Quantum Scattering that I will have some other project that takes advantage of spherical symmetry, so here it is. A very simple one, since I don’t have time and motivation right now for complex projects.

Program in action

As usual, I have a video showing the program in action. This is with an older version, the newer one is much better for heavy atoms:

This time I decided to not chart anything, it just shows some computed values. I had it first writing in the console, then I changed the code to display it in a rich text control. For that I simply redirected cout.

Theory

Some description

I won’t write much here, the references should be enough to figure the code out. For solving an atom, one might want to go relativistic (relativity does matter, especially for heavy atoms, the electrons are fast enough to need relativity). In such a case you either go with the Dirac equation or you go with relativistic corrections. In order to keep the project simple, I chose to ignore this. If you want to see what is needed, check out the Martin book or this link. Also this link should be very helpful: How to build an atom². I wanted a simple project, so I didn’t go that path. To benefit from the spherical symmetry, I also went with Local Density Approximation in 1D. For an atom, despite having spherical symmetry for the nucleus potential, there is no spherical symmetry in general, only for fully occupied shells, in the case of LDA. In general the symmetry is cylindrical, but with Local Spin Density Approximation you also have it for half-filled shells. It would be reasonably simple to implement LSDA, you just go with two densities instead of one, but I chose to go with LDA, for the purpose of this blog it should work well enough even for atoms that are not in the noble gases column. For solving the Kohn-Sham equation I chose the shooting method with the Numerov method. Since I already used that for Quantum Scattering, I thought that it would be too simple to also use it to solve the Poisson equation, so I went there with the full multigrid method. One additional complexity is that the code is using a non uniform grid, the benefit is quite good to be worth it although it adds a little complexity. One simply goes with something like equation 23 from here, writes the functions as functions of n instead of position, applies the chain rule for derivatives, may do some other change of variable to simplify things and computations go over a regular grid that’s simpler that the original one (delta is 1, the difference between n+1 and n). The method is exposed as a problem in Thijssen’s book. I think it’s worth mentioning that the code is using the shooting method from different points, the core wave-functions are more localized than the outer shell ones. The starting point is decided based on the value of the approximation used for far starting values.

Here is a succint description with the main points:

• A non uniform grid is used.
• The electron density is set to a constant one in the beginning, one could go with a better one, the Thomas-Fermi one, but it would spare only of a few iteration steps and maybe it will converge in some situations when this code wouldn’t, but for its purpose seems to be good enough.
• The shooting method with Numerov is used for computing the solution for Kohn-Sham.
• Locating the starting energy range is using the fact that the radial wavefunction has n-l-1 nodes (n counted from 1) for the ground state. An excited state with the same l has more nodes and this helps locate the interval where to search for solution for a particular n and l.
• In the located interval, the searching is done using the bisection method.
• Shooting is done from different distances, depending on the energy.
• Solution is finally computed by ‘shooting’ from both far away and the nucleus and matching the solution in between. I chose to use as the matching point the point after the first maximum. You could do better, matching not only the amplitude but also the derivatives in the ‘meeting’ point (from left and right). In that case you might want to use a different matching point (such as the classical turning point).
• For solving the Poisson equation, a full multigrid method is used.
• Density mixing is used to obtain the new electron density to be used in the next step.

On this blog

Related with the multigrid method, I have this blog post that exposes it at an entry level: Relaxation Method. In this case though it is full multigrid, not as simple as there. A project that uses self-consistency is the Hartree-Fock one. A blog post about DFT theory is here and a project that uses DFT, but with a plane-waves basis is the DFT for a Quantum Dot one. I already mentioned the Quantum Scattering project for the Numerov method that is also used there. If you would like to go with Runge-Kutta instead of Numerov, you might find the Electric Field Lines and Chaos posts relevant.

Books

I already mentioned the Computational Physics book by Jos Thijssen³. It’s worth mentioning again, chapter 5 is on the Density Functional Theory and the problem 5.1 is on switching from the uniform grid to the non-uniform one. Another great book is Electronic Structure, Basic Theory and Practical Methods by Richard M Martin⁴. Chapter 10 treats the electronic structure of atoms, but a lot of the book deals with relevant information here and some other projects described on this blog.

Articles and so on

I won’t put many links here, for more please visit the relevant posts on this blog. I’ll add only some relevant to this project.
One nice article for the theory is the already mentioned How to build an atom by M.S.S Brooks². You’ll also find there fortran code, if you like fortran more than C++. I also recommend looking over this lecture from Rutgers University: DFT lecture⁵ and the associated C++ code project for the lecture⁶. NIST has some theory on the site, with some references, so it’s worth checking not only for verifying the results: NIST⁷. Also this arxiv paper is worth checking, it will lead you to a more serious fortran project: dftatom: A robust and general Schrödinger and Dirac solver for atomic structure calculations by Ondrej Certik, John E. Pask d, Jiri Vackar⁸. You may find a lot of info on the internet about the full multigrid method, here is just one link: A multigrid tutorial by William L Briggs and others⁹.

Since I mentioned fortran above, twice, I think it’s worth mentioning this: Fortran is still a thing.

The code

The code is relatively simple once you get the theory, it’s well under 2000 lines of code, probably well under 1000 for the code that matters for the atom computations. The relevant classes are in the DFT namespace, the most important being the Numerov class for the algorithm used for the shooting method used to solve the Kohn-Sham equation, the PoissonSolver for solving the Poisson equation and the DFTAtom class, this one being the one that does the computations using the other ones. There are other helper classes, they are either similar as used in other projects on this blog (as for example VWNExchCor) or they are straightforward (as for example Integral). The program uses wxWidgets¹⁰ as many other projects described on this blog. To ease up understanding of the code, I also added code that is using an uniform grid, for comparison. The uniform grid computation is in DFTAtom::CalculateUniform, while the non uniform one is in DFTAtom::CalculateNonUniform.

Results

First, here is how the non uniform grid looks like, for the hydrogen wavefunction:

Delta was 0.01, radius 10 and multigrid levels 10, that is, 1025 grid points. You can easily see that the grid points get much denser closer to the nucleus.

Let’s see how it goes for a bad scenario, Argon:

Step: 30
Energy 1s: -3204.75628814 Num nodes: 0
Energy 2s: -546.577960661 Num nodes: 1
Energy 2p: -527.533025107 Num nodes: 0
Energy 3s: -133.369144873 Num nodes: 2
Energy 3p: -124.172862647 Num nodes: 1
Energy 3d: -106.945006737 Num nodes: 0
Energy 4s: -31.2308038208 Num nodes: 3
Energy 4p: -27.1089854743 Num nodes: 2
Energy 4d: -19.4499946904 Num nodes: 1
Energy 4f: -8.95331847594 Num nodes: 0
Energy 5s: -5.88968292298 Num nodes: 4
Energy 5p: -4.40870280587 Num nodes: 3
Energy 5d: -1.91132966098 Num nodes: 2
Energy 6s: -0.626570734867 Num nodes: 5
Energy 6p: -0.293180043718 Num nodes: 4
Etotal = -21861.3469029 Ekin = 21854.6726982 Ecoul = 8632.01604609 Eenuc = -51966.1203929 Exc = -381.915254274

Finished!

1s2 2s2 2p6 3s2 3p6 3d10 4s2 4p6 4d10 4f14 5s2 5p6 5d10 6s2 6p6


I used 17 for ‘multigrid levels’ (that means 131073 nodes) 0.0001 for delta, mixing 0.5 and the max radius 50.
The results are not perfect, I guess with some other parameters they might be improved somewhat. The energy levels usually get all decimals given by NIST right, but occasionally the last one is wrong.
The problem is for total energies, the total energy gets three decimals right, the partial ones get three or four decimals right.
Here are the NIST values for comparison: NIST Radon.

For a lighter noble gas I get better results, for Argon for example:

Step: 29
Energy 1s: -113.800134222 Num nodes: 0
Energy 2s: -10.7941723904 Num nodes: 1
Energy 2p: -8.44343924178 Num nodes: 0
Energy 3s: -0.88338408662 Num nodes: 2
Energy 3p: -0.382330129715 Num nodes: 1
Etotal = -525.946199815 Ekin = 524.969812556 Ecoul = 231.45812437 Eenuc = -1253.13198253 Exc = -29.2421542116

Finished!

1s2 2s2 2p6 3s2 3p6


I used 14 for ‘multigrid levels’ (that means 16385 nodes), 0.0005 for delta, mixing 0.5 and the max radius 25.
The energy levels results match all given decimals from NIST. The total energy gets five decimals right, the kinetic, coulomb and nuclear energies even six, the exchange correlation one seems to be the worse, with only four decimals (maybe there is room for improvement there?).

Here is the NIST data for comparison: NIST Argon.

Let’s see how good it can be for 1025 points only:

Step: 29
Energy 1s: -113.800105669 Num nodes: 0
Energy 2s: -10.7941413182 Num nodes: 1
Energy 2p: -8.44340935717 Num nodes: 0
Energy 3s: -0.883372166507 Num nodes: 2
Energy 3p: -0.382319656694 Num nodes: 1
Etotal = -525.945824275 Ekin = 524.965679469 Ecoul = 231.458119091 Eenuc = -1253.12751443 Exc = -29.2421084002

Finished!

1s2 2s2 2p6 3s2 3p6

Multigrid levels was set to 10, radius to 15, grid delta to 0.006.

Four decimals for energy levels and two for total energy and the partial ones, not bad at all for so small number of points. The non uniform grid helps a lot, you would need a lot more points if using an uniform one.

Conclusions

As always, please point out any issues/bugs or improvements/enhancements ideas. The program is far from perfect, I’m told you can get much better results, especially for heavy atoms, but for now the results seem of for the purpose of this blog. I might use it in the future in other project(s).

CodeProject