10  Ü3 · Dichtefunktionaltheorie

LDA, GGA, Hybrid-Funktionale und Elektronendichte

Autor:in

ZHAW Institut für Data Science

10.1 Lernziele

Nach dieser Übung kannst du:

  • DFT-Rechnungen mit LDA, GGA (PBE) und Hybrid-Funktionalen (B3LYP) durchführen
  • Energien und HOMO/LUMO-Gaps mit HF und unterschiedlichen Funktionalen vergleichen
  • Die Elektronendichte auf einem 3D-Gitter berechnen und schneiden
  • Geometrieoptimierungen mit DFT durchführen
  • Den Einfluss der Selbstwechselwirkungskorrektur (SIC) diskutieren
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, dft
    PYSCF_AVAILABLE = True
except ImportError:
    PYSCF_AVAILABLE = False
    print("PySCF nicht gefunden. Installiere mit: uv add pyscf")

10.2 Aufgabe 1 · Vergleich von Funktionalen für H₂O

Berechne die Gesamtenergie und den HOMO/LUMO-Gap von H₂O mit HF und verschiedenen DFT-Funktionalen.

if not PYSCF_AVAILABLE:
    raise SystemExit("PySCF benötigt")

mol_h2o = gto.Mole(
    atom="O 0 0 0; H 0.757 0.586 0; H -0.757 0.586 0",
    basis='6-31g*',
    verbose=0
)
mol_h2o.build()

n_electrons = mol_h2o.nelectron
homo_idx = n_electrons // 2 - 1  # (geschlossenschalig, 0-indiziert)

methods = {
    'HF':    lambda mol: scf.RHF(mol),
    'LDA':   lambda mol: dft.RKS(mol, xc='lda,vwn'),
    'PBE':   lambda mol: dft.RKS(mol, xc='pbe,pbe'),
    'BLYP':  lambda mol: dft.RKS(mol, xc='b88,lyp'),
    'B3LYP': lambda mol: dft.RKS(mol, xc='b3lyp'),
    'PBE0':  lambda mol: dft.RKS(mol, xc='pbe0'),
    'M06-2X': lambda mol: dft.RKS(mol, xc='m062x'),
}

data = []
for name, builder in methods.items():
    mf = builder(mol_h2o)
    mf.verbose = 0
    E = mf.kernel()
    
    eps = mf.mo_energy
    homo = eps[homo_idx] * 27.211
    lumo = eps[homo_idx + 1] * 27.211
    gap  = lumo - homo
    
    data.append({'Methode': name, 'E [Hartree]': E,
                 'HOMO [eV]': homo, 'LUMO [eV]': lumo, 'Gap [eV]': gap})
    print(f"{name:8s}: E={E:.4f} Eh  HOMO={homo:.2f} eV  LUMO={lumo:.2f} eV  Gap={gap:.2f} eV")

df = pd.DataFrame(data)
print()
print(df.to_string(index=False))
HF      : E=-76.0091 Eh  HOMO=-13.54 eV  LUMO=5.84 eV  Gap=19.38 eV
LDA     : E=-75.8409 Eh  HOMO=-6.26 eV  LUMO=1.23 eV  Gap=7.50 eV
PBE     : E=-76.3198 Eh  HOMO=-6.18 eV  LUMO=1.29 eV  Gap=7.48 eV
BLYP    : E=-76.3855 Eh  HOMO=-6.07 eV  LUMO=1.16 eV  Gap=7.23 eV
B3LYP   : E=-76.4068 Eh  HOMO=-7.91 eV  LUMO=1.87 eV  Gap=9.78 eV
PBE0    : E=-76.3238 Eh  HOMO=-8.28 eV  LUMO=2.36 eV  Gap=10.64 eV
M06-2X  : E=-76.3721 Eh  HOMO=-10.22 eV  LUMO=3.28 eV  Gap=13.50 eV

Methode  E [Hartree]  HOMO [eV]  LUMO [eV]  Gap [eV]
     HF   -76.009128 -13.539739   5.841790 19.381528
    LDA   -75.840921  -6.264228   1.232827  7.497055
    PBE   -76.319775  -6.182715   1.294792  7.477507
   BLYP   -76.385493  -6.069906   1.161357  7.231262
  B3LYP   -76.406789  -7.913954   1.865413  9.779367
   PBE0   -76.323827  -8.281060   2.358992 10.640052
 M06-2X   -76.372096 -10.222013   3.277934 13.499947
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))

