import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh
# Einheitliches Plot-Styling
plt.rcParams.update({'figure.dpi': 120, 'font.size': 11})8 Ü1 · QM Grundlagen
Wellenfunktionen, Operatoren und die Schrödinger-Gleichung
8.1 Lernziele
Nach dieser Übung kannst du:
- Wellenfunktionen für einfache Systeme (Teilchen im Kasten) numerisch berechnen und visualisieren
- Erwartungswerte von Observablen (Ort, Impuls, Energie) berechnen
- Die Heisenbergsche Unschärferelation numerisch überprüfen
- Eine erste Hartree-Fock-Rechnung mit PySCF ausführen und die Ausgabe interpretieren
8.2 Aufgabe 1 · Teilchen im Kasten (Partikel-in-einer-Box)
Das eindimensionale Teilchen-im-Kasten-Modell (Länge \(L\), unendlich hohe Wände) hat exakte analytische Lösungen:
\[\psi_n(x) = \sqrt{\frac{2}{L}} \sin\!\left(\frac{n\pi x}{L}\right), \quad E_n = \frac{n^2 \pi^2 \hbar^2}{2mL^2}\]
8.2.1 1a) Wellenfunktionen visualisieren
# Parameter (atomare Einheiten: hbar = m_e = 1)
L = 10.0 # Kastenlänge in Bohr
N_grid = 500 # Gitterpunkte
x = np.linspace(0, L, N_grid)
def psi_box(n, x, L):
"""Wellenfunktion des Teilchens im Kasten für Quantenzahl n."""
return np.sqrt(2 / L) * np.sin(n * np.pi * x / L)
def E_box(n, L):
"""Energie in atomaren Einheiten (hbar=m=1)."""
return (n * np.pi)**2 / (2 * L**2)
# Plotte n = 1, 2, 3, 4
fig, axes = plt.subplots(2, 2, figsize=(10, 7))
axes = axes.ravel()
for i, n in enumerate([1, 2, 3, 4]):
psi = psi_box(n, x, L)
prob = psi**2
En = E_box(n, L)
axes[i].plot(x, psi, label=r'$\psi_n$', color='steelblue', lw=2)
axes[i].plot(x, prob, label=r'$|\psi_n|^2$', color='darkorange', lw=2, ls='--')
axes[i].set_title(f'n = {n}, E = {En:.4f} Hartree')
axes[i].set_xlabel('x / Bohr')
axes[i].legend(fontsize=9)
axes[i].axhline(0, color='k', lw=0.5)
plt.suptitle('Wellenfunktionen und Wahrscheinlichkeitsdichten: Teilchen im Kasten', fontsize=13)
plt.tight_layout()
plt.show()
8.2.2 1b) Normierung überprüfen
Aufgabe: Überprüfe numerisch, dass \(\langle\psi_n|\psi_n\rangle = 1\) für \(n = 1, 2, 3\). Berechne ausserdem \(\langle\psi_1|\psi_2\rangle\) und interpretiere das Ergebnis.
dx = x[1] - x[0] # Gitterabstand
print("Normierungsprüfung (numerische Integration):")
for n in [1, 2, 3]:
psi = psi_box(n, x, L)
norm = np.trapezoid(psi**2, x)
print(f" <psi_{n}|psi_{n}> = {norm:.8f}")
print()
print("Orthogonalität:")
for (m, n) in [(1, 2), (1, 3), (2, 3)]:
overlap = np.trapezoid(psi_box(m, x, L) * psi_box(n, x, L), x)
print(f" <psi_{m}|psi_{n}> = {overlap:.2e}")Normierungsprüfung (numerische Integration):
<psi_1|psi_1> = 1.00000000
<psi_2|psi_2> = 1.00000000
<psi_3|psi_3> = 1.00000000
Orthogonalität:
<psi_1|psi_2> = 0.00e+00
<psi_1|psi_3> = 5.55e-17
<psi_2|psi_3> = -5.55e-178.2.3 1c) Erwartungswerte
Aufgabe: Berechne \(\langle x \rangle\), \(\langle x^2 \rangle\) und \(\sigma_x = \sqrt{\langle x^2\rangle - \langle x\rangle^2}\) für \(n = 1, 2, 3\).
print(f"{'n':>3} | {'<x>':>10} | {'<x²>':>10} | {'σ_x':>10}")
print("-" * 42)
for n in [1, 2, 3, 4]:
psi = psi_box(n, x, L)
prob = psi**2
x_mean = np.trapezoid(prob * x, x)
x2_mean = np.trapezoid(prob * x**2, x)
sigma_x = np.sqrt(x2_mean - x_mean**2)
print(f"{n:>3} | {x_mean:>10.4f} | {x2_mean:>10.4f} | {sigma_x:>10.4f}")
print()
print(f"Analytisch: <x> = L/2 = {L/2:.4f}") n | <x> | <x²> | σ_x
------------------------------------------
1 | 5.0000 | 28.2673 | 1.8076
2 | 5.0000 | 32.0668 | 2.6583
3 | 5.0000 | 32.7704 | 2.7876
4 | 5.0000 | 33.0167 | 2.8314
Analytisch: <x> = L/2 = 5.00008.3 Aufgabe 2 · Numerische Lösung der Schrödinger-Gleichung
Löse die zeitunabhängige Schrödinger-Gleichung für ein harmonisches Potential \(V(x) = \frac{1}{2}kx^2\) numerisch via Finite-Differenzen-Diskretisierung.
Der diskrete kinetische Energie-Operator lautet: \[T_{ij} = -\frac{\hbar^2}{2m}\cdot\frac{-2\delta_{ij} + \delta_{i,j-1} + \delta_{i,j+1}}{(\Delta x)^2}\]
def solve_schrodinger_1d(V_func, x_min=-8, x_max=8, N=600, n_states=6):
"""
Löst die 1D stationäre Schrödinger-Gleichung via Finite Differenzen.
Atomare Einheiten: hbar = m = 1.
Returns: x, energies (n_states lowest), wavefunctions (columns)
"""
x = np.linspace(x_min, x_max, N)
dx = x[1] - x[0]
# Kinetische Energie: Tridiagonal-Matrix
diag = np.ones(N) / dx**2 # Hauptdiagonale
off_diag = -0.5 * np.ones(N-1) / dx**2 # Nebendiagonalen
T = np.diag(diag) + np.diag(off_diag, 1) + np.diag(off_diag, -1)
# Potentialmatrix (diagonal)
V = np.diag(V_func(x))
# Hamilton-Operator H = T + V
H = T + V
# Eigenwertproblem lösen (symmetrisch)
energies, psi = eigh(H, subset_by_index=[0, n_states - 1])
# Normieren
for i in range(n_states):
psi[:, i] /= np.sqrt(np.trapezoid(psi[:, i]**2, x))
return x, energies, psi
# Harmonisches Potential (k=1, analytische Lösung: E_n = n + 1/2)
k = 1.0
x, E, psi = solve_schrodinger_1d(lambda x: 0.5 * k * x**2)
print("Harmonischer Oszillator – Energievergleich:")
print(f"{'n':>3} | {'E_num':>10} | {'E_exakt':>10} | {'Fehler':>10}")
print("-" * 42)
for n in range(6):
E_exact = n + 0.5 # in atomaren Einheiten
print(f"{n:>3} | {E[n]:>10.6f} | {E_exact:>10.6f} | {abs(E[n]-E_exact):>10.2e}")Harmonischer Oszillator – Energievergleich:
n | E_num | E_exakt | Fehler
------------------------------------------
0 | 0.499978 | 0.500000 | 2.23e-05
1 | 1.499889 | 1.500000 | 1.11e-04
2 | 2.499710 | 2.500000 | 2.90e-04
3 | 3.499442 | 3.500000 | 5.58e-04
4 | 4.499086 | 4.500000 | 9.14e-04
5 | 5.498640 | 5.500000 | 1.36e-03# Visualisierung
fig, ax = plt.subplots(figsize=(9, 6))
colors = plt.cm.viridis(np.linspace(0.1, 0.9, 6))
V_plot = 0.5 * k * x**2
ax.plot(x, V_plot, 'k-', lw=2, label='$V(x) = \\frac{1}{2}x^2$')
for n in range(5):
scale = 0.8 # Skalierung für Darstellung
ax.axhline(E[n], color=colors[n], ls=':', lw=1, alpha=0.5)
ax.plot(x, scale * psi[:, n] + E[n], color=colors[n], lw=2,
label=f'n={n}, E={E[n]:.3f}')
ax.set_xlim(-6, 6)
ax.set_ylim(-0.3, 5.5)
ax.set_xlabel('x / Bohr')
ax.set_ylabel('Energie / Hartree')
ax.set_title('Harmonischer Oszillator: numerische Eigenzustände')
ax.legend(fontsize=9, loc='upper center')
plt.tight_layout()
plt.show()
8.4 Aufgabe 3 · Heisenbergsche Unschärferelation
Die Heisenbergsche Unschärferelation lautet: \(\sigma_x \cdot \sigma_p \geq \frac{\hbar}{2}\)
Der Impuls-Operator in Ortsdarstellung ist \(\hat{p} = -i\hbar\,\frac{d}{dx}\), dessen Erwartungswert via numerische Ableitung berechnet wird.
Aufgabe: Überprüfe die Unschärferelation für die ersten 4 Zustände des harmonischen Oszillators. Der Grundzustand des harmonischen Oszillators erfüllt die Gleichung exakt – überprüfe das!
dx = x[1] - x[0]
print("Heisenbergsche Unschärferelation (hbar = 1):")
print(f"{'n':>3} | {'σ_x':>8} | {'σ_p':>8} | {'σ_x·σ_p':>10} | {'≥ 0.5?':>8}")
print("-" * 50)
for n in range(4):
p = psi[:, n]
prob = p**2
# Ortsunschärfe
x_mean = np.trapezoid(prob * x, x)
x2_mean = np.trapezoid(prob * x**2, x)
sigma_x = np.sqrt(x2_mean - x_mean**2)
# Impulsunschärfe: numerische Ableitung (finite Differenz)
dpsi = np.gradient(p, dx) # d/dx psi
d2psi = np.gradient(dpsi, dx) # d²/dx² psi
# <p²> = integral psi* (-d²/dx²) psi dx (atomare Einheiten)
p2_mean = -np.trapezoid(p * d2psi, x)
# <p> = integral psi* (-i d/dx) psi dx -> imaginär, für Reelfkt. = 0
p_mean = 0.0 # Reelle Eigenfunktionen: <p> = 0
sigma_p = np.sqrt(max(p2_mean - p_mean**2, 0))
product = sigma_x * sigma_p
ok = "✓" if product >= 0.499 else "✗"
print(f"{n:>3} | {sigma_x:>8.4f} | {sigma_p:>8.4f} | {product:>10.4f} | {ok:>8}")
print()
print("Analytisch (harm. Osz.): σ_x·σ_p = (n + 1/2)·ℏ")
for n in range(4):
print(f" n={n}: {n + 0.5:.1f}")Heisenbergsche Unschärferelation (hbar = 1):
n | σ_x | σ_p | σ_x·σ_p | ≥ 0.5?
--------------------------------------------------
0 | 0.7071 | 0.7070 | 0.4999 | ✓
1 | 1.2247 | 1.2245 | 1.4996 | ✓
2 | 1.5810 | 1.5806 | 2.4988 | ✓
3 | 1.8705 | 1.8699 | 3.4978 | ✓
Analytisch (harm. Osz.): σ_x·σ_p = (n + 1/2)·ℏ
n=0: 0.5
n=1: 1.5
n=2: 2.5
n=3: 3.58.5 Aufgabe 4 · Erste PySCF-Rechnung: H₂O
Führe eine Hartree-Fock-Rechnung für H₂O durch und interpretiere die Ausgabe.
Voraussetzung:
pyscfmuss installiert sein (uv run pip install pyscfim Projektverzeichnis).
try:
from pyscf import gto, scf
PYSCF_AVAILABLE = True
except ImportError:
PYSCF_AVAILABLE = False
print("PySCF nicht gefunden. Installiere mit: uv add pyscf")
if PYSCF_AVAILABLE:
# Molekül definieren (experimentelle Geometrie, Angström)
mol = gto.Mole()
mol.atom = """
O 0.000 0.000 0.000
H 0.757 0.586 0.000
H -0.757 0.586 0.000
"""
mol.basis = 'sto-3g'
mol.verbose = 3
mol.build()
# RHF-Rechnung
mf = scf.RHF(mol)
E_hf = mf.kernel()converged SCF energy = -74.9629466565388if PYSCF_AVAILABLE:
import pandas as pd
mo_energies = mf.mo_energy # in Hartree
mo_occ = mf.mo_occ
df = pd.DataFrame({
'MO': range(1, len(mo_energies) + 1),
'E [Hartree]': mo_energies,
'E [eV]': mo_energies * 27.211,
'Besetzt': mo_occ.astype(int)
})
print(f"\nHF-Energie (STO-3G): {E_hf:.6f} Hartree")
print(f"HOMO-Energie: {mo_energies[mo_occ > 0][-1]*27.211:.2f} eV")
print(f"LUMO-Energie: {mo_energies[mo_occ == 0][0]*27.211:.2f} eV")
print(f"HOMO-LUMO-Gap: {(mo_energies[mo_occ == 0][0] - mo_energies[mo_occ > 0][-1])*27.211:.2f} eV")
print()
print(df.to_string(index=False))
HF-Energie (STO-3G): -74.962947 Hartree
HOMO-Energie: -10.65 eV
LUMO-Energie: 16.48 eV
HOMO-LUMO-Gap: 27.12 eV
MO E [Hartree] E [eV] Besetzt
1 -20.241762 -550.798579 2
2 -1.268360 -34.513356 2
3 -0.617863 -16.812664 2
4 -0.452999 -12.326557 2
5 -0.391242 -10.646096 2
6 0.605576 16.478334 0
7 0.742245 20.197237 08.5.1 Reflexionsfragen
- Warum sind alle MO-Energien für besetzte Orbitale negativ? Was bedeutet das physikalisch?
- Das HOMO von H₂O liegt bei ca. −12.6 eV (experimentell: erstes Ionisierungspotential). Wie gut ist STO-3G?
- Wieso spricht man beim Koopmans-Theorem von einer Näherung? Was wird vernachlässigt?
8.6 Zusammenfassung
In dieser Übung hast du:
- Wellenfunktionen für das Teilchen im Kasten und den harmonischen Oszillator berechnet und visualisiert
- Erwartungswerte numerisch via numerischer Integration bestimmt
- Die Heisenbergsche Unschärferelation am Beispiel des harmonischen Oszillators überprüft
- Eine erste HF/STO-3G-Rechnung mit PySCF für H₂O durchgeführt