Chapter 4: Density Functional Theory (DFT) | Introduction to Quantum Chemistry

🎯 Learning Objectives

Understand the Hohenberg-Kohn theorems and the uniqueness of the electron density
Learn the derivation and physical meaning of the Kohn-Sham equations
Understand the role of exchange-correlation functionals and their approximations
Learn the local density approximation (LDA) and the generalized gradient approximation (GGA)
Master the concepts of hybrid functionals (B3LYP, PBE0)
Understand plane-wave basis sets and pseudopotentials
Compute the band structure and density of states of solids
Learn applications to real materials calculations

📖 The Hohenberg-Kohn Theorems

Fundamental Principle of DFT

The ground state of an N-electron system is determined solely by the electron density \(\rho(\mathbf{r})\).

First theorem (uniqueness theorem):

The external potential \(v_{ext}(\mathbf{r})\) is uniquely determined by the electron density \(\rho(\mathbf{r})\) (up to a constant). Therefore, all ground-state properties can be expressed as functionals of \(\rho(\mathbf{r})\).

Second theorem (variational principle):

The energy functional \(E[\rho]\) takes its minimum value at the ground-state density:

\[ E[\rho] = T[\rho] + V_{ext}[\rho] + V_{ee}[\rho] \geq E_0 \]

Here \(E_0\) is the exact ground-state energy.

The Kohn-Sham Equations

The real system is mapped onto a non-interacting system that reproduces the same electron density:

\[ \left[-\frac{1}{2}\nabla^2 + v_{KS}(\mathbf{r})\right]\phi_i(\mathbf{r}) = \varepsilon_i \phi_i(\mathbf{r}) \]

Kohn-Sham potential:

\[ v_{KS}(\mathbf{r}) = v_{ext}(\mathbf{r}) + v_H(\mathbf{r}) + v_{xc}(\mathbf{r}) \]

\(v_{ext}\): external potential (nucleus-electron attraction)
\(v_H\): Hartree potential (classical electron-electron repulsion)
\(v_{xc}\): exchange-correlation potential (quantum effects)

Electron density:

\[ \rho(\mathbf{r}) = \sum_{i=1}^{N/2} 2|\phi_i(\mathbf{r})|^2 \]

💻 Example 4.1: Kohn-Sham DFT Calculation (1D)

Kohn-Sham Calculation for a 1D Harmonic Well

We solve the Kohn-Sham equations numerically for a simple 1D system:

External potential: \(v_{ext}(x) = \frac{1}{2}\omega^2 x^2\)

LDA approximation (1D): \(\varepsilon_{xc}[\rho] = C_x \rho^{4/3}\)

Python implementation: 1D Kohn-Sham DFT

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh
from scipy.integrate import simps