colors = plt.cm.Set2(np.linspace(0, 1, len(df)))

# HOMO/LUMO-Energien
x = np.arange(len(df))
ax1.bar(x - 0.2, df['HOMO [eV]'], 0.35, label='HOMO', color='steelblue', alpha=0.8)
ax1.bar(x + 0.2, df['LUMO [eV]'], 0.35, label='LUMO', color='darkorange', alpha=0.8)
ax1.axhline(-12.62, ls='--', color='blue', alpha=0.5, label='IP exp. (H₂O)')
ax1.set_xticks(x)
ax1.set_xticklabels(df['Methode'], rotation=30, ha='right')
ax1.set_ylabel('Orbitalenergie / eV')
ax1.set_title('HOMO/LUMO-Energien')
ax1.legend()

# HOMO-LUMO-Gap
ax2.bar(x, df['Gap [eV]'], color=[c for c in colors], alpha=0.85)
ax2.set_xticks(x)
ax2.set_xticklabels(df['Methode'], rotation=30, ha='right')
ax2.set_ylabel('HOMO-LUMO-Gap / eV')
ax2.set_title('HOMO-LUMO-Gap')
for i, v in enumerate(df['Gap [eV]']):
    ax2.text(i, v + 0.1, f'{v:.2f}', ha='center', fontsize=9)

plt.suptitle('H₂O/6-31G* – Vergleich von HF und DFT-Funktionalen')
plt.tight_layout()
plt.show()

10.3 Aufgabe 2 · Elektronendichte visualisieren

Berechne die Elektronendichte \(\rho(\mathbf{r})\) auf einem 2D-Schnitt durch das Molekül.

from pyscf import dft as pydft

mol = mol_h2o
mf_b3lyp = dft.RKS(mol, xc='b3lyp')
mf_b3lyp.verbose = 0
mf_b3lyp.kernel()

# Dichtematrix im AO-Basis
dm = mf_b3lyp.make_rdm1()

# 2D-Gitter in der Molekülebene (xz-Ebene, y=0)
N_grid = 100
x_vals = np.linspace(-4, 4, N_grid)   # in Bohr
z_vals = np.linspace(-3, 5, N_grid)

# Gitterpunkte als (N², 3)-Array
xx, zz = np.meshgrid(x_vals, z_vals)
grid_coords = np.column_stack([xx.ravel(), np.zeros(N_grid**2), zz.ravel()])

# AO-Werte auf dem Gitter berechnen
ao_vals = mol.eval_gto('GTOval_sph', grid_coords)  # (N_points, N_ao)

# Elektronendichte: rho = sum_ij P_ij * phi_i * phi_j
rho = np.einsum('pi,ij,pj->p', ao_vals, dm, ao_vals)
rho = rho.reshape(N_grid, N_grid)

# Kernpositionen
atom_coords_bohr = mol.atom_coords()  # in Bohr

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))

# Lineare Skala
cf1 = ax1.contourf(xx, zz, rho, levels=30, cmap='Blues')
cs1 = ax1.contour(xx, zz, rho, levels=15, colors='white', alpha=0.4, linewidths=0.5)
plt.colorbar(cf1, ax=ax1, label='ρ(r) / e Bohr⁻³')
for i, (sym, coord) in enumerate([(mol.atom_symbol(i), atom_coords_bohr[i]) for i in range(mol.natm)]):
    ax1.plot(coord[0], coord[2], 'ro', ms=8)
    ax1.text(coord[0]+0.1, coord[2]+0.1, sym, color='red', fontsize=10, fontweight='bold')
ax1.set_xlabel('x / Bohr')
ax1.set_ylabel('z / Bohr')
ax1.set_title('Elektronendichte ρ(r) – lineare Skala')
ax1.set_aspect('equal')

