7  Dichtefunktionaltheorie

Hartree-Fock sucht die optimale Wellenfunktion – eine antisymmetrische Funktion von \(3N_e\) Koordinaten. Für H₂O sind das 30 Koordinaten. Für ein Protein mit 10’000 Elektronen wären es 30’000. Das skaliert schlecht.

Die Dichtefunktionaltheorie (DFT) stellt die Frage radikal anders: Was wenn es genügt, die Elektronendichte \(\rho(\mathbf{r})\) zu kennen – eine Funktion von nur 3 Koordinaten?

---

7.1 Die Elektronendichte

Die Elektronendichte \(\rho(\mathbf{r})\) ist die Wahrscheinlichkeitsdichte, irgendein Elektron am Punkt \(\mathbf{r}\) zu finden:

\[\rho(\mathbf{r}) = N_e \int |\psi(\mathbf{r}, \mathbf{r}_2, \ldots, \mathbf{r}_{N_e})|^2\, d\mathbf{r}_2 \cdots d\mathbf{r}_{N_e} \tag{5.1}\]

Sie ist normiert auf die Elektronenzahl: \(\int \rho(\mathbf{r})\,d\mathbf{r} = N_e\).

In einer HF-Rechnung (oder KS-DFT) mit besetzten Orbitalen \(\psi_m\) ist sie:

\[\rho(\mathbf{r}) = 2\sum_{m=1}^{N_\text{occ}} |\psi_m(\mathbf{r})|^2 \tag{5.2}\]

Die Dichte ist eine Observable – sie lässt sich experimentell messen (Röntgenbeugung) und hängt nur von 3 Koordinaten ab, unabhängig von \(N_e\).

HinweisFunktional vs. Funktion

Ein Funktion \(f: \mathbb{R}^n \to \mathbb{R}\) weist einem Vektor eine Zahl zu. Ein Funktional \(F[\rho]\) weist einer Funktion eine Zahl zu.

Die Energie als Funktion der Dichte, \(E[\rho]\), ist ein Funktional: Für jede Dichte \(\rho(\mathbf{r})\) (eine Funktion von \(\mathbf{r}\)) liefert \(E[\rho]\) genau eine Zahl (die Energie).

---

7.2 Das Hohenberg-Kohn-Theorem (1964)

7.2.1 Existenztheorem (HK1)

WichtigHohenberg-Kohn-Existenztheorem

Die Grundzustandsenergie und alle anderen Grundzustandseigenschaften eines elektronischen Systems sind eindeutig durch die Grundzustandselektronendichte \(\rho_0(\mathbf{r})\) bestimmt.

Beweisidee (reductio ad absurdum): Angenommen, es gäbe zwei verschiedene externe Potentiale \(v(\mathbf{r})\) und \(v'(\mathbf{r})\) (d.h. zwei verschiedene Moleküle), die zur selben Grundzustandsdichte \(\rho_0\) führen. Verwendet man die Wellenfunktion des einen Systems als Testfunktion für das andere (Variationsprinzip), erhält man:

\[E_0 < E'_0 + \int \rho_0(\mathbf{r})\,[v(\mathbf{r}) - v'(\mathbf{r})]\,d\mathbf{r}\]

und symmetrisch:

\[E'_0 < E_0 - \int \rho_0(\mathbf{r})\,[v(\mathbf{r}) - v'(\mathbf{r})]\,d\mathbf{r}\]

Addition beider Ungleichungen ergibt \(E_0 + E'_0 < E'_0 + E_0\) – ein Widerspruch. Also war die Annahme falsch: Die Dichte bestimmt das Potential eindeutig, und damit das gesamte Molekül und alle seine Eigenschaften.

Was das bedeutet: Die Grundzustandsdichte enthält die vollständige Information über das System. Man braucht die hochdimensionale Wellenfunktion \(\psi(\mathbf{r}_1, \ldots, \mathbf{r}_{N_e})\) prinzipiell nicht.

7.2.2 Variationstheorem (HK2)

