< hartree_fock_h2_v1.py >

From Fortran to Python - Running inside your browser.

input deck (h2.in)

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

system output (stdout)

Waiting for execution...

Von Fortran zu Python

Dieser Rechner ist ein direkter Port eines ursprünglich in Fortran 95 geschriebenen Programms. Das Ziel des Ports war es, die Quantenmechanik quasi direkt und transparent im Browser laufen zu lassen. Der zugrundeliegende Code ist etwas unkonventionell strukturiert, da der Python-Code 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. 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

Dabei handelt es sich um eine alte Schreibweise für Gleitkommazahlen mit doppelter Genauigkeit. Der Buchstabe D steht für den Exponenten zur Basis 10. D+02 bedeutet 102. Die Zahl entspricht somit 18.73113696. Da moderne Programmiersprachen wie Python 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

Ein zentraler Aspekt ist die Verwaltung des Arbeitsspeichers. In Fortran 77 mussten Array-Größen oft zur Kompilierzeit feststehen. Mit Fortran 90 wurde die dynamische Speicherverwaltung (ALLOCATABLE) eingeführt. In diesem Python-Port geschieht dies 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. Die Subroutinen: Berechnung der 1- und 2-Elektronenintegrale

Bevor wir den SCF-Zyklus betrachten, betrachten wir die Integrale, die im Verlauf über die primitiven Gauss-Funktionen berechnet werden. Dies wird von der Funktion integrale_berechnen gehandhabt, welche ihrerseits auf vier praktische Hilfsfunktionen zurückgreift:

  • overlap (S): Berechnet ⟨ μ | ν ⟩. Da Atomorbitale nicht orthogonal sind, ist das Überlappungsintegral zwischen zwei Gauss-Funktionen ungleich Null.
  • kinetic (T): Berechnet die kinetische Energie der Elektronen: ⟨ μ | -½ ∇2 | ν ⟩.
  • nuc_att (V): Berechnet die Anziehung der Elektronen an die Kerne: ⟨ μ | -Z/r | ν ⟩. Dies erfordert die Auswertung der Boys-Funktion F0(x), welche mathematisch auf der Fehlerfunktion basiert: F0(x) = erf(√x) / √x.
  • eri_calc (2e): Berechnet die Elektron-Elektron-Abstoßung (μν|λσ). Dieser Schritt ist der rechenintensivste, da er formal mit der vierten Potenz der Basisgröße skaliert.

4. Der Hartree-Fock Algorithmus in 12 Schritten

Der eigentliche Programmablauf folgt exakt dem Roothaan-Hall-Schema, unterteilt in 12 logische Schritte (strikt gemäß der Vorlesung!):

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

Schritt 2: Integralberechnung
Alle Ein-Elektronen- (S, T, V) und Zwei-Elektronen-Integrale (ERI) werden berechnet und im Arbeitsspeicher abgelegt. Hcore 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-½ 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 Hcore 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' = XT F X.

Schritt 8: Diagonalisierung von F'
Hier wird die Schrödingergleichung F'C' = C'ε gelöst. Im Fortran-Original wurde dies durch eine manuelle Jacobi-Rotation gelöst. In Python übernimmt das netterweise numpy.linalg.eigh (basierend auf LAPACK), was mathematisch effizienter ist. Dies liefert die Orbitalenergien ε 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 Pneu konstruiert.

Schritt 11: Konvergenzprüfung
Pneu wird mit der Dichte aus dem vorherigen Schritt verglichen. Wenn die Änderung kleiner als 10-8 ist, gilt das Feld als konvergiert. Wenn nicht, beginnt der Zyklus erneut bei 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 Enuc, erhält man die totale Energie des Systems in Hartree.
packages = ["numpy", "scipy"]