import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
plt.rcParams.update({'figure.dpi': 120, 'font.size': 11})
try:
from pyscf import gto, scf, mp, ao2mo
PYSCF_AVAILABLE = True
except ImportError:
PYSCF_AVAILABLE = False
print("PySCF nicht gefunden. Installiere mit: uv add pyscf")9 Ü2 · Hartree-Fock
SCF-Konvergenz, Roothaan-Gleichungen und Elektronenkorrelation
9.1 Lernziele
Nach dieser Übung kannst du:
- SCF-Konvergenzverhalten analysieren und visualisieren
- Überlappmatrix \(S\) und Fock-Matrix \(F\) direkt inspizieren
- Den Einfluss des Basissatzes auf die HF-Energie systematisch untersuchen
- Elektronenkorrelation über den Vergleich HF vs. MP2 quantifizieren
- Das Koopmans-Theorem zur Vorhersage von Ionisierungspotentialen anwenden
9.2 Aufgabe 1 · SCF-Konvergenzverlauf
Das SCF-Verfahren iteriert bis zur Selbstkonsistenz. Beobachte den Verlauf der Gesamtenergie über die SCF-Iterationen für Wasser mit verschiedenen Basissätzen.
if not PYSCF_AVAILABLE:
raise SystemExit("PySCF benötigt")
def h2o_mol(basis='sto-3g'):
"""Erzeugt ein H2O-Molekülobjekt mit der gewünschten Basis."""
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 = basis
mol.verbose = 0
mol.build()
return mol
def run_scf_tracked(mol, max_cycle=50):
"""RHF mit Energieprotokoll pro Iteration."""
mf = scf.RHF(mol)
mf.max_cycle = max_cycle
mf.conv_tol = 1e-10
energies = []
def callback(envs):
energies.append(envs['e_tot'])
mf.callback = callback
mf.kernel()
return mf, np.array(energies)
# Vergleich verschiedener Basissätze
bases = ['sto-3g', '6-31g', '6-31g*', 'cc-pvdz']
results = {}
for basis in bases:
mol = h2o_mol(basis)
mf, energies = run_scf_tracked(mol)
results[basis] = {'mf': mf, 'energies': energies, 'E': mf.e_tot}
print(f"{basis:12s}: E = {mf.e_tot:.6f} Hartree ({len(energies)} Iterationen)")sto-3g : E = -74.962947 Hartree (7 Iterationen)
6-31g : E = -75.983993 Hartree (9 Iterationen)
6-31g* : E = -76.009128 Hartree (8 Iterationen)
cc-pvdz : E = -76.026794 Hartree (9 Iterationen)fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
for basis, res in results.items():
E_arr = res['energies']
ax1.plot(range(1, len(E_arr)+1), E_arr, marker='o', ms=4, label=basis)
# Konvergenz: Betrag der Energie-Änderung
if len(E_arr) > 1:
delta = np.abs(np.diff(E_arr))
ax2.semilogy(range(1, len(delta)+1), delta, marker='s', ms=4, label=basis)
ax1.set_xlabel('SCF-Iteration')
ax1.set_ylabel('Gesamtenergie / Hartree')
ax1.set_title('SCF-Energieverlauf')
ax1.legend()
ax2.set_xlabel('SCF-Iteration')
ax2.set_ylabel('|ΔE| / Hartree')
ax2.set_title('Konvergenz (log-Skala)')
ax2.axhline(1e-10, ls='--', color='gray', label='Konvergenzkriterium')
ax2.legend(fontsize=9)
plt.suptitle('H₂O RHF – SCF-Konvergenz')
plt.tight_layout()
plt.show()
9.3 Aufgabe 2 · Überlappmatrix und Fock-Matrix
Inspiziere die Matrizen des SCF-Verfahrens direkt.
mol_sto3g = h2o_mol('sto-3g')
mf_sto3g, _ = run_scf_tracked(mol_sto3g)
# Matrizen aus dem konvergierten SCF
S = mol_sto3g.intor('int1e_ovlp') # Überlappmatrix
F = mf_sto3g.get_fock() # Fock-Matrix (konvergiert)
P = mf_sto3g.make_rdm1() # Dichtematrix
ao_labels = [f"{ao[1]} {ao[2]}{ao[3]}" for ao in mol_sto3g.ao_labels(fmt=None)]
fig, axes = plt.subplots(1, 3, figsize=(14, 4))
for ax, mat, title in zip(axes, [S, F, P], ['Überlappmatrix S', 'Fock-Matrix F', 'Dichtematrix P']):
im = ax.imshow(mat, cmap='RdBu_r', aspect='auto')
plt.colorbar(im, ax=ax, shrink=0.8)
ax.set_title(title)
ax.set_xticks(range(len(ao_labels)))
ax.set_yticks(range(len(ao_labels)))
ax.set_xticklabels(ao_labels, rotation=45, ha='right', fontsize=8)
ax.set_yticklabels(ao_labels, fontsize=8)
plt.suptitle('H₂O/STO-3G – Matrizen des SCF-Verfahrens')
plt.tight_layout()
plt.show()
print("\nBasisorbitale (AO-Labels):")
for i, l in enumerate(ao_labels):
print(f" {i}: {l}")
Basisorbitale (AO-Labels):
0: O 1s
1: O 2s
2: O 2px
3: O 2py
4: O 2pz
5: H 1s
6: H 1s9.4 Aufgabe 3 · Basissatz-Konvergenzstudie
Untersuche systematisch, wie die HF-Energie mit der Basissatzgrösse konvergiert.
basis_series = [
('STO-3G', 'sto-3g', 7),
('3-21G', '3-21g', 13),
('6-31G', '6-31g', 13),
('6-31G*', '6-31g*', 19),
('6-31G**', '6-31g**', 23),
('cc-pVDZ', 'cc-pvdz', 24),
('cc-pVTZ', 'cc-pvtz', 58),
]
# HF-Limit für H2O (Referenz aus Literatur): ~-76.0675 Hartree
E_hf_limit = -76.0675
data = []
for name, basis, n_ao in basis_series:
mol = h2o_mol(basis)
mf = scf.RHF(mol)
mf.verbose = 0
E = mf.kernel()
data.append({'Basissatz': name, 'N_AO': n_ao, 'E_HF [Hartree]': E,
'ΔE_HF [mHartree]': (E - E_hf_limit)*1000})
print(f"{name:12s}: {E:.6f} Hartree (N_AO={n_ao})")
df_basis = pd.DataFrame(data)
print("\n", df_basis.to_string(index=False))STO-3G : -74.962947 Hartree (N_AO=7)
3-21G : -75.585395 Hartree (N_AO=13)
6-31G : -75.983993 Hartree (N_AO=13)
6-31G* : -76.009128 Hartree (N_AO=19)
6-31G** : -76.022641 Hartree (N_AO=23)
cc-pVDZ : -76.026794 Hartree (N_AO=24)
cc-pVTZ : -76.057161 Hartree (N_AO=58)
Basissatz N_AO E_HF [Hartree] ΔE_HF [mHartree]
STO-3G 7 -74.962947 1104.553343
3-21G 13 -75.585395 482.104758
6-31G 13 -75.983993 83.506772
6-31G* 19 -76.009128 58.372212
6-31G** 23 -76.022641 44.859063
cc-pVDZ 24 -76.026794 40.706355
cc-pVTZ 58 -76.057161 10.339319fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
ax1.plot(df_basis['N_AO'], df_basis['E_HF [Hartree]'], 'o-', color='steelblue', ms=7)
ax1.axhline(E_hf_limit, ls='--', color='red', label=f'HF-Limit ≈ {E_hf_limit:.4f} Eh')
for _, row in df_basis.iterrows():
ax1.annotate(row['Basissatz'], (row['N_AO'], row['E_HF [Hartree]']),
textcoords='offset points', xytext=(5, 3), fontsize=8)
ax1.set_xlabel('Anzahl Basisorbitale')
ax1.set_ylabel('E_HF / Hartree')
ax1.set_title('HF-Energie vs. Basissatzgrösse')
ax1.legend()
ax2.semilogy(df_basis['N_AO'], np.abs(df_basis['ΔE_HF [mHartree]']), 'o-', color='darkorange', ms=7)
ax2.set_xlabel('Anzahl Basisorbitale')
ax2.set_ylabel('|E - E_HF_limit| / mHartree')
ax2.set_title('Abweichung vom HF-Limit (log)')
plt.suptitle('H₂O – Basissatz-Konvergenzstudie')
plt.tight_layout()
plt.show()
9.5 Aufgabe 4 · Koopmans-Theorem und Ionisierungspotentiale
Das Koopmans-Theorem besagt: Das negative der Orbitalenergie des HOMO entspricht dem (ersten) Ionisierungspotential: \[\text{IP} \approx -\varepsilon_{\text{HOMO}}\]
Vergleiche die Koopmans-Vorhersage mit dem Ergebnis aus einer direkten \(\Delta\)SCF-Rechnung.
molecules = {
'H₂O': """O 0 0 0; H 0.757 0.586 0; H -0.757 0.586 0""",
'NH₃': """N 0 0 0.116; H 0 0.939 -0.269; H 0.813 -0.469 -0.269; H -0.813 -0.469 -0.269""",
'CH₄': """C 0 0 0; H 0.629 0.629 0.629; H -0.629 -0.629 0.629; H -0.629 0.629 -0.629; H 0.629 -0.629 -0.629""",
}
# Experimentelle erste Ionisierungspotentiale in eV (NIST)
IP_exp = {'H₂O': 12.62, 'NH₃': 10.82, 'CH₄': 12.61}
results_ip = []
for name, atom_str in molecules.items():
# Neutrales Molekül
mol = gto.Mole(atom=atom_str, basis='6-31g*', verbose=0)
mol.build()
mf = scf.RHF(mol)
E_neutral = mf.kernel()
# Koopmans-IP: -ε_HOMO
homo_idx = int(mol.nelectron // 2) - 1
IP_koopmans = -mf.mo_energy[homo_idx] * 27.211 # Hartree → eV
# ΔSCF: Kation berechnen
mol_cat = gto.Mole(atom=atom_str, basis='6-31g*', charge=1, spin=1, verbose=0)
mol_cat.build()
mf_cat = scf.UHF(mol_cat)
E_cation = mf_cat.kernel()
IP_delta = (E_cation - E_neutral) * 27.211 # Hartree → eV
results_ip.append({
'Molekül': name,
'IP Koopmans [eV]': IP_koopmans,
'IP ΔSCF [eV]': IP_delta,
'IP exp. [eV]': IP_exp[name],
'Fehler Koopmans [eV]': IP_koopmans - IP_exp[name],
'Fehler ΔSCF [eV]': IP_delta - IP_exp[name],
})
print(f"{name}: Koopmans={IP_koopmans:.2f} eV | ΔSCF={IP_delta:.2f} eV | exp={IP_exp[name]:.2f} eV")
df_ip = pd.DataFrame(results_ip)
print()
print(df_ip.to_string(index=False))H₂O: Koopmans=13.54 eV | ΔSCF=10.85 eV | exp=12.62 eV
NH₃: Koopmans=11.47 eV | ΔSCF=9.29 eV | exp=10.82 eV
CH₄: Koopmans=14.82 eV | ΔSCF=13.38 eV | exp=12.61 eV
Molekül IP Koopmans [eV] IP ΔSCF [eV] IP exp. [eV] Fehler Koopmans [eV] Fehler ΔSCF [eV]
H₂O 13.539739 10.848314 12.62 0.919739 -1.771686
NH₃ 11.473014 9.293784 10.82 0.653014 -1.526216
CH₄ 14.823016 13.377483 12.61 2.213016 0.7674839.6 Aufgabe 5 · Elektronenkorrelation: HF vs. MP2
Berechne die Korrelationsenergie \(E_\text{corr} = E_\text{MP2} - E_\text{HF}\) für eine Reihe kleiner Moleküle.
Die MP2-Methode (Møller-Plesset 2. Ordnung) ist die einfachste post-HF-Methode und erfasst dynamische Elektronenkorrelation.
from pyscf import mp
molecules_corr = {
'H₂': "H 0 0 0; H 0 0 0.741",
'H₂O': "O 0 0 0; H 0.757 0.586 0; H -0.757 0.586 0",
'NH₃': "N 0 0 0.116; H 0 0.939 -0.269; H 0.813 -0.469 -0.269; H -0.813 -0.469 -0.269",
'HF': "F 0 0 0; H 0 0 0.917",
}
corr_data = []
for name, atom_str in molecules_corr.items():
mol = gto.Mole(atom=atom_str, basis='cc-pvdz', verbose=0)
mol.build()
# HF
mf = scf.RHF(mol)
E_hf = mf.kernel()
# MP2
mp2 = mp.MP2(mf)
E_corr, _ = mp2.kernel()
E_mp2 = E_hf + E_corr
corr_data.append({
'Molekül': name,
'E_HF [Eh]': E_hf,
'E_MP2 [Eh]': E_mp2,
'E_corr [mEh]': E_corr * 1000,
'E_corr [kJ/mol]': E_corr * 2625.5,
})
print(f"{name:6s}: E_HF={E_hf:.4f}, E_corr={E_corr*1000:.1f} mEh ({E_corr*2625.5:.1f} kJ/mol)")
df_corr = pd.DataFrame(corr_data)
# Visualisierung
fig, ax = plt.subplots(figsize=(8, 5))
x = np.arange(len(df_corr))
ax.bar(x, np.abs(df_corr['E_corr [kJ/mol]']), color='steelblue', alpha=0.8)
ax.set_xticks(x)
ax.set_xticklabels(df_corr['Molekül'], fontsize=12)
ax.set_ylabel('|E_corr| / kJ mol⁻¹')
ax.set_title('Korrelationsenergie (MP2/cc-pVDZ)')
for i, v in enumerate(np.abs(df_corr['E_corr [kJ/mol]'])):
ax.text(i, v + 5, f'{v:.0f}', ha='center', fontsize=10)
plt.tight_layout()
plt.show()H₂ : E_HF=-1.1287, E_corr=-26.4 mEh (-69.3 kJ/mol)
H₂O : E_HF=-76.0268, E_corr=-204.0 mEh (-535.5 kJ/mol)
NH₃ : E_HF=-56.1956, E_corr=-189.1 mEh (-496.6 kJ/mol)
HF : E_HF=-100.0194, E_corr=-203.8 mEh (-535.0 kJ/mol)
9.6.1 Reflexionsfragen
- Welches Molekül zeigt die grösste absolute Korrelationsenergie? Warum?
- Warum ist die Korrelationsenergie immer negativ?
- Das Koopmans-Theorem überschätzt typischerweise das IP. Warum systematisch?
- Welche Näherung wird im HF-Verfahren durch die Slater-Determinante gemacht, und was folgt daraus für die Elektronenkorrelation?
9.7 Zusammenfassung
In dieser Übung hast du:
- SCF-Konvergenzkurven analysiert und den Einfluss des Basissatzes quantifiziert
- Überlapp-, Fock- und Dichtematrix direkt inspiziert
- Ionisierungspotentiale via Koopmans-Theorem und ΔSCF berechnet
- Elektronenkorrelationsenergien auf MP2-Niveau berechnet und verglichen