Goal of this lecture

DFT introduction #2


Team Rutgers

(Michele, Ezekiel & Valeria)

Rutgers University-Newark
An interactive session

Retrieve this presentation at:


https://github.com/asesma-org/miniASESMA2026

mini ASESMA 2026 \(\bullet\) Accra, Ghana \(\bullet\) June 17, 2026

Goal of this lecture#

  • Understand what basis sets are

  • Introduce basics about basis sets (PWs and local orbitals like GTOs)

  • Discretize the KS equations so they can be solved by a computer

  • Focus on PW (GTOs done on your own, numerical orbitals introduced by Javier later)

Split into N groups#

  • Assign number 1-N to each student

  • Groups sit together

  • Possibly have 1 instructor per group

from IPython.display import IFrame
import numpy as np
import matplotlib.pyplot as plt
from ase.build import molecule

My assumptions#

  • Basic linear algebra (vectors, matrices, operators, wavefunctions)

  • Undergraduate QM

  • Dirac notation (i.e., i-th vector \(\to |i\rangle\) or \(|\psi_i\rangle\), scalar product between i and j \(\to \langle \psi_i | \psi_j\rangle\))

  • Kohn-Sham equations

Why are basis sets needed?#

  • They provide a representation

  • Allow us to discretize the problem

Let’s expand a function in a basis#

\[ \psi_k(r) = \sum_i C_{ik} \chi_i(r) \]

  • \(\{\chi_i(r)\}\) are basis functions forming a basis set

  • \(\{C_{ik}\}\) are expansion coefficients, \(C_{ik}\in \mathbb{C}\)

Examples of basis sets#

  • Plane waves (PW)

\[ \chi_{G}(r) = \frac{1}{\sqrt{\Omega}} e^{iGr} \]

  • Gaussian-Type Orbitals (GTOs)

\[ \chi_{\mu,\alpha}(r) = \left(\frac{2\alpha}{\pi}\right)^{d/2} e^{-\alpha |r-\mu|^2} \]

  • Numerical Orbitals (Javier will introduce later this week)

  • …

What do they look like? Let’s use Python#

L=10; X = np.linspace(start=0,stop=L,num=1000,endpoint=False)
def pw(x,G):
    return 1.0/np.sqrt(L)*np.exp(1j*2*np.pi*G*x)
def gto(x,center,alpha):
    return np.sqrt(alpha/np.pi)*np.exp(-alpha*(x-center)**2)

fig = plt.figure()#(figsize=(8,5))
plt.plot(X,np.real(pw(X,G=1/L)), label=r"Re[PW] with $G=1/L$")
plt.plot(X,np.imag(pw(X,G=1/L)), label=r"Im[PW] with $G=1/L$")
plt.plot(X,gto(X,center=5,alpha=1),label=r"GTO centered in $x=5$, with $\alpha=1$")
plt.plot(5, 0, 'ko', markersize=10); plt.text(x=4.8,y=-0.1,s="ion position")
plt.ylim([-.5,1]); plt.legend()

<matplotlib.legend.Legend at 0x13a53cdf0>
_images/671f8a89799826ef8318a2fe6d995711062033da4d574616b49f149aaecb8526.png

Challenge 1#

  • Are GTOs an orthonormal basis set? Justify.

  • Are PWs an orthonormal basis set? Justify.

The condition for orthonormality \(S_{ij}=\langle \chi_i | \chi_j\rangle = \delta_{ij}\).

How many plane waves are needed?#

  • You will attack this question in your Homework

  • The homework will introduce the PW cutoff needed by most PW codes, like QE

  • For QEpy, you will see something like (all values in Ry):

qe_options = {}
qe_options["&electrons"] = {}
qe_options["&electrons"]["ecutwf"] = 80
qe_options["&electrons"]["ecutrho"] = 320

  • ecutwf is the kinetic energy of the PW with the largest \(G\) point in Ry

  • ecutrho is the kinetic energy of the PW with the largest \(G\) point used to Fourier transform the electron density \(n(r)\) in Ry

  • Let me know during the break if you wish to know why these two cutoffs are different.

Challenge 1.2#

  • Convert 1 Ha into Ry

  • Convert 1 Ry into eV

How do basis sets “discretize” the KS equations?#


