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.
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.
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.
I think this is an ambitious task. I suggest you to divide the writing DFT program into two subproblems:
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!
