< hartree_fock.py >

From Fortran to Python - Running inside your browser.

input (h2.in)

Format: Title / Dist(Bohr) / Gaussian Basis Format

output (h2.out)

Waiting for execution...

Was macht der Code eigentlich?

Dieser kleine Rechner ist ein direkter Port eines ursprünglich in Fortran 95 geschriebenen Programms, um das HF-Programm direkt im Browser laufen zu lassen. Dazu wird PyScript verwendet. Teile des Codes sind etwas unkonventionell strukturiert, weil der Python-Code direkt auf einem Fortran-Programm basiert.

1. Der Input: Gaussian-Format und Parser

Das Programm parst zunächst die Input-Datei (welche für das Fortran-Programm als h2.in gespeichert wurde). Diese beginnt mit einer Kommentar-Zeile (welche übersprungen wird), in der darauffolgenden Zeile steht die Bindungslänge (in atomaren Einheiten, also Bohr und NICHT in Ångström). Es folgt eine weitere Kommentarzeile (in Fortran-Logik auskommentiert durch ein "!"), dann wird der Basissatz im Gaussian-Format eingelesen, so, wie es beispielsweise auf basissetexchange.org zu finden ist. Eine typische Zeile in diesem Format sieht wie folgt aus:

0.1873113696D+02

Das ist eine alte Schreibweise für Gleitkommazahlen mit doppelter Genauigkeit; der Buchstabe D steht für den Exponenten zur Basis 10. D+02 bedeutet $10^2$. Die Zahl entspricht somit 18.73113696. Da moderne Programmiersprachen wie Python (selbst Excel!!) standardmäßig das "E" verwenden (z.B. 1.8E+01), führt der Parser eine Text-Ersetzung durch, bevor die Zahlen als numerische Werte in die dynamischen Arrays geladen werden.

2. Speicherverwaltung: Statisch vs. Dynamisch

Wichtig im ursprünglichen Code war die Verwaltung des Arbeitsspeichers. In Fortran 77 mussten Array-Größen quasi zur Kompilierzeit feststehen (wir befinden uns hier noch in der Lochkarten-Ära). Mit Fortran 90 wurde die dynamische Speicherverwaltung (ALLOCATABLE) eingeführt. Python macht das praktischerweise automatisch: Listen und NumPy-Arrays passen sich dynamisch der Größe des eingelesenen Basissatzes an, egal ob STO-3G oder 6-311++G verwendet wird. Anzumerken ist, das die aktuelle Version des Codes nur s-Orbitale einlesen kann.

3. It's Submarine Subroutine Time!

Die Berechnung der 1- und 2-Elektronenintegrale ist der eigentlich spannende Teil im HF. Da die Basissätze (wie z.B. STO-3G oder 6-311++G) aus kontrahierten Gauss-Funktionen bestehen, können die Integrale nicht direkt berechnet werden. Stattdessen wird in einer Schleifenstruktur über die Basisfunktionen iteriert und innerhalb dieser über die zugrundeliegenden primitiven Gauss-Funktionen ($\phi_{prim} \propto e^{-\alpha r^2}$). Die Berechnung erfolgt rein analytisch unter Ausnutzung des Gaussian Product Theorems.

Der Trick:
Das Produkt zweier Gauss-Funktionen an den Zentren A und B ist mathematisch äquivalent zu einer einzelnen neuen Gauss-Funktion an einem Punkt P auf der Verbindungslinie zwischen A und B. Dies reduziert Mehrzentren-Integrale auf Integrale über einen einzigen Punkt (Kurzübersicht siehe hier).

Die 1-Elektronen-Integrale

Diese Integrale beschreiben die Bewegung eines einzelnen Elektrons im Feld der Kerne. Sie werden in der Funktion integrale_berechnen berechnet und füllen die Matrizen S, T und V. Jedes Matrixelement (z.B. $S_{\mu\nu}$) ist die Summe über alle Primitiven $i$ (aus $\mu$) und $j$ (aus $\nu$):

