How to write my own density functional theory (DFT) code in Python?

I am a DFT user and at some point in the future, I would like to write my own DFT code in Python to help gain a deeper understanding of DFT. As mentioned in a previous answer people have written their own DFT codes to understand more deeply how the theory and algorithms work.

Where do I start? What are the important aspects I should consider?

Understanding, deriving, writing, testing and debugging an ab initio code can be a lenghty and tedious task. I'd like to provide a starting point for you here. If you just to it for pedagocial reasons, it might be advisable to start with the atomic problem and try to solve it with DFT. The effort for that is not too big, but it covers nearly all the nesseary physical principles to understand how one arrives at total energies and electron densities. With such a program, you can predict for example, energy levels of a certain element. Also, If you want to extend your program to treat crystalline material i. e., the written routines for the atomic problem will be of fundamental use. Core states are essential calculated like this.

The following websites provide some beginner-friendly introduction to the topic of simple DFT calculations of atoms.

I want to list down some the steps that are nessary for a simple solver written in python.

Setting up the equations

The starting point is the Kohn–Sham equation, which is the one-electron Schrödinger-like equation of a fictitious system of non-interacting particles. If you take the Born-Oppenheimer approximation into account and assume that the Kohn-Sham potential is radial symmetric, the task reduces to solve the radial part of the Kohn-Sham equation under some restrictions. Basically, you need to construct a solver that calculates the radial wave function and the eigenenergy for a given Kohn-Sham potential. For this you will need a shooting method similar to [1]. Basically, one solves the an inital value problem starting from $$r=0$$ and $$r=\infty$$ and change the energy till the solution matches at some point.

Calculating the Kohn-Sham potential

The effective Kohn-Sham potential is the sum of the core contributions ($$1/r$$ coloumb potential), the exchange-correlation (input here is the electron density) and the Hartree potential. In order to get the Hartree potential, you will need to solve the radial Poisson equation.

In order to set up a self-consistent cycle, you will need to calculate the poisson equation, which can be treaten as an inital value problem or also as an boundary value problem.

Remark to computational efficency in Python

The computational effort is considerable for bigger atoms. You might want to take a look some acceleration methods for Python like Numba.

http://numba.pydata.org/

• +10. Great first answer, and welcome! Hopefully we see more of you here! May 24, 2020 at 14:52
• +10 for tinydft! May 25, 2020 at 16:53
• Shooting method is just one possibility; there's a number of methods that can be used to solve the radial equations that could be used, instead; see onlinelibrary.wiley.com/doi/full/10.1002/qua.25968 for a review. Aug 28, 2020 at 7:53

For molecular DFT, there are definitely several tutorials around for major pieces.

The Psi4 code has a set of Jupyter notebooks called Psi4Numpy, working through each element of a quantum chemical program.

This includes:

There are a variety of other tutorials for response properties, geometry optimization, etc. Certainly you can work through much using the tutorials, particularly if you build off of Psi4 for "missing pieces."

Moreover, most (almost all?) functionals have been implemented as part of libxc which also has an example of calculating the XC energy for varying density.

• Sorry, I made a mistake. I replaced the answer. May 25, 2020 at 15:01
• No worries - just wanted you to know May 25, 2020 at 15:30

Computational Physics by Thijssen

The book starts at the simple Hartree problem, builds up to the Hartree-Fock problem, and eventually walk's the reader through constructing a pseudopotential planewave DFT program for silicon (gamma point only, but the generalization to other K points is straighforward). It covers the theory, provides coding exercises, and sample programs are on the author's website. I have found this book indispensable for actually understanding what I am doing when I set parameters in a calculation, interpreting outputs, etc.

• +1. Recommendation questions are fine, we even have a software-recommendations tag, and therefore your answer is also perfectly good! Aug 27, 2020 at 16:29
• @NikeDattani and Ty I removed the meta-discussion and added a bit of formatting. Hopefully you don't mind and feel free to edit more if you wish. Welcome to the site. Aug 27, 2020 at 19:00
• @CodyAldaz I don't mind at all, thanks for helping me out. I'll get the hang of how this forum works sooner-or-later! Aug 27, 2020 at 19:27

I think this is an ambitious task. I suggest you to divide the writing DFT program into two subproblems:

1. Writing an integration library for 1e and 2e electronic integrals.
2. Construct a SCF procedure.

The first part of the task is hard you might start from calling existing integration library and write an SCF code on your own. And also you could consider two aspect: how to write an correct code and how to implement it efficiently. Good luck!