# Logarithmische Skala
rho_log = np.log10(np.maximum(rho, 1e-6))
cf2 = ax2.contourf(xx, zz, rho_log, levels=30, cmap='RdYlBu_r')
plt.colorbar(cf2, ax=ax2, label='log₁₀[ρ(r)]')
for i, (sym, coord) in enumerate([(mol.atom_symbol(i), atom_coords_bohr[i]) for i in range(mol.natm)]):
    ax2.plot(coord[0], coord[2], 'ko', ms=8)
    ax2.text(coord[0]+0.1, coord[2]+0.1, sym, color='black', fontsize=10, fontweight='bold')
ax2.set_xlabel('x / Bohr')
ax2.set_ylabel('z / Bohr')
ax2.set_title('Elektronendichte ρ(r) – log Skala')
ax2.set_aspect('equal')

plt.suptitle('H₂O B3LYP/6-31G* – Elektronendichte (xz-Ebene)')
plt.tight_layout()
plt.show()

# Gesamtanzahl Elektronen überprüfen (Integral der Dichte)
# Grobe Schätzung via Riemann-Summe
dx_val = x_vals[1] - x_vals[0]
dz_val = z_vals[1] - z_vals[0]
# Achtung: dies ist nur der Beitrag aus der xz-Schnittebene (×1 Bohr Tiefe)
print(f"Integral der Dichte in der xz-Ebene (×1 Bohr): {np.sum(rho)*dx_val*dz_val:.2f} Elektronen")
print(f"Gesamtelektronen (PySCF): {mol.nelectron}")

Integral der Dichte in der xz-Ebene (×1 Bohr): 11.49 Elektronen
Gesamtelektronen (PySCF): 10

10.4 Aufgabe 3 · Dispersionskorrektur (DFT-D3)

Standard-DFT-Funktionale beschreiben Van-der-Waals-Wechselwirkungen schlecht. Die Dispersionskorrektur nach Grimme (D3) behebt dies semiempirisch.

Aufgabe: Berechne die Wechselwirkungsenergie des Wasserdimers (H₂O···H₂O) mit und ohne D3-Korrektur.

# Wasserdimer (Geometrie nahe Minimum)
dimer_geom = """
O  0.000  0.000  0.000
H  0.757  0.586  0.000
H -0.757  0.586  0.000
O  0.000  0.000  2.900
H  0.757  0.586  2.900
H -0.757  0.586  2.900
"""

monomer_geom = "O 0 0 0; H 0.757 0.586 0; H -0.757 0.586 0"

# Ohne D3
mol_dim  = gto.Mole(atom=dimer_geom, basis='6-31g*', verbose=0); mol_dim.build()
mol_mon  = gto.Mole(atom=monomer_geom, basis='6-31g*', verbose=0); mol_mon.build()

mf_dim_b3lyp = dft.RKS(mol_dim, xc='b3lyp'); mf_dim_b3lyp.verbose = 0
E_dim   = mf_dim_b3lyp.kernel()

mf_mon_b3lyp = dft.RKS(mol_mon, xc='b3lyp'); mf_mon_b3lyp.verbose = 0
E_mon   = mf_mon_b3lyp.kernel()

E_int_b3lyp = (E_dim - 2 * E_mon) * 2625.5  # kJ/mol

# Mit D3 (falls dftd3 installiert, sonst manuelle Näherung)
try:
    from pyscf.dftd3 import dftd3
    mf_dim_d3 = dft.RKS(mol_dim, xc='b3lyp')
    mf_dim_d3 = dftd3(mf_dim_d3)
    mf_dim_d3.verbose = 0
    E_dim_d3 = mf_dim_d3.kernel()
    
    mf_mon_d3 = dft.RKS(mol_mon, xc='b3lyp')
    mf_mon_d3 = dftd3(mf_mon_d3)
    mf_mon_d3.verbose = 0
    E_mon_d3 = mf_mon_d3.kernel()
    
    E_int_d3 = (E_dim_d3 - 2 * E_mon_d3) * 2625.5
    d3_label = 'B3LYP-D3'
except Exception:
    E_int_d3 = None
    d3_label = 'B3LYP-D3 (nicht verfügbar)'

# Experimentell: ~21 kJ/mol (SAPT Benchmark)
E_ref = -21.0  # kJ/mol (Bindungsenergie, negativ weil stabilisierend)

