import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
plt.rcParams.update({'figure.dpi': 120, 'font.size': 11})
import matplotlib
matplotlib.use('Agg')
try:
from pyscf import gto, scf, dft, grad, hessian
PYSCF_AVAILABLE = True
except ImportError:
PYSCF_AVAILABLE = False
print("PySCF nicht gefunden. Installiere mit: uv add pyscf")11 Ü4 · Molekulare Eigenschaften
Dipolmoment, Geometrieoptimierung und Schwingungsfrequenzen
11.1 Lernziele
Nach dieser Übung kannst du:
- Elektrische Eigenschaften (Dipolmoment, Polarisierbarkeit) berechnen und interpretieren
- Geometrieoptimierungen durchführen und konvergierte Strukturen analysieren
- Potentialenergiekurven (PES) für einfache Systeme berechnen
- Mulliken- und Löwdin-Ladungen vergleichen
- Schwingungsfrequenzen aus numerischen Hessian-Matrizen bestimmen
11.2 Aufgabe 1 · Dipolmomente
Das Dipolmoment \(\boldsymbol{\mu}\) ist eine der wichtigsten molekularen Eigenschaften: \[\mu_\alpha = -\sum_i \langle\psi_i|r_\alpha|\psi_i\rangle + \sum_A Z_A R_{A\alpha}\]
Berechne Dipolmomente für eine Reihe polarer und unpolarer Moleküle.
if not PYSCF_AVAILABLE:
raise SystemExit("PySCF benötigt")
# 1 Debye = 0.393430 a.u.
AU_TO_DEBYE = 2.5418
molecules = {
'H₂O': ('O 0 0 0; H 0.757 0.586 0; H -0.757 0.586 0', 0, 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', 0, 0),
'HF': ('F 0 0 0; H 0 0 0.917', 0, 0),
'HCl': ('Cl 0 0 0; H 0 0 1.275', 0, 0),
'CO': ('C 0 0 0; O 0 0 1.128', 0, 0),
'H₂': ('H 0 0 -0.371; H 0 0 0.371', 0, 0),
'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', 0, 0),
}
# Experimentelle Dipolmomente in Debye
dip_exp = {'H₂O': 1.85, 'NH₃': 1.47, 'HF': 1.83, 'HCl': 1.11, 'CO': 0.11, 'H₂': 0.00, 'CH₄': 0.00}
results_dip = []
for name, (atom_str, charge, spin) in molecules.items():
mol = gto.Mole(atom=atom_str, basis='6-31g*', charge=charge, spin=spin, verbose=0)
mol.build()
mf = dft.RKS(mol, xc='b3lyp')
mf.verbose = 0
mf.kernel()
dip = mf.dip_moment(verbose=0) # in Debye
dip_norm = np.linalg.norm(dip)
results_dip.append({
'Molekül': name,
'|μ| calc [D]': dip_norm,
'|μ| exp [D]': dip_exp[name],
'Fehler [D]': dip_norm - dip_exp[name],
})
print(f"{name:6s}: μ = {dip_norm:.3f} D (exp: {dip_exp[name]:.2f} D)")
df_dip = pd.DataFrame(results_dip)
print()
print(df_dip.to_string(index=False))H₂O : μ = 2.073 D (exp: 1.85 D)
NH₃ : μ = 1.866 D (exp: 1.47 D)
HF : μ = 1.842 D (exp: 1.83 D)
HCl : μ = 1.461 D (exp: 1.11 D)
CO : μ = 0.097 D (exp: 0.11 D)
H₂ : μ = 0.000 D (exp: 0.00 D)
CH₄ : μ = 0.000 D (exp: 0.00 D)
Molekül |μ| calc [D] |μ| exp [D] Fehler [D]
H₂O 2.073118e+00 1.85 2.231180e-01
NH₃ 1.865791e+00 1.47 3.957909e-01
HF 1.842358e+00 1.83 1.235794e-02
HCl 1.460798e+00 1.11 3.507978e-01
CO 9.657370e-02 0.11 -1.342630e-02
H₂ 2.072775e-15 0.00 2.072775e-15
CH₄ 3.248640e-14 0.00 3.248640e-14fig, ax = plt.subplots(figsize=(9, 5))
x = np.arange(len(df_dip))
ax.bar(x - 0.2, df_dip['|μ| calc [D]'], 0.35, label='B3LYP/6-31G*', color='steelblue', alpha=0.85)
ax.bar(x + 0.2, df_dip['|μ| exp [D]'], 0.35, label='Experiment', color='darkorange', alpha=0.85)
ax.set_xticks(x)
ax.set_xticklabels(df_dip['Molekül'], fontsize=12)
ax.set_ylabel('Dipolmoment |μ| / Debye')
ax.set_title('Dipolmomente: B3LYP/6-31G* vs. Experiment')
ax.legend()
plt.tight_layout()
plt.show()11.3 Aufgabe 2 · Mulliken- und Löwdin-Ladungen
Populationsanalysen ordnen Elektronen formal den Atomen zu. Vergleiche Mulliken- und Löwdin-Ladungen.
from pyscf.lo import orth
def lowdin_charges(mol, dm):
"""Löwdin-Populationsanalyse: orthogonalisierte Basis via S^{-1/2}."""
S = mol.intor('int1e_ovlp')
# S^{1/2}
eigval, eigvec = np.linalg.eigh(S)
S_half = eigvec @ np.diag(np.sqrt(eigval)) @ eigvec.T
# Dichtematrix in orthogonaler Basis
dm_orth = S_half @ dm @ S_half
# Populationen pro AO
ao_pop = np.diag(dm_orth)
# Summierung über Schalen pro Atom
atom_charges = np.array([mol.atom_charge(i) for i in range(mol.natm)], dtype=float)
for i, label in enumerate(mol.ao_labels(fmt=None)):
atom_idx = label[0]
atom_charges[atom_idx] -= ao_pop[i]
return atom_charges
test_molecules = {
'H₂O': "O 0 0 0; H 0.757 0.586 0; H -0.757 0.586 0",
'HF': "F 0 0 0; H 0 0 0.917",
'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",
}
for mol_name, atom_str in test_molecules.items():
mol = gto.Mole(atom=atom_str, basis='6-31g*', verbose=0)
mol.build()
mf = dft.RKS(mol, xc='b3lyp')
mf.verbose = 0
mf.kernel()
dm = mf.make_rdm1()
# Mulliken
mull_pop, mull_charges = mf.mulliken_pop(verbose=0)
# Löwdin
low_charges = lowdin_charges(mol, dm)
print(f"\n{mol_name} (B3LYP/6-31G*):")
print(f" {'Atom':>6} | {'Mulliken':>10} | {'Löwdin':>10}")
print(" " + "-" * 32)
for i in range(mol.natm):
print(f" {mol.atom_symbol(i):>6} | {mull_charges[i]:>10.4f} | {low_charges[i]:>10.4f}")
H₂O (B3LYP/6-31G*):
Atom | Mulliken | Löwdin
--------------------------------
O | -0.7978 | -0.6158
H | 0.3989 | 0.3079
H | 0.3989 | 0.3079
HF (B3LYP/6-31G*):
Atom | Mulliken | Löwdin
--------------------------------
F | -0.4878 | -0.3800
H | 0.4878 | 0.3800
NH₃ (B3LYP/6-31G*):
Atom | Mulliken | Löwdin
--------------------------------
N | -0.8983 | -0.6704
H | 0.2993 | 0.2235
H | 0.2995 | 0.2235
H | 0.2995 | 0.223511.4 Aufgabe 3 · Potentialenergiekurve: H₂-Bindungsdissoziation
Berechne die Potentialenergiekurve (PES) für die H-H-Streckbewegung und vergleiche HF mit DFT.
r_vals = np.linspace(0.5, 5.0, 30) # H-H-Abstand in Angström
E_hf = np.zeros(len(r_vals))
E_b3lyp = np.zeros(len(r_vals))
for i, r in enumerate(r_vals):
mol = gto.Mole(
atom=f"H 0 0 0; H 0 0 {r}",
basis='cc-pvdz', verbose=0, spin=0
)
mol.build()
# RHF
mf = scf.RHF(mol)
mf.verbose = 0
E_hf[i] = mf.kernel()
# B3LYP
mf_dft = dft.RKS(mol, xc='b3lyp')
mf_dft.verbose = 0
E_b3lyp[i] = mf_dft.kernel()
# Dissoziation: Energie relativ zum Gleichgewichtsabstand
E_hf_rel = (E_hf - E_hf.min()) * 2625.5 # kJ/mol
E_b3lyp_rel = (E_b3lyp - E_b3lyp.min()) * 2625.5 # kJ/mol
# Gleichgewichtsabstand
r_eq_hf = r_vals[E_hf.argmin()]
r_eq_b3lyp = r_vals[E_b3lyp.argmin()]
print(f"Gleichgewichtsabstand RHF: {r_eq_hf:.3f} Å")
print(f"Gleichgewichtsabstand B3LYP: {r_eq_b3lyp:.3f} Å")
print(f"Experiment: 0.741 Å")
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(r_vals, E_hf_rel, 'o-', color='steelblue', label='RHF/cc-pVDZ', lw=2)
ax.plot(r_vals, E_b3lyp_rel, 's-', color='darkorange', label='B3LYP/cc-pVDZ', lw=2)
ax.axvline(0.741, ls='--', color='gray', alpha=0.7, label='r_exp = 0.741 Å')
ax.set_xlabel('H-H-Abstand / Å')
ax.set_ylabel('Relative Energie / kJ mol⁻¹')
ax.set_title('H₂-Potentialenergiekurve: RHF vs. B3LYP')
ax.set_ylim(-20, 500)
ax.legend()
plt.tight_layout()
plt.show()
print("\nNote: RHF bricht die Dissoziation nicht korrekt ab (ionische Grenzstruktur dominiert)")
print("DFT beschreibt die Dissoziation ebenfalls nicht exakt korrekt (statische Korrelation)")Gleichgewichtsabstand RHF: 0.810 Å
Gleichgewichtsabstand B3LYP: 0.810 Å
Experiment: 0.741 Å
Note: RHF bricht die Dissoziation nicht korrekt ab (ionische Grenzstruktur dominiert)
DFT beschreibt die Dissoziation ebenfalls nicht exakt korrekt (statische Korrelation)11.5 Aufgabe 4 · Geometrieoptimierung
Optimiere die Geometrie von H₂O mit B3LYP/6-31G* und vergleiche die Ergebnisse mit dem Experiment.
try:
from pyscf.geomopt.berny_solver import optimize
GEOMOPT = True
except ImportError:
try:
from pyscf.geomopt.geometric_solver import optimize
GEOMOPT = True
except ImportError:
GEOMOPT = False
print("Geometrieoptimierung benötigt: uv add geometric")
if GEOMOPT:
# Startgeometrie (leicht gestört vom Gleichgewicht)
mol = gto.Mole(
atom="O 0 0 0; H 0.800 0.600 0; H -0.800 0.600 0",
basis='6-31g*', verbose=0
)
mol.build()
mf = dft.RKS(mol, xc='b3lyp')
mf.verbose = 0
mol_opt = optimize(mf)
# Optimierte Koordinaten (in Angström)
coords_opt = mol_opt.atom_coords() * 0.529177 # Bohr → Angström
# O-H-Abstand und H-O-H-Winkel berechnen
O = coords_opt[0]
H1 = coords_opt[1]
H2 = coords_opt[2]
r_OH1 = np.linalg.norm(H1 - O)
r_OH2 = np.linalg.norm(H2 - O)
v1 = (H1 - O) / r_OH1
v2 = (H2 - O) / r_OH2
angle = np.degrees(np.arccos(np.clip(np.dot(v1, v2), -1, 1)))
print("Optimierte Geometrie (B3LYP/6-31G*):")
print(f" r(O-H): {r_OH1:.4f} Å (exp: 0.9572 Å)")
print(f" ∠(H-O-H): {angle:.2f}° (exp: 104.52°)")
else:
print("geometric nicht installiert – Geometrieoptimierung übersprungen.")
print("Installiere mit: uv add geometric")Geometrieoptimierung benötigt: uv add geometric
geometric nicht installiert – Geometrieoptimierung übersprungen.
Installiere mit: uv add geometric11.6 Aufgabe 5 · Schwingungsfrequenz: H₂ via numerischem Hessian
Berechne die Streckschwingungsfrequenz von H₂ via numerischer zweiter Ableitung der Energie.
# Numerischer Hessian für H₂ (nur eine Koordinate: H-H-Abstand)
# E(r) ≈ E(r0) + 1/2 * k * (r - r0)^2 → k = d²E/dr²
def E_h2(r_angstrom, xc='b3lyp', basis='cc-pvdz'):
mol = gto.Mole(
atom=f"H 0 0 0; H 0 0 {r_angstrom}",
basis=basis, verbose=0
)
mol.build()
mf = dft.RKS(mol, xc=xc)
mf.verbose = 0
return mf.kernel()
# Gleichgewichtsabstand (optimiert oder experimentell)
r0 = 0.742 # Angström
dr = 0.005 # Schrittweite
E_plus = E_h2(r0 + dr)
E_0 = E_h2(r0)
E_minus = E_h2(r0 - dr)
# Numerische zweite Ableitung (Hartree / Ų)
k_hartree_ang = (E_plus - 2*E_0 + E_minus) / dr**2
# Einheitenumrechnung → SI
# 1 Hartree/Bohr² = 1556 N/m
# 1 Bohr = 0.529177e-10 m
# 1 Hartree = 4.3597e-18 J
hartree_per_bohr2_to_Nm = 1556.0
ang_to_bohr = 1.0 / 0.529177
k_hartree_bohr2 = k_hartree_ang * ang_to_bohr**2
k_SI = k_hartree_bohr2 * hartree_per_bohr2_to_Nm # N/m
# Reduzierte Masse von H₂ (1 amu)
mu_amu = 0.5 * 1.00794 # amu
mu_kg = mu_amu * 1.66054e-27 # kg
# Kreisfrequenz → Frequenz in cm⁻¹
omega = np.sqrt(k_SI / mu_kg) # rad/s
nu_hz = omega / (2 * np.pi) # Hz
nu_cm = nu_hz / (2.998e10) # cm⁻¹
print(f"Kraftkonstante k: {k_SI:.1f} N/m")
print(f"Schwingungsfrequenz ν̃: {nu_cm:.0f} cm⁻¹")
print(f"Experimentell (H₂): 4161 cm⁻¹")Kraftkonstante k: 8164.8 N/m
Schwingungsfrequenz ν̃: 16582 cm⁻¹
Experimentell (H₂): 4161 cm⁻¹11.6.1 Reflexionsfragen
- Warum ist die Mulliken-Ladung basissatzabhängig, die Löwdin-Ladung dagegen stabiler?
- Warum beschreibt RHF die H₂-Dissoziation falsch? Welche Methode wäre besser geeignet?
- Welche Einflüsse werden bei der harmonischen Näherung für die Schwingungsfrequenz vernachlässigt?
- Das berechnete Dipolmoment von CO ist kleiner als erwartet. Was zeigt das über die Genauigkeit einfacher Basissätze?
11.7 Zusammenfassung
In dieser Übung hast du:
- Dipolmomente für polare und unpolare Moleküle berechnet und mit Experimenten verglichen
- Mulliken- und Löwdin-Ladungen verglichen und Stärken/Schwächen diskutiert
- Eine Potentialenergiekurve für H₂-Dissoziation mit RHF und B3LYP berechnet
- Geometrieoptimierung und Schwingungsfrequenzen numerisch bestimmt