# 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?

• (2) is answered here and (1) is essentially answered here. As for the rest, I worry that others might vote to close due "advice/suggestions" being "opinion-based" or needing better focus. I have not voted to close, but this is just what I worry others might do, based on my experience on SE :) – Nike Dattani May 24 at 5:35
• @NikeDattani. Thanks for pointing it out. I have rephrased the question. If it is not a worthy question, let the community close it :) – Thomas May 24 at 5:55
• I don't see any advantages to do that. – Camps May 24 at 18:20

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! – Nike Dattani May 24 at 14:52
• +10 for tinydft! – Charlie Crown May 25 at 16: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. – Paulie Bao May 25 at 15:01
• No worries - just wanted you to know – Geoff Hutchison May 25 at 15:30

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!