class KohnSham1D:
    """1D Kohn-Sham DFT calculation"""

    def __init__(self, N_electrons=2, omega=1.0, L=10.0, N_grid=200):
        """
        N_electrons: number of electrons
        omega: strength of the harmonic potential
        L: size of the computational domain
        N_grid: number of grid points
        """
        self.N_electrons = N_electrons
        self.omega = omega
        self.L = L
        self.N_grid = N_grid

        # Spatial grid
        self.x = np.linspace(-L/2, L/2, N_grid)
        self.dx = self.x[1] - self.x[0]

        # Kinetic energy operator (finite difference)
        self.T = self.kinetic_energy_matrix()

    def kinetic_energy_matrix(self):
        """Kinetic energy matrix (3-point finite difference)"""
        N = self.N_grid
        dx = self.dx

        # -1/2 d²/dx²
        T = np.zeros((N, N))
        for i in range(1, N-1):
            T[i, i-1] = -1 / (2 * dx**2)
            T[i, i] = 1 / dx**2
            T[i, i+1] = -1 / (2 * dx**2)

        # Boundary conditions (ψ=0)
        T[0, 0] = T[-1, -1] = 1e10  # Large value enforces ψ≈0

        return T

    def external_potential(self):
        """External potential (harmonic well)"""
        return 0.5 * self.omega**2 * self.x**2

    def hartree_potential(self, rho):
        """Hartree potential (1D)"""
        # Simplified Coulomb potential (1D)
        # v_H(x) = ∫ ρ(x') / |x - x'| dx'
        v_H = np.zeros_like(self.x)

        for i, xi in enumerate(self.x):
            # Avoid the singularity
            denominator = np.abs(self.x - xi) + 1e-10
            v_H[i] = simps(rho / denominator, self.x)

        return v_H

    def xc_potential_lda(self, rho):
        """Exchange-correlation potential (LDA, 1D)"""
        # 1D LDA: v_xc = d(ε_xc ρ)/dρ
        # ε_xc ∝ ρ^(4/3) → v_xc ∝ ρ^(1/3)
        C_x = -0.5  # Exchange energy constant (1D, simplified)

        v_xc = (4/3) * C_x * np.sign(rho) * np.abs(rho)**(1/3)
        v_xc[np.abs(rho) < 1e-10] = 0  # Numerical stability

        return v_xc

    def kohn_sham_potential(self, rho):
        """Kohn-Sham potential"""
        v_ext = self.external_potential()
        v_H = self.hartree_potential(rho)
        v_xc = self.xc_potential_lda(rho)

        return v_ext + v_H + v_xc

    def solve_kohn_sham(self, rho_init=None, max_iter=50, conv_threshold=1e-6):
        """Self-consistent solution of the Kohn-Sham equations"""
        # Initial density
        if rho_init is None:
            # Gaussian initial density
            rho = np.exp(-self.x**2) / np.sqrt(np.pi)
            rho = rho * self.N_electrons / simps(rho, self.x)

        energies = []
        converged = False

        for iteration in range(max_iter):
            # Kohn-Sham potential
            v_KS = self.kohn_sham_potential(rho)

            # Hamiltonian matrix
            V_diag = np.diag(v_KS)
            H = self.T + V_diag

            # Solve the eigenvalue problem
            eigenvalues, eigenvectors = eigh(H)

            # Occupied orbitals (lowest N/2)
            n_occupied = self.N_electrons // 2

            # New density
            rho_new = np.zeros_like(self.x)
            for i in range(n_occupied):
                rho_new += 2 * eigenvectors[:, i]**2  # Spin pair

            # Normalization
            rho_new = rho_new * self.N_electrons / simps(rho_new, self.x)

            # Energy calculation
            E_total = self.compute_energy(rho_new, eigenvalues[:n_occupied])
            energies.append(E_total)

            # Convergence check
            if iteration > 0:
                delta_rho = np.max(np.abs(rho_new - rho))
                if delta_rho < conv_threshold:
                    converged = True
                    break

            # Density mixing (stabilizes convergence)
            alpha_mix = 0.3
            rho = alpha_mix * rho_new + (1 - alpha_mix) * rho

        return {
            'converged': converged,
            'iterations': iteration + 1,
            'energy': E_total,
            'density': rho_new,
            'orbitals': eigenvectors[:, :n_occupied],
            'orbital_energies': eigenvalues[:n_occupied],
            'energy_history': energies
        }

    def compute_energy(self, rho, orbital_energies):
        """Total energy"""
        # Kohn-Sham energy
        E_KS = 2 * np.sum(orbital_energies)

        # Double-counting correction (simplified)
        v_H = self.hartree_potential(rho)
        E_H = 0.5 * simps(rho * v_H, self.x)

        v_xc = self.xc_potential_lda(rho)
        E_xc = simps(rho * v_xc, self.x) * (3/4)  # Simplified correction

        E_total = E_KS - E_H - E_xc * (1/4)

        return E_total

# Run the DFT calculation
dft = KohnSham1D(N_electrons=2, omega=1.0, L=10, N_grid=200)
result = dft.solve_kohn_sham()

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Electron density
ax1 = axes[0, 0]
ax1.plot(dft.x, result['density'], 'b-', linewidth=2, label='DFT density')
ax1.fill_between(dft.x, 0, result['density'], alpha=0.3)
ax1.set_xlabel('Position x')
ax1.set_ylabel('Electron density ρ(x)')
ax1.set_title('Electron Density Distribution')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Kohn-Sham orbitals
ax2 = axes[0, 1]
for i, phi in enumerate(result['orbitals'].T):
    # Check normalization
    norm = np.sqrt(simps(phi**2, dft.x))
    phi_normalized = phi / norm
    ax2.plot(dft.x, phi_normalized + i*0.5, linewidth=2, label=f'φ_{i+1}')

ax2.set_xlabel('Position x')
ax2.set_ylabel('Kohn-Sham orbitals φ_i(x)')
ax2.set_title('Kohn-Sham Orbitals')
ax2.legend()
ax2.grid(True, alpha=0.3)