Das zweite HK-Theorem überträgt das Variationsprinzip auf Dichten:

WichtigHohenberg-Kohn-Variationstheorem

Für jede Probedichte \(\rho'(\mathbf{r}) \geq 0\) mit \(\int \rho' d\mathbf{r} = N_e\) gilt: \[E[\rho'] \geq E_0\] Gleichheit gilt genau für die wahre Grundzustandsdichte \(\rho_0\).

Das Minimierungsproblem (mit Lagrange-Multiplikator \(\mu\) für die Normierungsbedingung, vgl. Kapitel 00) lautet:

\[\delta\left[E[\rho] - \mu\!\int\rho(\mathbf{r})\,d\mathbf{r}\right] = 0 \quad \Longrightarrow \quad \mu = \frac{\delta E[\rho]}{\delta \rho(\mathbf{r})} \tag{5.3}\]

Der Lagrange-Multiplikator \(\mu\) ist das chemische Potential der Elektronen.

Das Problem: HK1 und HK2 garantieren die Existenz des Funktionals \(E[\rho]\), sagen aber nichts darüber, wie es aussieht. Wie lautet \(E[\rho]\) explizit?

fig, ax = plt.subplots(figsize=(9, 3.5))
ax.axis('off')

boxes = [
    (0.12, 0.55, 'Externes Potential\n$v(\\mathbf{r})$\n(Kernladungen, Positionen)',
     '#4393c3', 'white', 0.18, 0.28),
    (0.45, 0.55, 'Wellenfunktion\n$\\psi(\\mathbf{r}_1,\\ldots,\\mathbf{r}_{N_e})$\n$3N_e$ Koordinaten',
     '#d6604d', 'white', 0.16, 0.28),
    (0.78, 0.55, 'Elektronendichte\n$\\rho_0(\\mathbf{r})$\n3 Koordinaten',
     '#4dac26', 'white', 0.16, 0.28),
]

for x, y, txt, col, tc, w, h in boxes:
    ax.text(x, y, txt, ha='center', va='center', fontsize=10,
            color=tc, fontweight='bold',
            bbox=dict(boxstyle='round,pad=0.5', facecolor=col,
                      edgecolor='gray', lw=1.2, alpha=0.9),
            transform=ax.transAxes)

# Pfeile
for x1, x2, lbl, col in [
    (0.22, 0.36, 'Schrödinger-Gl.\n(eindeutig)', '#555'),
    (0.54, 0.69, 'HK1: eindeutig\n(Beweis durch Widerspruch)', '#2166ac'),
    (0.69, 0.54, 'HK2: Variation\n$E[\\rho\']\geq E_0$', '#2166ac'),
]:
    ax.annotate('', xy=(x2, 0.55), xytext=(x1, 0.55),
                xycoords='axes fraction', textcoords='axes fraction',
                arrowprops=dict(arrowstyle='->', color=col, lw=2.0))
    mx = (x1+x2)/2
    dy = 0.15 if '→' not in lbl else 0.15
    ax.text(mx, 0.55 + (0.22 if 'HK1' in lbl else -0.22),
            lbl, ha='center', va='center', fontsize=8.5,
            color=col, transform=ax.transAxes)

ax.text(0.45, 0.08,
        'HF-Welt: minimiere $E$ über Wellenfunktionen  |  DFT-Welt: minimiere $E[\\rho]$ über Dichten',
        ha='center', va='center', fontsize=10, color='#444',
        style='italic', transform=ax.transAxes)

ax.set_xlim(0, 1); ax.set_ylim(0, 1)
plt.tight_layout()
plt.show()

Abbildung 7.1: Das Hohenberg-Kohn-Theorem schlägt eine Brücke: Statt der 3Nₑ-dimensionalen Wellenfunktion genügt prinzipiell die 3-dimensionale Elektronendichte.

---

7.3 Die Kohn-Sham-Gleichungen

7.3.1 Das Referenzsystem

HK garantiert, dass \(E[\rho]\) existiert, aber nicht wie es lautet. Das eigentliche Problem: Der kinetische Energieterm \(T[\rho]\) ist unbekannt. Die Thomas-Fermi-Näherung dafür ist zu ungenau für Chemie.

W. Kohn und L.J. Sham (1965) fanden einen eleganten Ausweg:

WichtigKohn-Sham-Idee

Führe ein fiktives Referenzsystem aus \(N_e\) nicht-wechselwirkenden Elektronen ein, das dieselbe Grundzustandsdichte \(\rho_0(\mathbf{r})\) wie das reale System hat. Für dieses System ist der kinetische Energieterm exakt berechenbar.

Das Referenzsystem hat den Hamiltonoperator:

\[\hat{H}_\text{ref} = \sum_i \hat{h}_i^\text{KS}, \qquad \hat{h}^\text{KS} = -\frac{\hbar^2}{2m_e}\nabla^2 + v_\text{eff}(\mathbf{r}) \tag{5.4}\]

Da die Elektronen nicht wechselwirken, ist die Lösung eine Slater-Determinante aus Kohn-Sham-Orbitalen \(\psi_m^\text{KS}\):

\[\rho(\mathbf{r}) = 2\sum_{m=1}^{N_\text{occ}} |\psi_m^\text{KS}(\mathbf{r})|^2 \tag{5.5}\]

7.3.2 Die Energiezerlegung

Die Gesamtenergie des realen Systems wird nun geschrieben als:

\[E[\rho] = T_\text{ref}[\rho] + J[\rho] + E_\text{xc}[\rho] + \int \rho(\mathbf{r})\,v(\mathbf{r})\,d\mathbf{r} \tag{5.6}\]

Term	Bedeutung	Bekannt?
\(T_\text{ref}[\rho]\)	Kinetische Energie des Referenzsystems	✓ aus KS-Orbitalen
\(J[\rho]\)	Klassische Coulomb-Abstossung Elektronen	✓ analytisch
\(\int \rho v\,d\mathbf{r}\)	Kern-Elektron-Anziehung	✓ analytisch
\(E_\text{xc}[\rho]\)	Austausch-Korrelationsenergie	✗ muss approximiert werden

Das Austausch-Korrelationsfunktional \(E_\text{xc}[\rho]\) fängt alles auf, was nicht exakt bekannt ist:

\[E_\text{xc}[\rho] = \underbrace{(T[\rho] - T_\text{ref}[\rho])}_{\text{Korrelation der kin. Energie}} + \underbrace{(V_{ee}[\rho] - J[\rho])}_{\text{Austausch + Korr. der Wechselwirkung}} \tag{5.7}\]

7.3.3 Die Kohn-Sham-Gleichungen

Minimierung von \(E[\rho]\) (Gl. 5.6) über die KS-Orbitale (mit Orthonormalitäts-Nebenbedingungen via Lagrange, Kapitel 00) führt auf das Kohn-Sham-Eigenwertproblem:

\[\hat{h}^\text{KS}\,\psi_m^\text{KS}(\mathbf{r}) = \varepsilon_m^\text{KS}\,\psi_m^\text{KS}(\mathbf{r}) \tag{5.8}\]

mit dem effektiven Einteilchenpotential:

\[v_\text{eff}(\mathbf{r}) = \underbrace{v(\mathbf{r})}_{\text{Kerne}} + \underbrace{e^2\!\int \frac{\rho(\mathbf{r}')}{4\pi\varepsilon_0|\mathbf{r}-\mathbf{r}'|}\,d\mathbf{r}'}_{\text{Hartree-Potential } v_H} + \underbrace{v_\text{xc}(\mathbf{r})}_{\text{xc-Potential}} \tag{5.9}\]

wobei das Austausch-Korrelationspotential die Funktionalableitung ist:

\[v_\text{xc}(\mathbf{r}) = \frac{\delta E_\text{xc}[\rho]}{\delta \rho(\mathbf{r})} \tag{5.10}\]

Genau wie bei HF ist das ein nichtlineares Problem: \(v_\text{eff}\) hängt von \(\rho\) ab, \(\rho\) hängt von den KS-Orbitalen ab, die Eigenfunktionen von \(\hat{h}^\text{KS}(v_\text{eff})\) sind. Deshalb wird auch hier iterativ (SCF) gelöst.

fig, axes = plt.subplots(1, 2, figsize=(11, 5))

for ax, title, terms, col in zip(
    axes,
    ['Hartree-Fock', 'Kohn-Sham DFT'],
    [
        [('$\\hat{h}_1$', 'Kern-Hamilton\n(kin. + Kern-el.)', '#4393c3'),
         ('$\\sum 2\\hat{J}_{m\'} - \\hat{K}_{m\'}$', 'Coulomb + Austausch\n(nichtlokal, orbital-abh.)', '#d6604d'),
         ('$\\hat{f}\\,\\psi_m = \\varepsilon_m\\,\\psi_m$', 'Fock-Gleichung', '#2ca25f')],
        [('$\\hat{h}_1$', 'Kern-Hamilton\n(kin. + Kern-el.)', '#4393c3'),
         ('$v_H(\\mathbf{r}) = e^2\\!\\int\\frac{\\rho}{r_{12}}d\\mathbf{r}\'$',
          'Hartree-Potential\n(klassisch, lokal)', '#f4a582'),
         ('$v_\\mathrm{xc}(\\mathbf{r}) = \\delta E_\\mathrm{xc}/\\delta\\rho$',
          'xc-Potential\n(lokal, approximiert)', '#d6604d'),
         ('$\\hat{h}^\\mathrm{KS}\\psi_m^\\mathrm{KS} = \\varepsilon_m^\\mathrm{KS}\\psi_m^\\mathrm{KS}$',
          'Kohn-Sham-Gleichung', '#2ca25f')],
    ],
    ['#2166ac', '#35978f']
):
    ax.axis('off')
    ax.set_title(title, fontsize=14, fontweight='bold', color=col, pad=12)

    y_pos = 0.85
    for term, desc, tc in terms:
        ax.text(0.15, y_pos, term, ha='center', va='center', fontsize=11,
                color='white', fontweight='bold',
                bbox=dict(boxstyle='round,pad=0.4', facecolor=tc,
                          edgecolor='gray', lw=1),
                transform=ax.transAxes)
        ax.text(0.62, y_pos, desc, ha='left', va='center', fontsize=9.5,
                color='#333', transform=ax.transAxes)
        if y_pos > 0.15:
            ax.annotate('', xy=(0.15, y_pos - 0.14),
                        xytext=(0.15, y_pos - 0.04),
                        xycoords='axes fraction', textcoords='axes fraction',
                        arrowprops=dict(arrowstyle='->', color='#aaa', lw=1.5))
        y_pos -= 0.22

plt.suptitle('Strukturvergleich: HF vs. Kohn-Sham', fontsize=13, y=1.01)
plt.tight_layout()
plt.show()

Abbildung 7.2: Strukturvergleich HF vs. Kohn-Sham. Beide lösen ein SCF-Problem; der Unterschied liegt im effektiven Potential: HF hat nichtlokalen Austausch, KS hat lokales xc-Potential.

---

7.4 Das Austausch-Korrelationsproblem

Das KS-Schema ist exakt – wenn man \(E_\text{xc}[\rho]\) exakt kennte. Man kennt es nicht. Verschiedene Näherungen bilden eine Hierarchie („Jacob’s Ladder”):

7.4.1 Lokale Dichtenäherung (LDA)

Die einfachste Approximation: \(E_\text{xc}\) hängt nur von der Dichte am jeweiligen Punkt ab:

\[E_\text{xc}^\text{LDA}[\rho] = \int \rho(\mathbf{r})\, \varepsilon_\text{xc}(\rho(\mathbf{r}))\,d\mathbf{r} \tag{5.11}\]

\(\varepsilon_\text{xc}(\rho)\) ist die xc-Energie pro Elektron eines homogenen Elektronengases mit Dichte \(\rho\) – bekannt aus Monte-Carlo-Simulationen.

7.4.2 Generalized Gradient Approximation (GGA)

Ergänzt die lokale Dichte um ihren Gradienten:

\[E_\text{xc}^\text{GGA}[\rho] = \int f(\rho(\mathbf{r}), |\nabla\rho(\mathbf{r})|)\,d\mathbf{r} \tag{5.12}\]

Populäre GGA-Funktionale: PBE (Perdew-Burke-Ernzerhof, 1996) – parameterfrei, weit verbreitet in der Festkörperphysik.

7.4.3 Hybrid-Funktionale

Mischen GGA-Austausch mit einem Anteil exaktem HF-Austausch:

\[E_\text{xc}^\text{hybrid} = \alpha\, E_x^\text{HF} + (1-\alpha)\, E_x^\text{GGA} + E_c^\text{GGA} \tag{5.13}\]

Das bekannteste Beispiel ist B3LYP (\(\alpha = 0.2\), mit drei Parametern aus Fits an Atomdaten) – jahrzehntelang das meistverwendete Funktional in der Chemie.

fig, ax = plt.subplots(figsize=(7, 5.5))
ax.axis('off')

rungs = [
    ('LDA', 'Lokale Dichte\n$\\varepsilon_\\mathrm{xc}(\\rho)$', '#ffffcc'),
    ('GGA', 'Gradient\n$f(\\rho, |\\nabla\\rho|)$', '#c7e9b4'),
    ('meta-GGA', 'Kinetische Energiedichte\n$f(\\rho, \\nabla\\rho, \\tau)$', '#7fcdbb'),
    ('Hybrid', 'HF-Austausch-Anteil\nz.B. B3LYP, PBE0', '#41b6c4'),
    ('Double Hybrid', 'MP2-Korrelation\nz.B. B2PLYP', '#1d91c0'),
]

for i, (name, desc, col) in enumerate(rungs):
    y = 0.1 + i * 0.165
    width = 0.55 + i * 0.06
    x = 0.5 - width/2
    rect = plt.Rectangle((x, y), width, 0.13,
                          facecolor=col, edgecolor='#666', lw=1.2,
                          transform=ax.transAxes, clip_on=False)
    ax.add_patch(rect)
    ax.text(0.5, y + 0.065, f'{name}: {desc}',
            ha='center', va='center', fontsize=9.5,
            color='#222', transform=ax.transAxes)

ax.text(0.5, 0.97, "\"Jacob's Ladder\" der Dichtefunktionale",
        ha='center', va='top', fontsize=12, fontweight='bold',
        transform=ax.transAxes)
ax.text(0.92, 0.5, '← höhere\nGenauigkeit',
        ha='center', va='center', fontsize=9, color='#555',
        rotation=90, transform=ax.transAxes)
ax.text(0.08, 0.5, 'chemisches\nPotential →',
        ha='center', va='center', fontsize=9, color='#888',
        rotation=90, transform=ax.transAxes)

ax.annotate('', xy=(0.88, 0.88), xytext=(0.88, 0.12),
            xycoords='axes fraction', textcoords='axes fraction',
            arrowprops=dict(arrowstyle='->', color='#555', lw=2))

ax.set_xlim(0, 1); ax.set_ylim(0, 1)
plt.tight_layout()
plt.show()

Abbildung 7.3: Jacobs Leiter der Dichtefunktionale. Mit jeder Stufe steigt die Genauigkeit (tendenziell), aber auch der Rechenaufwand.

---

7.5 DFT-Rechnung für H₂O mit PySCF

7.5.1 LDA, GGA (PBE) und Hybrid (B3LYP) im Vergleich

mol = gto.Mole()
mol.atom = '''
    O   0.000000   0.000000   0.117176
    H   0.000000   0.756950  -0.468703
    H   0.000000  -0.756950  -0.468703
'''
mol.basis   = '6-31g*'
mol.charge  = 0
mol.spin    = 0
mol.verbose = 0
mol.build()

# HF als Referenz
mf_hf = scf.RHF(mol).run()

results = {'HF': mf_hf.e_tot}
functionals = {
    'LDA (SVWN)': 'lda,vwn',
    'GGA (PBE)':  'pbe,pbe',
    'Hybrid (B3LYP)': 'b3lyp',
}

mf_dft = {}
for label, xc in functionals.items():
    mf = dft.RKS(mol)
    mf.xc      = xc
    mf.verbose = 0
    mf.kernel()
    results[label]  = mf.e_tot
    mf_dft[label]   = mf

print(f"{'Methode':<20}  {'E / Hartree':>14}  {'ΔE vs HF / kcal/mol':>20}")
print('-' * 58)
for k, v in results.items():
    delta = (v - results['HF']) * 627.509
    print(f"  {k:<18}  {v:>14.6f}  {delta:>20.2f}")

Methode                  E / Hartree   ΔE vs HF / kcal/mol
----------------------------------------------------------
  HF                      -76.009132                  0.00
  LDA (SVWN)              -75.840914                105.56
  GGA (PBE)               -76.319767               -194.93
  Hybrid (B3LYP)          -76.406784               -249.53

fig, ax = plt.subplots(figsize=(6, 6))

mf_b3lyp = mf_dft['Hybrid (B3LYP)']
eps_hf   = mf_hf.mo_energy
eps_b3lyp = mf_b3lyp.mo_energy
occ_hf   = mf_hf.mo_occ

n_occ = int(np.sum(occ_hf > 0))
n_show = min(n_occ + 3, len(eps_hf))

lw = 2.5
for i in range(n_show):
    occ = i < n_occ
    col_hf    = '#2166ac' if occ else '#92c5de'
    col_b3lyp = '#d6604d' if occ else '#f4a582'

    ax.hlines(eps_hf[i],    0.1, 0.4, color=col_hf,    lw=lw)
    ax.hlines(eps_b3lyp[i], 0.6, 0.9, color=col_b3lyp, lw=lw)

    ax.plot([0.4, 0.6], [eps_hf[i], eps_b3lyp[i]],
            color='#aaa', lw=0.8, ls='--')

ax.text(0.25, 0.02, 'HF',     ha='center', fontsize=12, fontweight='bold',
        color='#2166ac', transform=ax.transAxes)
ax.text(0.75, 0.02, 'B3LYP', ha='center', fontsize=12, fontweight='bold',
        color='#d6604d', transform=ax.transAxes)

ax.axhline(0, color='gray', lw=0.8, ls=':', alpha=0.5)
ax.set_xlim(0, 1)
ax.set_xticks([])
ax.set_ylabel('Orbitalenergie / Hartree', fontsize=12)
ax.set_title('H₂O Orbitalenergien: HF vs. B3LYP (6-31G*)', fontsize=12)
ax.grid(True, axis='y', alpha=0.25)

from matplotlib.lines import Line2D
ax.legend(handles=[
    Line2D([0],[0], color='#2166ac', lw=2.5, label='HF besetzt'),
    Line2D([0],[0], color='#92c5de', lw=2.5, label='HF virtuell'),
    Line2D([0],[0], color='#d6604d', lw=2.5, label='B3LYP besetzt'),
    Line2D([0],[0], color='#f4a582', lw=2.5, label='B3LYP virtuell'),
], fontsize=9, loc='lower right')
plt.tight_layout()
plt.show()

Abbildung 7.4: Orbitalenergien von H₂O: HF vs. B3LYP (6-31G*). B3LYP-Orbitalenergien liegen systematisch höher (weniger negativ) als HF-Orbitalenergien.

7.5.2 Elektronendichte visualisieren

from pyscf.tools import cubegen
import tempfile, os

# Dichtematrix aus B3LYP
dm_b3lyp = mf_b3lyp.make_rdm1()

# Gitter in der Molekülebene (xz-Ebene, y=0)
nx, nz = 120, 100
x_grid = np.linspace(-3.0, 3.0, nx)
z_grid = np.linspace(-2.5, 2.5, nz)
X, Z   = np.meshgrid(x_grid, z_grid)
Y      = np.zeros_like(X)

coords = np.column_stack([X.ravel(), Y.ravel(), Z.ravel()])

# Auswertung der Dichte auf dem Gitter via PySCF ao_value
from pyscf import df
from pyscf.dft import numint

ao_vals = numint.eval_ao(mol, coords)                    # AO-Werte auf Gitter
rho     = numint.eval_rho(mol, ao_vals, dm_b3lyp)        # Dichte = Σ P_μν φ_μ φ_ν
RHO     = rho.reshape(nz, nx)

# Kernpositionen (in Bohr)
atm_coords = mol.atom_coords()   # Angström → Bohr bereits in PySCF
atm_syms   = [mol.atom_symbol(i) for i in range(mol.natm)]

fig, ax = plt.subplots(figsize=(7, 6))
levels  = np.array([0.001, 0.01, 0.05, 0.1, 0.2, 0.4, 0.8, 1.5])
cs = ax.contourf(X, Z, RHO, levels=levels,
                 cmap='Blues', extend='max')
ax.contour(X, Z, RHO, levels=levels, colors='navy', linewidths=0.5, alpha=0.4)
cb = plt.colorbar(cs, ax=ax)
cb.set_label('$\\rho(\\mathbf{r})$ / $e\\,a_0^{-3}$', fontsize=11)

# Atome einzeichnen (Koordinaten in Bohr → Angström ×0.529)
cols_atm = {'O': '#d6604d', 'H': '#4393c3'}
for sym, c in zip(atm_syms, atm_coords):
    ax.plot(c[0]*0.529, c[2]*0.529, 'o',
            color=cols_atm[sym], ms=14, zorder=5,
            markeredgecolor='white', markeredgewidth=1.5)
    ax.text(c[0]*0.529 + 0.08, c[2]*0.529 + 0.08, sym,
            fontsize=11, fontweight='bold', color=cols_atm[sym])

ax.set_xlabel('$x$ / Å', fontsize=12)
ax.set_ylabel('$z$ / Å', fontsize=12)
ax.set_title('Elektronendichte H₂O (B3LYP/6-31G*),\nSchnitt durch die Molekülebene',
             fontsize=12)
ax.set_aspect('equal')
plt.tight_layout()
plt.show()

Abbildung 7.5: Elektronendichte ρ(r) von H₂O (B3LYP/6-31G*) in der Molekülebene. Die Dichte konzentriert sich stark am Sauerstoffkern; die Bindungen sind als erhöhte Dichte zwischen den Atomen sichtbar.

7.5.3 HOMO und LUMO

mf_use  = mf_b3lyp
C       = mf_use.mo_coeff
n_occ_b = int(np.sum(mf_use.mo_occ > 0))

fig, axes = plt.subplots(1, 2, figsize=(11, 5))
titles    = [f'HOMO (MO {n_occ_b})', f'LUMO (MO {n_occ_b+1})']
mo_idxs   = [n_occ_b - 1, n_occ_b]

for ax, mo_idx, title in zip(axes, mo_idxs, titles):
    # MO-Amplitude auf Gitter
    ao_v  = numint.eval_ao(mol, coords)
    mo_v  = ao_v @ C[:, mo_idx]
    MO    = mo_v.reshape(nz, nx)

    vmax = np.percentile(np.abs(MO), 99)
    im = ax.contourf(X, Z, MO,
                     levels=np.linspace(-vmax, vmax, 21),
                     cmap='RdBu_r')
    ax.contour(X, Z, MO, levels=[0], colors='k', linewidths=1.0)
    plt.colorbar(im, ax=ax)

    for sym, c in zip(atm_syms, atm_coords):
        ax.plot(c[0]*0.529, c[2]*0.529, 'o',
                color=cols_atm[sym], ms=12, zorder=5,
                markeredgecolor='white', markeredgewidth=1.5)
        ax.text(c[0]*0.529+0.07, c[2]*0.529+0.07, sym,
                fontsize=10, fontweight='bold', color=cols_atm[sym])

    ax.set_xlabel('$x$ / Å', fontsize=11)
    ax.set_ylabel('$z$ / Å', fontsize=11)
    ax.set_title(title, fontsize=12)
    ax.set_aspect('equal')

plt.suptitle('H₂O Molekülorbitale (B3LYP/6-31G*)', fontsize=13)
plt.tight_layout()
plt.show()

Abbildung 7.6: HOMO (links) und LUMO (rechts) von H₂O aus der B3LYP/6-31G*-Rechnung. Das HOMO entspricht einem freien Elektronenpaar am Sauerstoff; das LUMO hat antibindenden Charakter.

---

7.6 HF vs. DFT: konzeptueller Vergleich

methods  = list(results.keys())
energies = [results[k] for k in methods]
ref      = results['HF']
delta_E  = [(e - ref)*627.509 for e in energies]

colors_bar = ['#4393c3', '#fdae61', '#a6d96a', '#1a9641']
fig, ax = plt.subplots(figsize=(7, 4))
bars = ax.bar(methods, delta_E, color=colors_bar,
              edgecolor='white', lw=1.5)
ax.axhline(0, color='gray', lw=1.0, ls='--')
for bar, val in zip(bars, delta_E):
    va  = 'bottom' if val >= 0 else 'top'
    off = 0.3 if val >= 0 else -0.3
    ax.text(bar.get_x()+bar.get_width()/2, val+off,
            f'{val:+.1f}', ha='center', va=va, fontsize=10,
            fontweight='bold', color='#333')
ax.set_ylabel('$\\Delta E$ relativ zu HF / kcal/mol', fontsize=12)
ax.set_title('Energievergleich HF vs. DFT, H₂O (6-31G*)', fontsize=12)
ax.grid(True, axis='y', alpha=0.3)
plt.tight_layout()
plt.show()

Abbildung 7.7: Energievergleich HF vs. verschiedene DFT-Funktionale für H₂O (6-31G*). Alle Methoden lösen ein SCF-Problem; der Unterschied liegt in der Behandlung von Austausch und Korrelation.

HinweisWarum DFT tiefere Energien liefert

LDA und GGA schliessen Korrelationsenergie ein – die HF per Definition fehlt. Daher liegen DFT-Energien tiefer als HF-Energien. Das bedeutet aber nicht automatisch, dass DFT genauer ist: Das Variationsprinzip gilt für die exakte Wellenfunktion, nicht für approximate Dichten. DFT-Energien können die wahre Energie auch unterschreiten.

---

7.7 Zusammenfassung: HF vs. KS-DFT

	Hartree-Fock	Kohn-Sham DFT
Grundgrösse	Wellenfunktion \(\Psi\) (\(3N_e\)-dim.)	Elektronendichte \(\rho(\mathbf{r})\) (3-dim.)
Näherung	Einzige Slater-Determinante	Einzige Slater-Det. des Referenzsystems
Austausch	Exakt (nichtlokal, \(\hat{K}\))	Approximiert in \(E_\text{xc}[\rho]\)
Korrelation	Fehlt vollständig	Approximiert in \(E_\text{xc}[\rho]\)
Skalierung	\(O(N_b^4)\)	\(O(N_b^3)\) (formell)
Genauigkeit	Systematisch verbesserbar (post-HF)	Abhängig vom Funktional
SCF-Gleichung	\(\mathbf{F}\mathbf{c} = \mathbf{S}\mathbf{c}\boldsymbol{\varepsilon}\)	Gleiche Struktur, anderes \(\mathbf{F}\)

Das Letzte in der Tabelle ist entscheidend: Die Roothaan-Gleichungen aus Kapitel 03 gelten für DFT unverändert – nur die Fock-Matrix \(\mathbf{F}\) enthält statt des nichtlokalen Austauschoperators \(\hat{K}\) das lokale xc-Potential \(v_\text{xc}(\mathbf{r})\). Der SCF-Algorithmus ist identisch.