$S_{\mu\nu} = \sum_{i,j} d_i d_j N_i N_j \cdot \text{overlap}( \alpha_i, \alpha_j, ... )$

  • overlap_func (Überlappung S):
    Berechnet $\langle \phi_a | \phi_b \rangle$. Eingabe: Zwei Exponenten ($\alpha, \beta$) und zwei Kernkoordinaten ($R_A, R_B$). Da s-Orbitale kugelsymmetrisch sind, hängt das Integral nur vom Abstand der Zentren und der Summe der Exponenten ab. Das Ergebnis ist ein Skalar zwischen 0 und 1.
  • kinetic_func (Kinetische Energie T):
    Berechnet $\langle \phi_a | -\frac{1}{2}\nabla^2 | \phi_b \rangle$. Hier wird der Laplace-Operator angewendet. Für s-Gauss-Funktionen lässt sich dieses Integral algebraisch auf eine Kombination aus dem Überlappungsintegral und den Exponenten zurückführen. Es repräsentiert die Energie der Elektronenbewegung.
  • nuc_att_func (Kernanziehung V):
    Berechnet $\langle \phi_a | \frac{-Z_C}{|r - R_C|} | \phi_b \rangle$. Dies ist das komplizierteste 1-Elektronen-Integral. Es beschreibt die Coulomb-Anziehung zwischen der Elektronenwolke (gebildet aus dem Produkt von $\phi_a$ und $\phi_b$) und einem Atomkern $C$. Der Code muss über alle Atome im Molekül summieren.
    Implementierung: Da der Operator $1/r$ eine Singularität bei $r=0$ hat, wird die Boys-Funktion $F_0(x)$ benötigt. Diese Funktion glättet die Singularität mathematisch durch Nutzung der Fehlerfunktion: $F_0(x) = \text{erf}(\sqrt{x})/\sqrt{x}$.

Die 2-Elektronen-Integrale (ERI)

Das Integral eri_calc_func berechnet die Abstoßung zwischen zwei Elektronenwolken: $(\mu\nu | \lambda\sigma) = \int \int \phi_\mu(1)\phi_\nu(1) \frac{1}{r_{12}} \phi_\lambda(2)\phi_\sigma(2) d\tau_1 d\tau_2$.

  • Struktur: Hier liegt der Bottleneck der Rechnung. Der Code benötigt vier ineinander verschachtelte Schleifen über die Basisfunktionen (mu, nu, lam, sig). Innerhalb jeder dieser Schleifen laufen weitere Schleifen über die Primitiven. Das führt zu einer formalen Skalierung von $N^4$.
  • Berechnung:
    1. Das Produkt $\phi_\mu \phi_\nu$ wird zu einer Gauss-Verteilung am Punkt $P$ zusammengefasst.
    2. Das Produkt $\phi_\lambda \phi_\sigma$ wird zu einer Gauss-Verteilung am Punkt $Q$ zusammengefasst.
    3. Das Problem reduziert sich auf die Coulomb-Wechselwirkung zwischen zwei Ladungsverteilungen an $P$ und $Q$.
  • Die Boys-Funktion: Auch hier kommt die Hilfsfunktion boys(x) zum Einsatz, um den $1/r_{12}$ Operator zu integrieren. Der Code nutzt scipy.special.erf für maximale numerische Präzision, fängt aber den Fall $x \approx 0$ (kleine Abstände) durch eine Taylor-Entwicklung ($1 - x/3$) ab, um Divisionen durch Null zu vermeiden.

4. Der Hartree-Fock Algorithmus in 12 Schritten

Der eigentliche Programmablauf folgt exakt den 12 logischen Schritten aus der Vorlesung:

Schritt 1: Spezifikation des Moleküls
Die Geometrie wird festgelegt (H bei 0.0 und H bei R). Die Kernabstoßungsenergie $E_{nuc} = 1/R$ wird berechnet.

Schritt 2: Integralberechnung
Alle Ein-Elektronen- (S, T, V) und Zwei-Elektronen-Integrale (ERI) werden berechnet und im Arbeitsspeicher abgelegt. $H_{core}$ wird als Summe von T und V gebildet.

Schritt 3: Diagonalisierung der Überlappungsmatrix S
Um das Problem in ein orthogonales Eigenwertproblem zu überführen, wird die Transformationsmatrix $X = S^{-1/2}$ berechnet.

Schritt 4: Initial Guess
Die Dichtematrix P wird initialisiert (hier $P=0$).

--- START SCF ZYKLUS ---

Schritt 5 & 6: Aufbau der Fock-Matrix
Die Fock-Matrix F wird aus dem Ein-Elektronen-Teil $H_{core}$ und dem Zwei-Elektronen-Teil $G(P)$ zusammengesetzt. G hängt von der aktuellen Dichtematrix P ab.

