5  Die Hartree-Fock-Methode

In den vorangegangenen Kapiteln haben wir das Problem klar formuliert: Wir suchen die Grundzustandswellenfunktion und -energie des elektronischen Hamiltonoperators für H₂O. Die Wellenfunktion muss antisymmetrisch sein – als Slater-Determinante ausgedrückt. Kapitel 03 zeigt, wie man die optimalen Orbitale in dieser Determinante findet.

Die Antwort ist das Hartree-Fock-Verfahren (HF): ein iterativer Algorithmus, der das Vielteilchenproblem in eine Folge von Einelektronenproblemen zerlegt, und dessen Kern – nach dem LCAO-Ansatz – ein verallgemeinertes Matrixeigenwertproblem ist, das du aus der linearen Algebra kennst.

---

5.1 Die Mean-field-Idee

Der Grund, warum die Schrödinger-Gleichung für Mehrelektronensysteme nicht exakt lösbar ist, ist der Zweielectronenterm:

\[\hat{V}_{ee} = \frac{1}{2}\sum_{i \neq j} \frac{e^2}{4\pi\varepsilon_0\, r_{ij}}\]

Er koppelt die Bewegung aller Elektronen gleichzeitig – man kann die Wellenfunktion nicht in ein Produkt von Einelektronenfunktionen separieren.

5.1.1 Die Näherung: Jedes Elektron im mittleren Feld

Der Ausweg: Man ersetzt die instantane, fluktuierende Wechselwirkung durch ein mittleres, effektives Potential. Statt zu fragen „Wie wechselwirkt Elektron 1 gerade mit Elektron 2?“, fragt man: „In welchem gemittelten elektrostatischen Feld der anderen Elektronen bewegt sich Elektron 1?”

Das ist die Mean-field-Näherung: Jedes Elektron bewegt sich unabhängig in einem effektiven Ein-Elektron-Potential, das die gemittelte Wirkung aller anderen Elektronen kodiert.

HinweisAnalogie aus dem ML

Die Mean-field-Idee ist konzeptuell verwandt mit Mean-field-Methoden in der Statistik: Man approximiert eine komplexe, hochdimensionale Verteilung über viele Variablen durch ein Produkt von Einzelverteilungen, optimiert diese Einzelverteilungen iterativ aneinander an.

---

5.2 Der Fock-Operator

Das Variationsprinzip fordert: Minimiere \(\langle E \rangle = \langle \Psi | \hat{H} | \Psi \rangle\) über alle normierten Slater-Determinanten \(\Psi = |\phi_1\cdots\phi_N\rangle\).

Das führt (nach Lagrange-Multiplikatoren für die Orthonormalitätsbedingungen, Details in Anhang A) auf eine Menge von gekoppelten Einelektronen-Gleichungen, eine pro besetztem Orbital:

\[\hat{f}\, \psi_m(\mathbf{r}) = \varepsilon_m\, \psi_m(\mathbf{r}) \qquad m = 1, \ldots, N_\text{occ} \tag{3.1}\]

Das ist wieder ein Eigenwertproblem – aber für den Fock-Operator \(\hat{f}\), nicht für den Hamiltonoperator \(\hat{H}\).

5.2.1 Aufbau des Fock-Operators

Der Fock-Operator besteht aus drei Teilen:

\[\hat{f} = \underbrace{\hat{h}}_{\text{Kern}} + \underbrace{\sum_{m'} 2\hat{J}_{m'}}_{\text{Coulomb}} - \underbrace{\sum_{m'} \hat{K}_{m'}}_{\text{Austausch}} \tag{3.2}\]

1. Der Kern-Hamiltonoperator \(\hat{h}\):

\[\hat{h}_1 = -\frac{\hbar^2}{2m_e}\nabla_1^2 - \sum_I \frac{Z_I e^2}{4\pi\varepsilon_0\, r_{1I}}\]