print("Wechselwirkungsenergie des Wasserdimers:")
print(f"  B3LYP/6-31G*: {E_int_b3lyp:.2f} kJ/mol")
if E_int_d3 is not None:
    print(f"  B3LYP-D3:     {E_int_d3:.2f} kJ/mol")
else:
    print(f"  {d3_label}")
print(f"  Referenz (SAPT): ~{E_ref} kJ/mol")
Wechselwirkungsenergie des Wasserdimers:
  B3LYP/6-31G*: 12.42 kJ/mol
  B3LYP-D3 (nicht verfügbar)
  Referenz (SAPT): ~-21.0 kJ/mol

10.5 Aufgabe 4 · Austausch-Korrelations-Potential

Visualisiere das Austausch-Korrelations-Potential \(v_{xc}(r)\) entlang einer Linie durch das H₂-Molekül.

from pyscf.dft import numint

# H₂ entlang der z-Achse
mol_h2 = gto.Mole(
    atom="H 0 0 -0.371; H 0 0 0.371",
    basis='cc-pvdz', verbose=0
)
mol_h2.build()

z_line = np.linspace(-4, 4, 300)
grid_h2 = np.column_stack([np.zeros(300), np.zeros(300), z_line])

# Nur LDA-Funktionale: brauchen nur rho (kein Gradient nötig)
# GGA/Hybrid benötigen rho + Gradienten (aufwändigere Berechnung)
xc_functionals = {'LDA (VWN)': 'lda,vwn', 'LDA (PW)': 'lda,pw', 'LDA (PZ)': 'lda,pz'}

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

ni = numint.NumInt()
for xc_name, xc in xc_functionals.items():
    mf = dft.RKS(mol_h2, xc=xc)
    mf.verbose = 0
    mf.kernel()
    dm = mf.make_rdm1()

    # AO-Werte auf Linie → Dichte
    ao_v = mol_h2.eval_gto('GTOval_sph', grid_h2)
    rho_line = np.einsum('pi,ij,pj->p', ao_v, dm, ao_v)
    rho_line = np.maximum(rho_line, 1e-12)  # numerische Stabilität

    # LDA eval_xc: rho muss 1D sein
    exc, vxc, _, _ = ni.eval_xc(xc, rho_line, spin=0)
    vxc_vals = vxc[0]  # v_xc für spin-up = spin-down (RKS)

    axes[0].plot(z_line, rho_line, label=xc_name, lw=2)
    axes[1].plot(z_line, vxc_vals, label=xc_name, lw=2)

for ax in axes:
    ax.axvline(-0.371*1.889, ls=':', color='gray', alpha=0.6, label='H-Kern')
    ax.axvline(+0.371*1.889, ls=':', color='gray', alpha=0.6)
    ax.set_xlabel('z / Bohr')
    ax.legend(fontsize=9)

axes[0].set_title('Elektronendichte ρ(z)')
axes[0].set_ylabel('ρ / e Bohr⁻³')
axes[1].set_title('xc-Potential v_xc(z) – LDA-Varianten')
axes[1].set_ylabel('v_xc / Hartree')
plt.suptitle('H₂ – Dichte und xc-Potential (LDA) entlang Bindungsachse')
plt.tight_layout()
plt.show()

10.5.1 Reflexionsfragen

  1. Warum beschreibt LDA die Elektronendichte von Molekülen dennoch so gut, obwohl es nur vom lokalen Wert \(\rho(\mathbf{r})\) abhängt?
  2. Was ist der konzeptuelle Unterschied zwischen HF und DFT bezüglich der gesuchten Grösse (Wellenfunktion vs. Elektronendichte)?
  3. Warum ist der HOMO-LUMO-Gap bei DFT-Funktionalen systematisch kleiner als bei HF?
  4. Welchen Vorteil bieten Hybrid-Funktionale wie B3LYP gegenüber reinen GGAs, und warum sind sie rechenaufwändiger?

10.6 Zusammenfassung

In dieser Übung hast du:

  • LDA, GGA (PBE/BLYP) und Hybrid-Funktionale (B3LYP, PBE0) für H₂O verglichen
  • Die Elektronendichte auf einem 2D-Gitter visualisiert
  • Van-der-Waals-Wechselwirkungen im Wasserdimer untersucht
  • Das xc-Potential verschiedener Funktionale verglichen