Schritt 7: Transformation
Die Fock-Matrix wird in die orthogonale Basis transformiert: $F' = X^T F X$.

Schritt 8: Diagonalisierung von F'
Hier wird die Schrödingergleichung $F'C' = C'\epsilon$ gelöst. In Fortran 95 wurde das durch eine manuelle Jacobi-Rotation gelöst. In Python übernimmt das netterweise numpy.linalg.eigh (basierend auf LAPACK), was effizienter ist (und mir beim porten Kopfschmerzen erspart hat). Dies liefert die Orbitalenergien $\epsilon$ und die Koeffizienten $C'$.

Schritt 9: Rücktransformation
Die Koeffizienten werden in die ursprüngliche Atomorbital-Basis zurücktransformiert: $C = X C'$.

Schritt 10: Bildung der neuen Dichte
Aus den besetzten Orbitalen in C wird eine neue Dichtematrix $P_{neu}$ konstruiert.

Schritt 11: Konvergenzprüfung
$P_{neu}$ wird mit der Dichte aus dem vorherigen Schritt verglichen. Wenn die Änderung kleiner als $10^{-11}$ ist, gilt das Feld als konvergiert. Wenn nicht, springen wir zurück zu Schritt 5.

--- ENDE SCF ZYKLUS ---

Schritt 12: Finale Energieberechnung
Die elektronische Gesamtenergie wird aus der konvergierten Dichte und den Integralen berechnet. Addiert man die Kernabstoßung $E_{nuc}$, erhält man die totale Energie des Systems in Hartree.

Zurück zu Bleistift-und-Papier Quantenmechanik

Da es sich komisch anfühlt, den rohen Code dieser Website noch mal auf die Website selbst zu setzen, hier stattdessen der Versuch, den Fortran-turned-Python-Code wieder zu Bleistift-und-Papier Quantenmechanik zu übersetzen. Wir betrachten das System, das oben via dem Beispiel-Input gegeben ist: H$_2$ ($R=1.4$ Bohr), aber mit der minimalen STO-3G-Basis. Die Indizes $\mu, \nu, \lambda, \sigma$ laufen über die kontrahierten Basisfunktionen (AO) $\chi$.

1. System-Initialisierung

E_nuc = 1.0 / R_abstand
$$ V_{nn} = \sum_{A>B} \frac{Z_A Z_B}{|\boldsymbol{R}_A - \boldsymbol{R}_B|} = \frac{1 \cdot 1}{1.4} \approx 0.7142 \, E_h $$
basis_alpha, basis_d (Initialisierung STO-3G)

Konstruktion der Basisfunktionen $\chi_\mu$ durch Kontraktion von 3 primitiven Gauß-Funktionen ($K=3$). Für das Wasserstoffatom $A$ bei $\boldsymbol{R}_A$:

$$ \chi_1(\boldsymbol{r}) = \sum_{k=1}^{3} d_{k1} \left(\frac{2\alpha_{k1}}{\pi}\right)^{3/4} e^{-\alpha_{k1} |\boldsymbol{r} - \boldsymbol{R}_A|^2} $$

Hier entsprechen basis_alpha den Exponenten $\alpha$ und basis_d den Kontraktionskoeffizienten $d$.

2. Berechnung der Integrale (Pre-Loop)

S[m,n] = val_S (in integrale_berechnen)
$$ S_{\mu\nu} = \langle \chi_\mu | \chi_\nu \rangle = \int \chi_\mu^*(\boldsymbol{r}) \chi_\nu(\boldsymbol{r}) \, d\boldsymbol{r} $$
T_int[m,n] = val_T
$$ T_{\mu\nu} = \langle \chi_\mu | -\frac{1}{2} \nabla^2 | \chi_\nu \rangle $$
V_int[m,n] = val_V
$$ V_{\mu\nu} = \langle \chi_\mu | \hat{V}_{ne} | \chi_\nu \rangle = \sum_A \langle \chi_\mu | -\frac{Z_A}{|\boldsymbol{r} - \boldsymbol{R}_A|} | \chi_\nu \rangle $$
H_core[m,n] = val_T + val_V
$$ H_{\mu\nu}^{\text{core}} = T_{\mu\nu} + V_{\mu\nu} $$
ERI[m,n,la,si] (4-fach Schleife)
$$ (\mu\nu|\lambda\sigma) = \iint \chi_\mu^*(1)\chi_\nu(1) \frac{1}{r_{12}} \chi_\lambda^*(2)\chi_\sigma(2) \, d\boldsymbol{r}_1 d\boldsymbol{r}_2 $$