# Potential components
ax3 = axes[1, 0]
v_ext = dft.external_potential()
v_H = dft.hartree_potential(result['density'])
v_xc = dft.xc_potential_lda(result['density'])
v_KS = v_ext + v_H + v_xc

ax3.plot(dft.x, v_ext, 'b-', linewidth=2, label='External')
ax3.plot(dft.x, v_H, 'r-', linewidth=2, label='Hartree')
ax3.plot(dft.x, v_xc, 'g-', linewidth=2, label='XC (LDA)')
ax3.plot(dft.x, v_KS, 'k--', linewidth=2, label='Total KS')

ax3.set_xlabel('Position x')
ax3.set_ylabel('Potential')
ax3.set_title('Kohn-Sham Potential Components')
ax3.legend()
ax3.grid(True, alpha=0.3)
ax3.set_ylim([0, 10])

# Convergence history
ax4 = axes[1, 1]
ax4.plot(result['energy_history'], 'go-', linewidth=2, markersize=6)
ax4.set_xlabel('SCF iteration')
ax4.set_ylabel('Total energy')
ax4.set_title('DFT-SCF Convergence')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('qchem_kohn_sham_dft.png', dpi=300, bbox_inches='tight')
plt.show()

# Numerical results
print("=== Kohn-Sham DFT calculation (1D) ===\n")
print(f"Number of electrons: {dft.N_electrons}")
print(f"Converged: {result['converged']} ({result['iterations']} iterations)")
print(f"Total energy: {result['energy']:.6f}")
print(f"\nKohn-Sham orbital energies:")
for i, eps in enumerate(result['orbital_energies']):
    print(f"  φ_{i+1}: ε = {eps:.6f}")

print(f"\nNormalization of the electron density:")
total_electrons = simps(result['density'], dft.x)
print(f"  ∫ρ(x)dx = {total_electrons:.6f} (target: {dft.N_electrons})")

💻 Example 4.2: Exchange-Correlation Functionals

Hierarchy of Exchange-Correlation Functionals

Local density approximation (LDA):

\[ E_{xc}^{LDA}[\rho] = \int \rho(\mathbf{r}) \varepsilon_{xc}(\rho(\mathbf{r})) d^3r \]

Uses the energy of the uniform electron gas
Exchange: \(\varepsilon_x(\rho) = -C_x \rho^{1/3}\), \(C_x = \frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3}\)
Correlation: parameterizations such as Vosko-Wilk-Nusair (VWN)

Generalized gradient approximation (GGA):

\[ E_{xc}^{GGA}[\rho] = \int \rho(\mathbf{r}) \varepsilon_{xc}(\rho, |\nabla\rho|) d^3r \]

Accounts for the density gradient, improving inhomogeneous systems
Representative examples: PBE, BLYP, PW91

Hybrid functionals:

\[ E_{xc}^{hybrid} = aE_x^{HF} + (1-a)E_x^{DFT} + E_c^{DFT} \]

Mix in a fraction of Hartree-Fock exchange
B3LYP: \(a=0.2\), PBE0: \(a=0.25\)
High accuracy for band gaps and reaction energies

Python implementation: comparison of exchange-correlation functionals

import numpy as np
import matplotlib.pyplot as plt

def exchange_lda(rho):
    """LDA exchange energy density"""
    C_x = (3/4) * (3/np.pi)**(1/3)
    return -C_x * rho**(4/3)

def exchange_pbe(rho, grad_rho, kappa=0.804, mu=0.2195):
    """PBE exchange energy density (simplified)"""
    # Reduced density gradient
    s = np.abs(grad_rho) / (2 * (3*np.pi**2)**(1/3) * rho**(4/3))

    # Enhancement factor
    F_x = 1 + kappa - kappa / (1 + mu * s**2 / kappa)

    # LDA × enhancement
    epsilon_x_lda = exchange_lda(rho) / rho
    return rho * epsilon_x_lda * F_x

def correlation_vwn(rho):
    """VWN correlation energy density (simplified)"""
    # Vosko-Wilk-Nusair parameters (paramagnetic)
    A = 0.0310907
    x0 = -0.10498
    b = 3.72744
    c = 12.9352

    r_s = (3 / (4 * np.pi * rho))**(1/3)
    x = np.sqrt(r_s)

    X = x**2 + b*x + c
    X0 = x0**2 + b*x0 + c
    Q = np.sqrt(4*c - b**2)

    epsilon_c = A * (
        np.log(x**2 / X) +
        2*b/Q * np.arctan(Q / (2*x + b)) -
        b*x0/X0 * (
            np.log((x - x0)**2 / X) +
            2*(b + 2*x0)/Q * np.arctan(Q / (2*x + b))
        )
    )

    return rho * epsilon_c