Er beschreibt die kinetische Energie von Elektron 1 und seine Anziehung durch alle Kerne. Er hängt nicht von den anderen Elektronen ab.

2. Der Coulomb-Operator \(\hat{J}_{m'}\) (klassisch):

\[\hat{J}_{m'}(\mathbf{r}_1)\,\psi_m(\mathbf{r}_1) = \left[\int \psi_{m'}^*(\mathbf{r}_2)\,\frac{e^2}{4\pi\varepsilon_0\, r_{12}}\, \psi_{m'}(\mathbf{r}_2)\,d\mathbf{r}_2\right] \psi_m(\mathbf{r}_1) \tag{3.3}\]

Die eckige Klammer ist ein skalares Potential – die klassische elektrostatische Energie, die Elektron 1 durch die gemittelte Ladungsverteilung von Elektron 2 im Orbital \(\psi_{m'}\) spürt.

3. Der Austausch-Operator \(\hat{K}_{m'}\) (nichtklassisch):

\[\hat{K}_{m'}(\mathbf{r}_1)\,\psi_m(\mathbf{r}_1) = \left[\int \psi_{m'}^*(\mathbf{r}_2)\,\frac{e^2}{4\pi\varepsilon_0\, r_{12}}\, \psi_m(\mathbf{r}_2)\,d\mathbf{r}_2\right] \psi_{m'}(\mathbf{r}_1) \tag{3.4}\]

Er hat kein klassisches Analogon. Er entsteht direkt aus der Antisymmetrie der Slater-Determinante (Spinkorrelation aus Kapitel 02). Das Bemerkenswerte: Im Ergebnis vertauscht \(\hat{K}_{m'}\) die Orbitale \(\psi_m\) und \(\psi_{m'}\) im Integranden.

WichtigKernpunkt

Der Fock-Operator \(\hat{f}\) hängt selbst von seinen Eigenfunktionen \(\{\psi_{m'}\}\) ab – über \(\hat{J}_{m'}\) und \(\hat{K}_{m'}\). Das macht Gl. (3.1) zu einem nichtlinearen Eigenwertproblem, das nur iterativ gelöst werden kann.

5.2.2 Gesamtenergie aus den Orbitalenergien

Wichtig: Die Summe der Orbitalenergien \(\sum_m \varepsilon_m\) ist nicht die Gesamtenergie des Systems, weil die Elektron-Elektron-Wechselwirkungen sonst doppelt gezählt würden. Die korrekte HF-Gesamtenergie ist:

\[E_\text{HF} = 2\sum_{m=1}^{N_\text{occ}} \varepsilon_m - \sum_{m,m'} (2J_{mm'} - K_{mm'}) \tag{3.5}\]

---

5.3 Von Orbitalen zu Matrizen: Die Roothaan-Gleichungen

5.3.1 Der LCAO-Ansatz

Die HF-Gleichungen (3.1) sind Differentialgleichungen in \(\mathbb{R}^3\) – direkt numerisch zu lösen ist für Moleküle zu aufwändig. Der Ausweg: Man entwickelt die Molekülorbitale \(\psi_m\) in einer endlichen Basis von Atomorbitalen \(\{\chi_\mu\}\) (LCAO = Linear Combination of Atomic Orbitals):

\[\psi_m(\mathbf{r}) = \sum_{\mu=1}^{N_b} c_{\mu m}\, \chi_\mu(\mathbf{r}) \tag{3.6}\]

Die \(N_b\) Basisfunktionen \(\chi_\mu\) sind bekannte, fest gewählte Funktionen (Details in Kapitel 04). Die unbekannten Koeffizienten \(c_{\mu m}\) sind das, was wir bestimmen wollen. Das Eigenwertproblem im Funktionenraum wird dadurch zu einem endlichdimensionalen Matrixproblem.

5.3.2 Die Fock- und Überlappmatrix

Setzt man (3.6) in (3.1) ein, multipliziert von links mit \(\chi_\nu^*\) und integriert, erhält man:

\[\sum_\mu F_{\nu\mu}\, c_{\mu m} = \varepsilon_m \sum_\mu S_{\nu\mu}\, c_{\mu m} \tag{3.7}\]

In Matrixform – die Roothaan-Gleichungen:

\[\mathbf{F}\,\mathbf{c} = \mathbf{S}\,\mathbf{c}\,\boldsymbol{\varepsilon} \tag{3.8}\]

Matrix	Definition	Bedeutung
\(\mathbf{F}\)	\(F_{\nu\mu} = \int \chi_\nu^*\, \hat{f}\, \chi_\mu\, d\mathbf{r}\)	Fock-Matrix
\(\mathbf{S}\)	\(S_{\nu\mu} = \int \chi_\nu^*\, \chi_\mu\, d\mathbf{r}\)	Überlapp-Matrix
\(\mathbf{c}\)	\((N_b \times N_b)\)-Matrix der Koeffizienten \(c_{\mu m}\)	MO-Koeffizienten
\(\boldsymbol{\varepsilon}\)	Diagonalmatrix der \(\varepsilon_m\)	Orbitalenergien

Das ist ein verallgemeinertes Eigenwertproblem \(\mathbf{F}\mathbf{c} = \mathbf{S}\mathbf{c}\boldsymbol{\varepsilon}\). Im Unterschied zum Standard-Eigenwertproblem steht rechts nicht \(\mathbf{c}\boldsymbol{\varepsilon}\), sondern \(\mathbf{S}\mathbf{c}\boldsymbol{\varepsilon}\) – weil die Basis nicht orthogonal ist (\(\mathbf{S} \neq \mathbf{I}\)).

5.3.3 Die Fock-Matrix explizit

Die Matrixelemente der Fock-Matrix sind:

\[F_{\nu\mu} = h_{\nu\mu} + \sum_{\lambda,\sigma} P_{\lambda\sigma} \left[(\nu\mu|\sigma\lambda) - \frac{1}{2}(\nu\lambda|\sigma\mu)\right] \tag{3.9}\]

mit der Dichtematrix:

\[P_{\lambda\sigma} = 2\sum_{m=1}^{N_\text{occ}} c_{\lambda m}^*\, c_{\sigma m} \tag{3.10}\]

und den Zweielectronenintegralen:

\[(\nu\mu|\sigma\lambda) = \int\!\int \chi_\nu^*(\mathbf{r}_1)\,\chi_\mu(\mathbf{r}_1)\, \frac{e^2}{4\pi\varepsilon_0\, r_{12}}\, \chi_\sigma^*(\mathbf{r}_2)\,\chi_\lambda(\mathbf{r}_2)\, d\mathbf{r}_1\,d\mathbf{r}_2 \tag{3.11}\]

Die Zahl dieser Integrale wächst wie \(O(N_b^4)\) – für einen Basissatz mit \(N_b = 100\) Funktionen sind das bis zu \(10^8\) Integrale. Das ist der dominante Rechenkostenpunkt einer HF-Rechnung.

---

5.4 Der SCF-Algorithmus

Da \(\mathbf{F}\) von \(\mathbf{P}\) abhängt und \(\mathbf{P}\) von den Koeffizienten \(\mathbf{c}\) – die wiederum Eigenvektoren von \(\mathbf{F}\) sind – muss man iterativ vorgehen:

SCF-Algorithmus:
─────────────────────────────────────────────────────────
1. INIT    Wähle Startkoeffizienten c⁽⁰⁾ (z.B. aus Hückel)
2. DICHTE  Bilde P⁽ⁿ⁾ = 2 Σ_m c^(n)* c^(n)   [Gl. 3.10]
3. FOCK    Berechne F⁽ⁿ⁾(P⁽ⁿ⁾)               [Gl. 3.9]
4. DIAG    Löse F⁽ⁿ⁾ c⁽ⁿ⁺¹⁾ = S c⁽ⁿ⁺¹⁾ ε⁽ⁿ⁺¹⁾  [Gl. 3.8]
5. CHECK   ||P⁽ⁿ⁺¹⁾ - P⁽ⁿ⁾|| < Schwellenwert?
           → Ja:  KONVERGIERT → gib E_HF, ε, c aus
           → Nein: n ← n+1, weiter zu Schritt 2
─────────────────────────────────────────────────────────

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

boxes = [
    (0.5, 0.92, 'START\nBasissatz + Geometrie', '#4393c3', 'white'),
    (0.5, 0.76, 'Startkoeffizienten $\\mathbf{c}^{(0)}$\n(Hückel / Core-Hamiltonian)', '#92c5de', 'black'),
    (0.5, 0.60, 'Dichtematrix $\\mathbf{P}^{(n)}$\n$P_{\\lambda\\sigma} = 2\\sum_m c_{\\lambda m}^* c_{\\sigma m}$', '#f4a582', 'black'),
    (0.5, 0.44, 'Fock-Matrix $\\mathbf{F}^{(n)}(\\mathbf{P}^{(n)})$\n$F_{\\nu\\mu} = h_{\\nu\\mu} + \\sum P[2(\\nu\\mu|\\sigma\\lambda)-(\\nu\\lambda|\\sigma\\mu)]$', '#f4a582', 'black'),
    (0.5, 0.28, 'Verallg. Eigenwertproblem\n$\\mathbf{F}\\mathbf{c} = \\mathbf{S}\\mathbf{c}\\boldsymbol{\\varepsilon}$', '#f4a582', 'black'),
    (0.5, 0.12, 'KONVERGIERT\n$E_\\mathrm{HF}$, $\\varepsilon_m$, $\\mathbf{c}$', '#4dac26', 'white'),
]

for x, y, txt, col, tcol in boxes:
    ax.text(x, y, txt, ha='center', va='center', fontsize=10,
            color=tcol, fontweight='bold' if col in ['#4393c3','#4dac26'] else 'normal',
            bbox=dict(boxstyle='round,pad=0.5', facecolor=col,
                      edgecolor='gray', lw=1.2),
            transform=ax.transAxes)

# Pfeile
arrow_kw = dict(transform=ax.transAxes, lw=1.5,
                arrowstyle='->', color='#555')
for y_start, y_end in [(0.86, 0.80), (0.71, 0.65),
                        (0.55, 0.49), (0.38, 0.33), (0.22, 0.17)]:
    ax.annotate('', xy=(0.5, y_end), xytext=(0.5, y_start),
                xycoords='axes fraction', textcoords='axes fraction',
                arrowprops=dict(arrowstyle='->', color='#555', lw=1.8))

# Rückpfeil
ax.annotate('', xy=(0.85, 0.60), xytext=(0.85, 0.28),
            xycoords='axes fraction', textcoords='axes fraction',
            arrowprops=dict(arrowstyle='->', color='#d6604d',
                            lw=2, connectionstyle='arc3,rad=0.0'))
ax.text(0.93, 0.44, 'nicht\nkonvergiert', ha='center', va='center',
        fontsize=9, color='#d6604d', transform=ax.transAxes)

ax.annotate('', xy=(0.5, 0.28+0.05), xytext=(0.85, 0.28+0.05),
            xycoords='axes fraction', textcoords='axes fraction',
            arrowprops=dict(arrowstyle='->', color='#d6604d', lw=2))

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

Abbildung 5.1: SCF-Iterationsschema. Die Schleife läuft, bis die Dichtematrix zwischen zwei Iterationen konvergiert.

---

5.5 HF-Rechnung für H₂O mit PySCF

5.5.1 Aufsetzen und Ausführen

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  = 'sto-3g'
mol.charge = 0
mol.spin   = 0
mol.verbose = 0
mol.build()

# Restricted Hartree-Fock (RHF)
mf = scf.RHF(mol)
mf.conv_tol = 1e-10    # Konvergenzkriterium für Energie
E_hf = mf.kernel()

print(f"HF-Gesamtenergie (STO-3G):  {E_hf:.8f} Hartree")
print(f"                          = {E_hf * 627.509:.4f} kcal/mol")
print(f"\nAnzahl SCF-Iterationen:     {mf.scf_summary.get('niter', '?')}")
print(f"Konvergiert:                {mf.converged}")

HF-Gesamtenergie (STO-3G):  -74.96292783 Hartree
                          = -47039.9119 kcal/mol

Anzahl SCF-Iterationen:     ?
Konvergiert:                True

# Orbitalenergien (in Hartree)
eps   = mf.mo_energy
occ   = mf.mo_occ
n_occ = int(np.sum(occ > 0))

print("Orbitalenergien (Hartree)  |  Besetzung  |  Typ")
print("-" * 52)
for i, (e, o) in enumerate(zip(eps, occ)):
    typ = "besetzt (HOMO)" if i == n_occ-1 else \
          "virtuell (LUMO)" if i == n_occ else \
          "besetzt" if o > 0 else "virtuell"
    print(f"  MO {i+1:2d}:  {e:10.5f} Hartree   |   {int(o)}       |  {typ}")

Orbitalenergien (Hartree)  |  Besetzung  |  Typ
----------------------------------------------------
  MO  1:   -20.24174 Hartree   |   2       |  besetzt
  MO  2:    -1.26841 Hartree   |   2       |  besetzt
  MO  3:    -0.61794 Hartree   |   2       |  besetzt
  MO  4:    -0.45299 Hartree   |   2       |  besetzt
  MO  5:    -0.39124 Hartree   |   2       |  besetzt (HOMO)
  MO  6:     0.60568 Hartree   |   0       |  virtuell (LUMO)
  MO  7:     0.74240 Hartree   |   0       |  virtuell

fig, ax = plt.subplots(figsize=(5, 7))

colors_mo = {True: '#2166ac', False: '#aaaaaa'}
lw, llen = 2.5, 0.35

for i, (e, o) in enumerate(zip(eps, occ)):
    occupied = o > 0
    col = colors_mo[occupied]
    ax.hlines(e, 0.5-llen/2, 0.5+llen/2, color=col, lw=lw)
    if occupied:
        for dx in [-0.06, 0.06]:
            dy = 0.03
            ax.annotate('', xy=(0.5+dx, e+dy), xytext=(0.5+dx, e-dy),
                        arrowprops=dict(arrowstyle='->' if dx < 0 else '<-',
                                        color=col, lw=1.5))
    lbl = ''
    if i == n_occ - 1: lbl = ' HOMO'
    if i == n_occ:     lbl = ' LUMO'
    ax.text(0.5+llen/2+0.03, e, f'MO {i+1}{lbl}',
            va='center', fontsize=9,
            color='#333' if occupied else '#999')

ax.axhline(0, color='gray', lw=0.8, ls=':', alpha=0.5)
ax.set_xlim(0.0, 1.3)
ax.set_ylabel('Orbitalenergie / Hartree', fontsize=12)
ax.set_title('H₂O Orbitalenergien\n(RHF / STO-3G)', fontsize=12)
ax.set_xticks([])
ax.grid(True, axis='y', alpha=0.25)

from matplotlib.patches import Patch
ax.legend(handles=[Patch(color='#2166ac', label='besetzt'),
                   Patch(color='#aaa',    label='virtuell')],
          fontsize=10, loc='upper right')
plt.tight_layout()
plt.show()

Abbildung 5.2: Orbitalenergien von H₂O (RHF/STO-3G). Besetzte Orbitale in Blau, virtuelle in Grau. HOMO und LUMO sind markiert.

5.5.2 Die Dichte- und Fock-Matrix ansehen

# Dichtematrix P = 2 * C_occ @ C_occ.T  (nur besetzte MOs)
C     = mf.mo_coeff          # MO-Koeffizienten (N_b × N_b)
C_occ = C[:, :n_occ]         # nur besetzte Spalten
P     = 2 * C_occ @ C_occ.T  # Dichtematrix

# Fock-Matrix
F = mf.get_fock()

print("Dichtematrix P (7×7, STO-3G):")
print(np.round(P, 4))
print(f"\nSpur von P = {np.trace(P):.4f}  (= Anzahl Elektronen = {mol.nelectron})")

Dichtematrix P (7×7, STO-3G):
[[ 2.1063 -0.446  -0.      0.     -0.1086 -0.0284 -0.0284]
 [-0.446   1.9673  0.     -0.      0.6177 -0.034  -0.034 ]
 [-0.      0.      2.      0.      0.      0.      0.    ]
 [ 0.     -0.      0.      0.7355 -0.      0.5397 -0.5397]
 [-0.1086  0.6177  0.     -0.      1.2403 -0.4727 -0.4727]
 [-0.0284 -0.034   0.      0.5397 -0.4727  0.6009 -0.1912]
 [-0.0284 -0.034   0.     -0.5397 -0.4727 -0.1912  0.6009]]

Spur von P = 9.2511  (= Anzahl Elektronen = 10)

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

for ax, mat, title, cmap in zip(
    axes, [P, F],
    ['Dichtematrix $\\mathbf{P}$', 'Fock-Matrix $\\mathbf{F}$'],
    ['Blues', 'RdBu_r']
):
    im = ax.imshow(mat, cmap=cmap, aspect='auto')
    plt.colorbar(im, ax=ax, fraction=0.046)
    ax.set_title(title, fontsize=12)
    ax.set_xlabel('Basisfunktions-Index $\\mu$', fontsize=10)
    ax.set_ylabel('Basisfunktions-Index $\\nu$', fontsize=10)
    nb = mat.shape[0]
    ax.set_xticks(range(nb))
    ax.set_yticks(range(nb))

# Basisfunktionsbeschriftungen
# ao_labels(fmt=None) → Liste von Tupeln (atom_idx, atom_sym, nl, m)
bf_labels = [f'{ao[1]}\n{ao[2]}' for ao in mol.ao_labels(fmt=None)]
for ax in axes:
    ax.set_xticklabels(bf_labels, fontsize=7)
    ax.set_yticklabels(bf_labels, fontsize=7)

plt.tight_layout()
plt.show()

Abbildung 5.3: Dichtematrix P (links) und Fock-Matrix F (rechts) für H₂O (STO-3G). Die Dichtematrix kodiert, wie stark jede Basisfunktion zur Gesamtelektronendichte beiträgt.

5.5.3 Konvergenzverlauf beobachten

# SCF mit gespeichertem Verlauf
mol2 = mol.copy()
mol2.verbose = 0
mf2 = scf.RHF(mol2)

E_history = []
def callback(envs):
    E_history.append(envs['e_tot'])
mf2.callback = callback
mf2.kernel()

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

# Absolute Energie
axes[0].plot(range(1, len(E_history)+1), E_history,
             'o-', color='#2166ac', lw=2, ms=7)
axes[0].set_xlabel('SCF-Iteration', fontsize=12)
axes[0].set_ylabel('Energie / Hartree', fontsize=12)
axes[0].set_title('Absolute Energie pro Iteration', fontsize=12)
axes[0].grid(True, alpha=0.3)

# Energieänderung (log)
dE = np.abs(np.diff(E_history))
if len(dE) > 0:
    axes[1].semilogy(range(2, len(E_history)+1), dE,
                     's-', color='#d6604d', lw=2, ms=7)
    axes[1].axhline(1e-10, color='gray', ls='--', lw=1.2,
                    label='Konvergenzkriterium $10^{-10}$')
    axes[1].set_xlabel('SCF-Iteration', fontsize=12)
    axes[1].set_ylabel('$|\\Delta E|$ / Hartree', fontsize=12)
    axes[1].set_title('Energieänderung (log)', fontsize=12)
    axes[1].legend(fontsize=10)
    axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Abbildung 5.4: Konvergenz der SCF-Energie für H₂O (RHF/STO-3G). Die Energie fällt monoton und konvergiert nach wenigen Iterationen.

---

5.6 Visualisierung der Molekülorbitale

fig, axes = plt.subplots(1, n_occ, figsize=(12, 4), sharey=True)

bf_lbls = [f'{ao[1]}-{ao[2]}' for ao in mol.ao_labels(fmt=None)]
mo_names = ['MO 1\n(O 1s)', 'MO 2\n(O 2s)', 'MO 3\n(nbO)', 'MO 4\n(σ OH)', 'MO 5\n(HOMO)']

for idx, ax in enumerate(axes):
    coeffs = C[:, idx]
    colors = ['#2166ac' if c >= 0 else '#d6604d' for c in coeffs]
    bars = ax.barh(range(len(coeffs)), coeffs, color=colors, edgecolor='white')
    ax.axvline(0, color='black', lw=0.8)
    ax.set_title(mo_names[idx], fontsize=10)
    ax.set_xlim(-1.1, 1.1)
    if idx == 0:
        ax.set_yticks(range(len(bf_lbls)))
        ax.set_yticklabels(bf_lbls, fontsize=9)
    ax.set_xlabel('$c_{\\mu m}$', fontsize=9)
    ax.grid(True, axis='x', alpha=0.3)

fig.suptitle('MO-Koeffizienten $c_{\\mu m}$: besetzte Orbitale von H₂O (STO-3G)',
             fontsize=12)
plt.tight_layout()
plt.show()

Abbildung 5.5: Orbitalkoeffizienten der fünf besetzten MOs von H₂O (STO-3G). Jede Spalte zeigt, wie stark jede Basisfunktion zum jeweiligen MO beiträgt.

HinweisWas die Koeffizienten bedeuten

MO 1 ist fast ausschliesslich aus der O-1s-Basisfunktion aufgebaut – es ist ein “core”-Orbital, das tief im Inneren des Sauerstoffatoms liegt und kaum an der Bindung beteiligt ist. HOMO (MO 5) hingegen hat grosse Beiträge von den O-2p-Orbitalen – es entspricht den freien Elektronenpaaren des Sauerstoffs.

---

5.7 Grenzen von Hartree-Fock: Elektronenkorrelation

Die fundamentale Näherung von HF ist die einzige Slater-Determinante. Damit bewegt sich jedes Elektron in einem gemittelten Feld – die dynamische Korrelation, dass Elektronen sich instantan ausweichen, fehlt.

Die Korrelationsenergie ist definiert als:

\[E_\text{corr} = E_\text{exakt} - E_\text{HF} \tag{3.12}\]

Sie ist stets negativ (HF überschätzt die Energie) und typischerweise \(-0.5\) bis \(-1.5\) eV pro Elektronenpaar.

# Vergleich verschiedener Basissätze und Methoden für H2O
from pyscf import mp, cc

results = {}
for basis in ['sto-3g', '6-31g', '6-31g*']:
    m = gto.Mole(atom=mol.atom, basis=basis, verbose=0)
    m.build()
    mf_tmp = scf.RHF(m).run()
    results[f'HF/{basis}'] = mf_tmp.e_tot

# MP2 und CCSD für 6-31g*
m631gs = gto.Mole(atom=mol.atom, basis='6-31g*', verbose=0)
m631gs.build()
mf631gs = scf.RHF(m631gs).run()

mp2 = mp.MP2(mf631gs).run()
results['MP2/6-31g*']  = mf631gs.e_tot + mp2.e_corr

ccsd = cc.CCSD(mf631gs).run()
results['CCSD/6-31g*'] = mf631gs.e_tot + ccsd.e_corr

print(f"{'Methode':<18}  {'E / Hartree':>14}  {'E / kcal/mol':>14}")
print("-" * 52)
for k, v in results.items():
    print(f"  {k:<16}  {v:>14.6f}  {v*627.509:>14.2f}")

print(f"\n  Korrelationsenergie (MP2-HF, 6-31g*):  "
      f"{(results['MP2/6-31g*'] - results['HF/6-31g*'])*627.509:.2f} kcal/mol")
print(f"  Korrelationsenergie (CCSD-HF, 6-31g*): "
      f"{(results['CCSD/6-31g*'] - results['HF/6-31g*'])*627.509:.2f} kcal/mol")

Methode                E / Hartree    E / kcal/mol
----------------------------------------------------
  HF/sto-3g             -74.962928       -47039.91
  HF/6-31g              -75.983998       -47680.64
  HF/6-31g*             -76.009132       -47696.41
  MP2/6-31g*            -76.195313       -47813.24
  CCSD/6-31g*           -76.204317       -47818.89

  Korrelationsenergie (MP2-HF, 6-31g*):  -116.83 kcal/mol
  Korrelationsenergie (CCSD-HF, 6-31g*): -122.48 kcal/mol

ref = results['HF/sto-3g']
methods = list(results.keys())
energies_rel = [(results[k] - ref) * 627.509 for k in methods]

colors_bar = ['#4393c3']*3 + ['#f4a582', '#d6604d']
fig, ax = plt.subplots(figsize=(8, 4))
bars = ax.bar(methods, energies_rel, color=colors_bar, edgecolor='white', lw=1.2)
ax.axhline(0, color='gray', lw=0.8, ls='--')
ax.set_ylabel('$\\Delta E$ relativ zu HF/STO-3G / kcal/mol', fontsize=11)
ax.set_title('Relative Energien: H₂O mit verschiedenen Methoden', fontsize=12)
for bar, val in zip(bars, energies_rel):
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() - 1.5,
            f'{val:.1f}', ha='center', va='top', fontsize=9, color='white',
            fontweight='bold')
ax.tick_params(axis='x', rotation=15)
ax.grid(True, axis='y', alpha=0.3)
plt.tight_layout()
plt.show()

Abbildung 5.6: Energien von H₂O mit verschiedenen Methoden und Basissätzen (relativ zur HF/STO-3G-Energie). Mit grösserem Basissatz und Korrelationsmethoden sinkt die Energie.

---

5.8 Zusammenfassung

\[\boxed{ \underbrace{\mathbf{F}(\mathbf{P})}_{\text{hängt von } \mathbf{c} \text{ ab}} \mathbf{c} = \mathbf{S}\,\mathbf{c}\,\boldsymbol{\varepsilon} \quad \Longrightarrow \quad \text{SCF-Iteration bis Konvergenz} }\]

Schritt	Mathematik	Implementierung
LCAO-Ansatz	\(\psi_m = \sum_\mu c_{\mu m}\chi_\mu\)	Basissatzwahl (Kap. 04)
Fock-Matrix	\(F_{\nu\mu} = h_{\nu\mu} + \sum P_{\lambda\sigma}(\ldots)\)	\(O(N_b^4)\) Integrale
Roothaan-Gl.	\(\mathbf{F}\mathbf{c} = \mathbf{S}\mathbf{c}\boldsymbol{\varepsilon}\)	verallg. Eigenwertproblem
SCF-Schleife	Iteration bis \(\|\Delta\mathbf{P}\| < \epsilon\)	~10–50 Iterationen
Gesamtenergie	\(E = 2\sum\varepsilon_m - \sum(2J-K)\)	Korrekte Doppelzählung

TippVorschau

Kapitel 04 befasst sich mit der entscheidenden praktischen Frage: Welche Basisfunktionen \(\chi_\mu\) soll man wählen? Die Antwort – Gaussfunktionen in verschiedenen Kontraktionen – beeinflusst massgeblich Genauigkeit und Rechenaufwand jeder HF- oder DFT-Rechnung.