Die Funktion eri_calc_func nutzt hier explizit Boys-Funktionen, um das $1/r_{12}$ Integral analytisch zu lösen.

3. Orthogonalisierung (Löwdin)

epsilon_S, U = np.linalg.eigh(S)

Lösen des Eigenwertproblems der Überlappungsmatrix:

$$ \boldsymbol{S}\boldsymbol{U} = \boldsymbol{U}\boldsymbol{s} \quad (\boldsymbol{s} \text{ ist diagonal}) $$
X = U @ diag_inv_sqrt_s @ U.T

Konstruktion der Transformationsmatrix $\boldsymbol{X} = \boldsymbol{S}^{-1/2}$:

$$ \boldsymbol{X} = \boldsymbol{U} \boldsymbol{s}^{-1/2} \boldsymbol{U}^T $$

Dies transformiert die nicht-orthogonale AO-Basis in eine orthogonale Basis.

4. Der SCF-Zyklus (RHF)

for iteration in range(1, 201):

Beginn des iterativen Verfahrens zur Lösung der Roothaan-Hall-Gleichungen (und gleichzeitig Begrenzung der Maximalen Zyklenanzahl).

F[mu,nu] += P[lam,sig] * (ERI... - 0.5 * ERI...)

Konstruktion der Fock-Matrix $\boldsymbol{F}$ aus dem Einteilchen-Hamiltonian und dem mittleren Feld der Elektronen (Coulomb $\boldsymbol{J}$ und Austausch $\boldsymbol{K}$):

$$ F_{\mu\nu} = H_{\mu\nu}^{\text{core}} + \sum_{\lambda\sigma} P_{\lambda\sigma} \left[ (\mu\nu|\lambda\sigma) - \frac{1}{2} (\mu\lambda|\nu\sigma) \right] $$

Dies entspricht in Matrixform $\boldsymbol{F} = \boldsymbol{H}^{\text{core}} + \boldsymbol{G}(\boldsymbol{P})$.

F_strich = X.T @ F @ X

Transformation in die orthogonale Basis:

$$ \boldsymbol{F}' = \boldsymbol{X}^T \boldsymbol{F} \boldsymbol{X} $$
eps, C_strich = np.linalg.eigh(F_strich)

Diagonalisierung der Fock-Matrix (Lösung der Pseudo-Eigenwertgleichung):

$$ \boldsymbol{F}' \boldsymbol{C}' = \boldsymbol{C}' \boldsymbol{\epsilon} $$

$\boldsymbol{\epsilon}$ enthält die Orbitalenergien, $\boldsymbol{C}'$ die MO-Koeffizienten in der orthogonalen Basis.

C = X @ C_strich

Rücktransformation in die ursprüngliche AO-Basis:

$$ \boldsymbol{C} = \boldsymbol{X} \boldsymbol{C}' $$

Das Molekülorbital $\psi_i$ ist nun definiert als $\psi_i = \sum_\mu C_{\mu i} \chi_\mu$.

P_neu[mu,nu] = 2.0 * C[mu,0] * C[nu,0]

Aufbau der neuen Dichtematrix $\boldsymbol{P}$ für einen geschlossenschaligen Singulett-Zustand (RHF). Summation über die besetzten Orbitale (hier nur $a=1$):

$$ P_{\mu\nu} = 2 \sum_{a}^{N_{elec}/2} C_{\mu a} C_{\nu a}^* $$

5. Finale Energieberechnung

E_total += 0.5 * P[i,j] * (H_core[i,j] + F[i,j])

Berechnung der elektronischen Energie als Erwartungswert. Man beachte, dass $\boldsymbol{F} = \boldsymbol{H}^{\text{core}} + \boldsymbol{G}$ ist, wodurch die Doppelzählung der Elektron-Elektron-Wechselwirkung korrigiert wird:

$$ E_{elec} = \frac{1}{2} \text{Spur}\left( \boldsymbol{P} (\boldsymbol{H}^{\text{core}} + \boldsymbol{F}) \right) $$
E_total += E_nuc
$$ E_{tot} = E_{elec} + V_{nn} $$
packages = ["numpy", "scipy"]