# Electron density range
rho_range = np.logspace(-2, 1, 100)

# Gradient strength (GGA)
grad_rho_weak = 0.1 * rho_range**(4/3)
grad_rho_strong = 1.0 * rho_range**(4/3)

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# LDA exchange energy density
ax1 = axes[0, 0]
epsilon_x_lda = exchange_lda(rho_range) / rho_range

ax1.plot(rho_range, epsilon_x_lda, 'b-', linewidth=2, label='LDA exchange')
ax1.set_xlabel('Electron density ρ')
ax1.set_ylabel('ε_x (energy per electron)')
ax1.set_title('LDA Exchange Energy Density')
ax1.set_xscale('log')
ax1.legend()
ax1.grid(True, alpha=0.3)

# PBE vs LDA (weak gradient)
ax2 = axes[0, 1]
epsilon_x_pbe_weak = exchange_pbe(rho_range, grad_rho_weak) / rho_range

ax2.plot(rho_range, epsilon_x_lda, 'b-', linewidth=2, label='LDA')
ax2.plot(rho_range, epsilon_x_pbe_weak, 'r-', linewidth=2, label='PBE (weak ∇ρ)')
ax2.set_xlabel('Electron density ρ')
ax2.set_ylabel('ε_x')
ax2.set_title('GGA Effect (Weak Gradient)')
ax2.set_xscale('log')
ax2.legend()
ax2.grid(True, alpha=0.3)

# PBE vs LDA (strong gradient)
ax3 = axes[1, 0]
epsilon_x_pbe_strong = exchange_pbe(rho_range, grad_rho_strong) / rho_range

ax3.plot(rho_range, epsilon_x_lda, 'b-', linewidth=2, label='LDA')
ax3.plot(rho_range, epsilon_x_pbe_strong, 'r-', linewidth=2, label='PBE (strong ∇ρ)')
ax3.set_xlabel('Electron density ρ')
ax3.set_ylabel('ε_x')
ax3.set_title('GGA Effect (Strong Gradient)')
ax3.set_xscale('log')
ax3.legend()
ax3.grid(True, alpha=0.3)

# Correlation energy (VWN)
ax4 = axes[1, 1]
epsilon_c_vwn = correlation_vwn(rho_range) / rho_range

ax4.plot(rho_range, epsilon_c_vwn, 'g-', linewidth=2, label='VWN correlation')
ax4.set_xlabel('Electron density ρ')
ax4.set_ylabel('ε_c')
ax4.set_title('LDA Correlation Energy (VWN)')
ax4.set_xscale('log')
ax4.legend()
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('qchem_xc_functionals.png', dpi=300, bbox_inches='tight')
plt.show()

# Numerical comparison
print("\n=== Comparison of exchange-correlation functionals ===\n")
print("Exchange energy at a representative density (ρ = 0.1):")
rho_test = 0.1
print(f"  LDA exchange: {exchange_lda(rho_test) / rho_test:.6f}")

grad_test_weak = 0.01
grad_test_strong = 0.1
print(f"  PBE (|∇ρ| = {grad_test_weak}): {exchange_pbe(rho_test, grad_test_weak) / rho_test:.6f}")
print(f"  PBE (|∇ρ| = {grad_test_strong}): {exchange_pbe(rho_test, grad_test_strong) / rho_test:.6f}")

print(f"\nCorrelation energy:")
print(f"  VWN correlation: {correlation_vwn(rho_test) / rho_test:.6f}")

print("\nGuidelines for functional selection:")
print("  LDA: metals, highly symmetric systems")
print("  GGA (PBE): molecules, surfaces, general solids")
print("  Hybrid (B3LYP, PBE0): molecular chemistry, band gaps")
print("  Meta-GGA (TPSS, SCAN): high-accuracy solid-state calculations")

💻 Example 4.3: Plane-Wave Basis and Pseudopotentials

DFT Calculations for Periodic Systems

Bloch's theorem:

Wavefunctions in a periodic potential:

\[ \psi_{n\mathbf{k}}(\mathbf{r}) = e^{i\mathbf{k}\cdot\mathbf{r}} u_{n\mathbf{k}}(\mathbf{r}) \]

Here \(u_{n\mathbf{k}}\) has the periodicity of the lattice.

Plane-wave expansion:

\[ \psi_{n\mathbf{k}}(\mathbf{r}) = \sum_{\mathbf{G}} c_{n\mathbf{k}}(\mathbf{G}) e^{i(\mathbf{k}+\mathbf{G})\cdot\mathbf{r}} \]

\(\mathbf{G}\) is a reciprocal lattice vector.

Pseudopotential approximation:

Replace core electrons with a pseudopotential
Treat only the valence electrons explicitly
Norm-conserving, ultrasoft, and PAW methods

Python implementation: band structure of a 1D periodic system

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh

class PlaneWave1D:
    """Plane-wave DFT for a 1D periodic system"""

    def __init__(self, a=1.0, V0=1.0, N_pw=11):
        """
        a: lattice constant
        V0: potential strength
        N_pw: number of plane waves (from -N_pw/2 to N_pw/2)
        """
        self.a = a
        self.V0 = V0
        self.N_pw = N_pw

        # Reciprocal lattice vectors
        self.G = 2 * np.pi / a * np.arange(-N_pw//2, N_pw//2 + 1)

    def periodic_potential(self, x):
        """Periodic potential V(x) = V0 cos(2πx/a)"""
        return self.V0 * np.cos(2 * np.pi * x / self.a)

    def hamiltonian_matrix(self, k):
        """Hamiltonian matrix at a k-point"""
        N = len(self.G)
        H = np.zeros((N, N), dtype=complex)

        for i, Gi in enumerate(self.G):
            for j, Gj in enumerate(self.G):
                if i == j:
                    # Kinetic energy
                    H[i, j] = 0.5 * (k + Gi)**2
                else:
                    # Fourier component of the potential
                    # V(x) = V0 cos(2πx/a) → V_G = V0/2 δ_{G,±G0}
                    G_diff = Gi - Gj
                    G0 = 2 * np.pi / self.a

                    if np.abs(G_diff - G0) < 1e-10:
                        H[i, j] = self.V0 / 2
                    elif np.abs(G_diff + G0) < 1e-10:
                        H[i, j] = self.V0 / 2

        return H

    def compute_bands(self, k_points):
        """Band structure calculation"""
        bands = []

        for k in k_points:
            H = self.hamiltonian_matrix(k)
            eigenvalues = eigh(H, eigvals_only=True)
            bands.append(eigenvalues)

        return np.array(bands)

# Band structure calculation
pw = PlaneWave1D(a=1.0, V0=2.0, N_pw=11)

# Brillouin zone: from -π/a to π/a
k_points = np.linspace(-np.pi/pw.a, np.pi/pw.a, 100)
bands = pw.compute_bands(k_points)

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Band structure
ax1 = axes[0, 0]
n_bands = min(5, bands.shape[1])  # First 5 bands

for i in range(n_bands):
    ax1.plot(k_points * pw.a / np.pi, bands[:, i], linewidth=2)

ax1.set_xlabel('k (units of π/a)')
ax1.set_ylabel('Energy')
ax1.set_title(f'Band Structure (V₀ = {pw.V0})')
ax1.axvline(0, color='k', linestyle='--', linewidth=1)
ax1.grid(True, alpha=0.3)

# V0 dependence of the band gap
ax2 = axes[0, 1]
V0_range = np.linspace(0, 5, 20)
band_gaps = []

k_gamma = 0  # Γ point
k_edge = np.pi / pw.a  # Brillouin zone boundary

for V0 in V0_range:
    pw_temp = PlaneWave1D(a=1.0, V0=V0, N_pw=11)

    # Bands at the Γ point
    H_gamma = pw_temp.hamiltonian_matrix(k_gamma)
    E_gamma = eigh(H_gamma, eigvals_only=True)

    # Bands at the zone boundary
    H_edge = pw_temp.hamiltonian_matrix(k_edge)
    E_edge = eigh(H_edge, eigvals_only=True)

    # Band gap (lowest excitation energy)
    gap = E_gamma[1] - E_gamma[0] if len(E_gamma) > 1 else 0
    band_gaps.append(gap)

ax2.plot(V0_range, band_gaps, 'ro-', linewidth=2, markersize=6)
ax2.set_xlabel('Potential strength V₀')
ax2.set_ylabel('Band gap')
ax2.set_title('V₀ Dependence of the Band Gap')
ax2.grid(True, alpha=0.3)

# Density of states (DOS)
ax3 = axes[1, 0]
# Dense k-point sampling
k_dense = np.linspace(-np.pi/pw.a, np.pi/pw.a, 500)
bands_dense = pw.compute_bands(k_dense)

# DOS calculation by the histogram method
E_min, E_max = bands_dense.min(), bands_dense.max()
E_bins = np.linspace(E_min, E_max, 100)
dos, _ = np.histogram(bands_dense.flatten(), bins=E_bins)

ax3.plot(dos, E_bins[:-1], 'b-', linewidth=2)
ax3.set_xlabel('Density of States')
ax3.set_ylabel('Energy')
ax3.set_title('Density of States (DOS)')
ax3.grid(True, alpha=0.3)

# Wavefunction (Γ point, lowest band)
ax4 = axes[1, 1]
k_gamma = 0
H_gamma = pw.hamiltonian_matrix(k_gamma)
eigenvalues, eigenvectors = eigh(H_gamma)

# Real-space reconstruction
x = np.linspace(0, pw.a, 200)
psi_0 = np.zeros_like(x, dtype=complex)

for i, G in enumerate(pw.G):
    psi_0 += eigenvectors[i, 0] * np.exp(1j * (k_gamma + G) * x)

ax4.plot(x / pw.a, np.real(psi_0), 'b-', linewidth=2, label='Re(ψ)')
ax4.plot(x / pw.a, np.abs(psi_0)**2, 'r-', linewidth=2, label='|ψ|²')

# Potential (normalized for display)
V_plot = pw.periodic_potential(x)
V_normalized = V_plot / np.max(np.abs(V_plot)) * np.max(np.abs(psi_0))
ax4.plot(x / pw.a, V_normalized, 'g--', linewidth=1, label='V(x) (scaled)')

ax4.set_xlabel('Position (x/a)')
ax4.set_ylabel('Wavefunction')
ax4.set_title('Γ-Point Wavefunction (Lowest Band)')
ax4.legend()
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('qchem_plane_wave_bands.png', dpi=300, bbox_inches='tight')
plt.show()

# Numerical results
print("\n=== Plane-wave basis and band structure ===\n")
print(f"Lattice constant: a = {pw.a}")
print(f"Potential strength: V₀ = {pw.V0}")
print(f"Number of plane waves: {pw.N_pw}")
print(f"\nEnergies at the Γ point (lowest 3 bands):")
H_gamma = pw.hamiltonian_matrix(0)
E_gamma = eigh(H_gamma, eigvals_only=True)
for i in range(min(3, len(E_gamma))):
    print(f"  Band {i+1}: E = {E_gamma[i]:.6f}")

print(f"\nBand gap: {E_gamma[1] - E_gamma[0]:.6f}")

📚 Summary

The Hohenberg-Kohn theorems establish that the ground state is determined solely by the electron density
The Kohn-Sham equations enable practical calculations through a mapping onto a non-interacting system
Exchange-correlation functionals approximately incorporate quantum many-body effects
LDA is accurate for uniform systems, while GGA improves inhomogeneous systems
Hybrid functionals mix in HF exchange and achieve high accuracy in molecular chemistry
The plane-wave basis is the natural basis for periodic systems, and fast FFT is available
Pseudopotentials allow core electrons to be handled efficiently
The band structure and density of states provide understanding of the electronic states of solids
DFT is the standard method for first-principles calculations in materials science
Appropriate functional selection determines the accuracy of the calculation

💡 Exercises

DFT calculation of the He atom: Compute the He atom (2 electrons) with Kohn-Sham DFT and compare the energies obtained with LDA and GGA.
Jacob's Ladder: Investigate the hierarchy of LDA, GGA, meta-GGA, hybrid, and double-hybrid functionals, and summarize the applicable range of each.
Band gap problem: Investigate why DFT (LDA/GGA) underestimates the band gap of semiconductors.
1D metal: Perform a plane-wave calculation with V₀=0 and confirm that the DOS of a metal is finite at the Fermi level.
k-point sampling: Investigate the k-point convergence of the Brillouin zone integration and implement the Monkhorst-Pack method.

🔬 Chapter 4: Density Functional Theory (DFT)