Working with your group mates, discretize the KS equations
\[-\frac{1}{2}\nabla^2 \phi_k(r) + v_s(r) \phi_k(r) = \varepsilon_k \phi_k(r)\]
Use: \(\{\chi_i(r)\}\) such that \(\phi_k(r) = \sum_i C_{ik}\chi_i(r)\)
Consider: KS Hamiltonian matrix elements are \(H_{ji}=\langle \chi_j \mid \hat H \mid \chi_i \rangle\), and the overlap elements among the basis functions are \(S_{ji}=\langle \chi_j \mid \chi_i \rangle\).

Discretized KS equations#

The final result should be

\[ \sum_i H_{ji}C_{ik} = \varepsilon_k \sum_i S_{ji} C_{ik} \]

or in matrix form

\[ \color{green}{\mathbb{H}\mathbb{C} = \mathbb{S}\mathbb{C}\boldsymbol{\varepsilon}} \]

  • \(\mathbb{H}\): The KS Hamiltonian, mean field Hamiltonian or Fock matrix

  • \(\mathbb{C}\): The matrix of the coefficients. Each column represents a KS orbital. In PW, these are simply the bands

  • \(\varepsilon\): A diagonal matrix with the orbital energies, \(\varepsilon_k\). Also KS eigenvalues, …

  • \(\mathbb{S}\): Overlap matrix. Diagonal in PW.

Matrix elements of \(\hat T\) in plane waves, needed to solve \(\mathbb{HC=SC\varepsilon}\)#

\[ \langle \chi_i | \hat T | \chi_j \rangle = \left\langle \chi_i \left| -\frac{1}{2} \nabla^2 \right| \chi_j \right\rangle =-\frac{1}{2\Omega} \int dr ~ -G_j^2 ~ e^{i(G_j-G_i)}r = \mathbf{\color{red}{\frac{G_j^2}{2}\delta_{ij}}} \]
It's diagonal! ... and cheap to compute ...

Matrix elements of \(v_s(r)\), needed to solve \(\mathbb{HC=SC\varepsilon}\)#

\[ \langle \chi_i | v_s | \chi_j\rangle = \frac{1}{\Omega}\int dr v_s(r) e^{i(G_j-G_i)r} \]
Not a diagonal matrix

But… we can use the \(G\) components of \(\tilde v_s(G)\):

\[ \langle \chi_i | v_s | \chi_j\rangle = \frac{1}{\Omega}\tilde v_s(G_j-G_i) \]

which gives the matrix elements of \(v_s\) almost directly.

Hartree potential, the (computationally) most difficult piece of \(v_s(r)\)#

\[ v_H(r) = \int dr' \frac{n(r')}{|r-r'|} \equiv \color{red}{n \star \frac{1}{|r|}} \qquad \text{Poisson equation} \]
Let's use the convolution theorem : \(f = h \star g \;\longrightarrow\; \tilde f = \tilde h \cdot \tilde g\)
\[ \tilde v_H(G) = \tilde n(G) \frac{4\pi}{G^2} \to \langle \chi_i | v_H | \chi_j\rangle = \frac{1}{\Omega} \tilde v_H(G_j-G_i) \]
Thanks to Fast Fourier Transforms (FFT), the \(v_H\) potential is computed with \(\mathcal{O}(N\ln N)\) operations!

Challenge 2#

  • Can you derive the equations we just presented?

Now we can solve \(\mathbb{HC=C\varepsilon}\)?#

\[\begin{split} H_{ij} = \left[ \frac{G_i^2}{2} \delta_{ij} + \frac{1}{\Omega}\tilde v_s(G_j - G_i) \right] \\ \text{ ... and ... } S_{ij} = \delta_{ij} \end{split}\]

Remember \(v_s(r) = v_{eN}(r) + v_H(r) + v_{xc}(r)\).

Is it enough?
\(v_H\) and \(v_{xc}\) depend on \(n(r)\). Do we know what \(n(r)\) is?
We'll get to this later

First, a note on PWs#

  • Many plane waves are needed to represent “spiky” or highly oscillatory features, such as \(v_{eN}(r)\) and \(n(r)\) in the ionic core region.

  • Thus, pseudopotentials are typically employed to:

    1. avoid the Coulomb singularity of the electron–nuclear interaction potential;

    2. remove the need to represent the atomic cores, where densities and orbitals are highly oscillatory. More on this from Javier tomorrow!

Why are PW basis so useful?#

  • Orthonormal!

  • Systematically improvable: higher \(N\) (or equivalently \(\max\left[ G \right]\)) leads to better results.

  • Evaluation of the kinetic-energy operator and solution of the Poisson equation are easy, with \(\mathcal{O}(N\ln N)\) complexity.

Break time!
Contents