Quantenmechanik/ Hartree-Fock

Zur Rechentechnik

Für einfachere Formeln in diesem Abschnitt sorgen dimensionslose Längen, Zeiten, Energien. Die Variablen sollen so skaliert werden, dass aus Schrödinger-Gleichungen die Faktoren $q^{2}{\text{ bzw. }}q^{2}/(4\pi \epsilon _{0}){\text{ in SI,}}\;\hbar ^{2}/m,\;\hbar \;$ wegfallen. Es kommen ja nur Elektronen vor mit Ladung q und Masse m.

Definiere $R=r/a;\;T=t/b;\;U=u/c.$ Dimensionsloses in Großbuchstaben.

$u\phi =i\hbar \partial _{t}\phi =(\hbar ^{2}/(2m))\partial _{r}^{2}\phi +(q^{2}/r)\phi \quad \to$

$U\phi =i\partial _{T}\phi =\partial _{R}^{2}\phi /2+\phi /R$

Einsetzen von $u=cU,\;\partial _{t}=\partial _{T}/b,\;r=Ra,\;\partial _{r}=\partial _{R}/a$ ergibt

$cU\phi =(i\hbar /b)\partial _{T}\phi =(\hbar ^{2}/m/a^{2})\partial _{R}^{2}\phi /2+(q^{2}/a)\phi /R$

$\Rightarrow \;c=\hbar /b=(\hbar ^{2}/m/a^{2})=(q^{2}/a)\;\Rightarrow$

Längeneinheit $a=\hbar ^{2}/(mq^{2}),$ Energie-Einheit $c=q^{2}/a=\hbar ^{2}/(ma^{2}),$ Zeiteinheit $b=\hbar /c.$

Mit diesen Bohrschen Einheiten sind Ladung und Masse der Elektronen und das Wirkungsquantum im Folgenden rausgeworfen. Symbole wie ${\vec {r}},E,t$ sind dimensionslos und bedeuten ${\vec {r}}/a$ und so weiter.

Zur Herleitung einiger Näherungsverfahren des Viel-Fermion-Problems werden Elemente der Funktionalrechnung gebraucht, die auch woanders nützlich sind.

Lagrange-Multiplikatoren

Das Problem, extreme Werte einer glatten Funktion f(x) von n Variablen $x=(x_{1},...,x_{n})$ zu finden, führt unweigerlich zur Aufgabe, den Gradienten, also alle partiellen Ableitungen, auf Null zu setzen: $\partial _{i}f=\partial f/\partial x_{i}=0.$
Nun sollen k Nebenbedingungen dazu kommen, etwa $g(x)=0,\;g=(g_{1},...,g_{k}).$ Glatte Nebenbedingungen natürlich, dazu nur mit einfachen Nullstellen. Wir Normalmenschen versuchen zuerst, die Sache auf (n-k) unabhängige Variablen umzurechnen und dann wie im uneingeschränkten Fall fortzufahren. Lagrange war schlauer als normale Menschen und fand einen besseren Trick mit (n+k) unabhängigen Variablen. Mehr statt weniger!

Angenommen, bei $x=x_{0}$ sitzt ein Extremum von f, auf der Teilmenge mit den k Bedingungen $g_{j}(x)=0;\;j=1...k.$ Dann gibt es mindestens (n-k) linear unabhängige Richtungsvektoren $\{u_{h}\}$ am Punkt x und (n-k) Geraden $z_{h}(t):=x_{0}+tu_{h}$ als Tangenten der Mannigfaltigkeit $\{x|g(x)=0\}.$ Für diese verschwinden folgende Ableitungen bei $t=0:$

$0=(d/dt)f(x_{0}+tu_{h})=\sum \nolimits _{i}u_{hi}\;\partial _{i}f(x_{0})$

$0=(d/dt)g_{j}(x_{0}+tu_{h})=\sum \nolimits _{i}u_{hi}\;\partial _{i}g_{j}(x_{0})$

Die Gradienten von f und von den k Bedingungen stehen also senkrecht auf den (n-k) unabhängigen Richtungsvektoren. Insgesamt sind das (k+1) Gradienten der Dimension n. Folglich müssen sie linear abhängig sein und es gibt Faktoren $\{z_{j}\},$ so dass $\partial _{i}f(x_{0})+\sum \nolimits _{j}z_{j}\;\partial _{i}g_{j}(x_{0})=0.$

Nun definiert man die Funktion von (n+k) Variablen ${\bar {f}}(x,y):=f(x)+\sum \nolimits _{j}y_{j}\;g_{j}(x)$
Am Punkt $(x,y)=(x_{0},z)$ ist ihr Gradient Null, denn:

$\partial {\bar {f}}/\partial x_{i}=\partial _{i}f+\sum \nolimits _{j}z_{j}\;\partial _{i}g_{j}=0;\quad \partial {\bar {f}}/\partial y_{j}=g_{j}(x)=0.$

Es folgt der Lagrange'sche Algorithmus zu Optimierungen und so weiter: Zur Funktion f(x) mit gesuchtem Extremum werden die Nebenbedingungen $g_{j}(x)$ mit Lagrange-Multiplikatoren $y_{j}$ addiert. Dann wird nach unbeschränkten Extrema im Definitionsbereich von $(x,y)$ gesucht.

Multiplikator-Methode mit Funktional-Ableitungen

Für nichtlineare Funktionale gibt es eine Differenzialrechnung und die häufige Suche nach stationären oder extremen Werten. Das Wirkungprinzip für Bahnkurven und für Felder jeder Art ist Anwendungsfall Nummer Eins in der Physik.

Ein Funktional A sei hier eine Abbildung von einem Definitionsbereich von glatten, stark abfallenden Funktionen in den Wertebereich der reellen Zahlen. Die Funktionen sollen vom Typ $f:\mathbb {R} ^{n}\to \mathbb {R} ^{m}$ sein und die Funktionale sind lokal oder bi-lokal oder komplizierter. Etwa so:

$A[f]=\int d^{n}x\;L(x,f(x),\partial _{i}f(x),\partial _{i}\partial _{k}f(x))$

$A[f]=\int d^{n}x\;d^{n}y\;M(x,y,f(x),\partial _{i}f(x),f(y),\partial _{i}f(y))$

Die Funktionalableitung wird definiert, so dass sie sich verhält wie die Gradienten einer Funktion y() von mehreren Variablen:

$(d/dt)|_{t=0}(y(x+t\eta ))=\sum \nolimits _{i}\partial _{i}y\;\eta _{i}$

$(d/dt)|_{t=0}(A[f+t\eta ])=\sum \nolimits _{j}\int d^{n}x\;{\frac {\delta A}{\delta f_{j}(x)}}\eta _{j}(x)$

Das heißt, $\delta A/\delta f_{j}$ ist ein Integralkern für die lineare Approximation eines Funktionals, wenn zum Argument 'kleine' Abweichungen $\delta f(x)$ kommen.

Das Extremalproblem, etwa ein Prinzip der stationären Wirkung, führt auf ein Gleichungssystem ${\frac {\delta A}{\delta f_{j}(x)}}=0.$ Fürs lokale Funktional des ersten Beispiels ist die Funktionalableitung

${\frac {\delta A}{\delta f_{j}(x)}}={\frac {\partial L}{\partial f_{j}}}-\sum \nolimits _{i}{\frac {\partial }{\partial x_{i}}}{\frac {\partial L}{\partial (\partial _{i}f_{j})}}+\sum \nolimits _{i,k}{\frac {\partial ^{2}}{\partial x_{i}\;\partial x_{k}}}{\frac {\partial ^{2}L}{\partial (\partial _{i}\;\partial _{k}f_{j})}}$

Die alternierenden Vorzeichen kommen über partielle Integration herein.

Sei $A[f]=\int d^{n}x\;d^{n}y\;M(f(x),f(y),x,y).$

Bezeichne $M_{1i}:=\partial M(f,g,x,y)/\partial f_{i}(x);\;M_{2i}:=\partial M(f,g,x,y)/\partial g_{i}(y).$

${\frac {\delta A}{\delta f_{i}(x)}}=\int d^{n}y\;(M_{1i}(f(x),f(y),x,y)+M_{2i}(f(y),f(x),y,x))$

Ist M symmetrisch bei Vertauschung $(x,f(x))\leftrightarrow (y,f(y)),$ dann bekommt man als Integrand $2M_{1i}().$

Nebenbedingungen bei der Suche nach Extremalen von A[f] könen zweierlei Form haben. Entweder handelt es sich ebenfalls um Funktionale $B_{j}[f]=0,$ oder aber es sind wesentlich strengere punktweise Kriterien $C_{j}(f(x),\partial _{i}f(x))=0.$ Der erste Fall ist beispielsweise eine Normalisierung der Lösung f, der zweite Fall ist eine Eichung. Lagrange-Multiplikatoren, die sich an A[f] als weitere Argumente anfügen, sind einfache reelle Variablen bei Funktional- Bedingungen oder neue Funktions-Freiheitsgrade bei punktweisen Bedingungen.

Der Beweis, dass der Multiplikator-Algorithmus bei Funktionalen richtig ist, wäre zum Beispiel als Grenzübergang zu führen: ersetze die Definitions- Mannigfaltigkeit der Funktionen durch N Gitterpunkte. Dann sind die Funktionen einfach hochdimensionale Tupel; das Funktional wird beliebig gut angenähert durch eine Funktion endlich vieler Variablen. Nur bei endlicher Dimension greift das ursprüngliche Argument aus der Linearen Algebra zur Existenz der Multiplikatoren.
Intuitives Argument: Ein Funktional A[f] ist extremal, wenn jegliche Variation $f\to f+\delta f$ nur dann den Wert von A linear verändert, wenn $\delta f$ aus der erlaubten Teilmenge herausführt, also mindestens eines der $\delta B_{j}[f]$ nicht verschwindet. Das heißt, wenn die Funktionalvariation $\delta F[f]$ eine Linearkombination ist aus den Variationen der Nebenbedingungen.

Funktional mit Nebenbedingungen: ${\bar {A}}[f,y]=A[f]+\sum \nolimits _{j}y_{j}B_{j}[f].$

A[f] ist stationär für die Funktion $f=\{f_{i}\},$ wenn es ein Tupel $y=\{y_{j}\}$ so gibt, dass:

${\frac {\delta A}{\delta f_{i}(x)}}+\sum \nolimits _{j}y_{j}{\frac {\delta B_{j}}{\delta f_{i}(x)}}=0\;\forall i\;\forall x;\quad B_{j}[f]=0\;\forall j.$

Bei komplexwertigen Funktionen im Funktional greift ein wichtiger Trick aus der Differenzialrechnung: Statt algebraische Ausdrücke mühsam in den partiellen Ableitungen nach dem Realteil und Imaginärteil eines komplexen Arguments $z=u+iv$ zu entwickeln, benutzt man die partiellen Ableitungen nach dem Symbol $z$ und nach seinem Konjugierten $z^{*},$ zu betrachten als zwei unabhängige Variablen -- mit allen gewöhnlichen Produktregeln und so weiter. Es ergeben sich äquivalente Gleichungssysteme.
Ist ein Funktional $\phi \mapsto F[\phi ]$ reellwertig und invariant unter einer Symmetrie $(\phi _{i}({\vec {r}})\leftrightarrow \phi _{i}^{*}({\vec {r}})),$ dann reicht es aus für Extremalprinzipien und ähnliches, das Tupel der komplexwertigen Funktionalableitungen $\eta _{i}({\vec {r}})=\delta F/\delta \phi _{i}^{*}({\vec {r}})$ auszurechnen.

Beispiel 1: Funktional $F=\int d^{3}r\;\phi ^{*}({\vec {r}})\nabla ^{2}\phi ({\vec {r}})$

Einfach nach den Regeln: $\delta F/\delta \phi ^{*}({\vec {r}})=\nabla ^{2}\phi ({\vec {r}})$

Beispiel 2: $F=\int d^{3}r\;d^{3}s\;v({\vec {r}},{\vec {s}})\phi _{i}^{*}({\vec {r}})\phi _{j}({\vec {r}})\phi _{k}^{*}({\vec {s}})\phi _{l}({\vec {s}})$

$\delta F/\delta \phi _{h}^{*}({\vec {r}})=\int d^{3}s\;[\delta _{ih}\phi _{j}({\vec {r}})v({\vec {r}},{\vec {s}})\phi _{k}^{*}({\vec {s}})\phi _{l}({\vec {s}})+\delta _{kh}\phi _{l}({\vec {r}})v({\vec {s}},{\vec {r}})\phi _{i}^{*}({\vec {s}})\phi _{j}({\vec {s}})]$

Hartree-Fock-Gleichungen

Die Näherung von Hartree und Fock ist ein Variationsansatz, um das Minimum des Energie-Erwartungswertes für ein System aus n Elektronen zu suchen. Als Wellenfunktionen erlaubt man nur reine Produktfunktionen aus n linear unabhängigen Einteilchen-Wellen, genauer, die Slater-Determinanten. Verschränkte Funktionen werden also ignoriert! Nebenbedingung ist, dass die n Faktorwellen $\phi _{1},...,\phi _{n}$ ein Orthonormalsystem bilden. Der Erwartungswert des Hamiltonoperators ist dann ein Funktional auf der Menge der n-tupel von Funktionen $\{\phi _{i}({\vec {x}})\}.$ Er hat lokale kinetische und externe Anteile zweiter Ordnung, sowie einen bilokalen vierter Ordnung mit dem Wechselwirkungs-Potenzial. Unter Funktional-Nebenbedingungen wird das Extremalproblem mit Lagrange-Multiplikatoren aufgesetzt. Dies liefert n Integro-Differenzial-Gleichungen für die n Einteilchenwellen.
Die Anwendung auf Atome und Moleküle geschieht in der Approximation von unbeweglich angeordneten Kernen, die das externe Potenzial beisteuern.

Es geht darum, einen reellen Energie-Erwartungswert zu minimieren, mit folgenden Randbedingungen:

• Die n-Teilchen Wellenfunktion $\Psi ({\vec {r}}_{1},...{\vec {r}}_{n})$ ist eine Slater-Determinante von Einteilchen-Funktionen $\{\phi _{i}({\vec {r}});\;i=i...n\}.$
• $\int d^{3}r\;\phi _{j}^{*}({\vec {r}})\phi _{k}({\vec {r}})=\delta _{jk}$ (Orthonormierung).

Der Erwartungswert sowie die Normierungsbedingungen sind reelle Funktionale auf einer Menge von n-komponentigen Funktion $\phi =\{\phi _{1},...,\phi _{n}\}$ und sind Integrale über algebraische Terme der $\phi _{i}$ und ihrer komplex-Konjugierten. Daher brauchen wir nur die Funktionalableitungen $\delta F/\delta \phi _{i}^{*}({\vec {r}}).$

Auszurechnen ist der Energie-Erwartungswert der Slater-Determinanten mit orthonormalen $\{\phi _{i}\}.$ Die n-Elektronen-Welle hat die Form

$\Psi ({\vec {r}}_{1},...,{\vec {r}}_{n})=\det[\phi _{i}({\vec {r}}_{k})]/{\sqrt {n!}}$

Die Determinante kann entwickelt werden als Summe über die Permutationen der Indizes von $\phi$ oder derjenigen von ${\vec {r}};$ hier wird gebraucht:

${\sqrt {n!}}\;\Psi ({\vec {r}}_{1},...,{\vec {r}}_{n})=\sum \nolimits _{p\in Sn}{\text{sgn}}(p)\phi _{p1}({\vec {r}}_{1})\cdots \phi _{pn}({\vec {r}}_{n}).$

Mit dem externen Atomkerne-Potenzial $v({\vec {r}})$ und dimensionslosen Variablen ist der Hamilton-Operator H die Summe von:

• Einteilchen-Kinetischer Energie $T\Psi ={\Big [}-\sum \nolimits _{i}{\vec {\nabla }}_{i}^{2}/2{\Big ]}\Psi$
• Anziehung durch die Kerne $V\Psi ={\Big [}\sum \nolimits _{i}v({\vec {r}}_{i}){\Big ]}\Psi$
• Coulomb-Abstoßung $C\Psi ={\Big [}\sum \nolimits _{jk;\;j\neq k}(1/|{\vec {r}}_{j}-{\vec {r}}_{k}|){\Big ]}\Psi$

Die Coulomb-Summe hier zählt alle Paare von Indizes nur einmal, ungeordnet.
Nun zur Kombinatorik des Erwartungswerts $\langle H\rangle =\langle \Psi |H|\Psi \rangle .$

Ein Einteilchen-Operator $Q({\vec {r}}_{j})$ sieht jedes $\phi _{i},$ wenn es als $\phi _{i}(r_{j})$ in einem Term der Determinanten auf beiden Seiten auftaucht und wenn beide Terme die exakt gleiche Permutation sind (Orthonormalität). Das macht (n-1)! Kopien von $\langle \phi _{i}|Q|\phi _{i}\rangle ,$ alle mit positivem Vorzeichen. Insgesamt, weil die Fakultät mit zwei Faktoren $1/{\sqrt {n!}}$ von $\Psi$ wegschmilzt:

$\langle Q\rangle :=\langle \Psi |Q|\Psi \rangle =(1/n)\sum _{i}\langle \phi _{i}|Q|\phi _{i}\rangle .$

Eine symmetrische Summe von Einteilchen-Operatoren $R=\sum \nolimits _{j}Q(r_{j})$ ergibt $\langle \Psi |R|\Psi \rangle =\sum \nolimits _{i}\langle \phi _{i}|Q|\phi _{i}\rangle .$ Damit folgt der Anteil T+V des Hamilton-Operators

$\langle T+V\rangle =\langle \Psi |T+V|\Psi \rangle =\sum \nolimits _{i}\langle \phi _{i}|-{\vec {\nabla }}^{2}/2+v({\vec {r}})|\phi _{i}\rangle$

Eine symmetrischer Zweipunkt-Operator $P=P({\vec {r}}_{j},{\vec {r}}_{k});\;j\neq k,$ fängt zwischen zwei Determinanten jedes Paar von Wellenfaktoren $(\phi _{l},\phi _{m})$ viermal ein, mit Minuszeichen für die vertauschten Terme:

$\phi _{l}^{*}(r_{j})\phi _{l}(r_{j})P\phi _{m}^{*}(r_{k})\phi _{m}(r_{k})+\phi _{m}^{*}(r_{j})\phi _{m}(r_{j})P\phi _{l}^{*}(r_{k})\phi _{l}(r_{k})$

$-\phi _{l}^{*}(r_{j})\phi _{m}(r_{j})P\phi _{m}^{*}(r_{k})\phi _{l}(r_{k})-\phi _{m}^{*}(r_{j})\phi _{l}(r_{j})P\phi _{l}^{*}(r_{k})\phi _{m}(r_{k});$

weil jede Paarvertauschung eine Permutation mit Vorzeichen Minus ist. Für jeden Term fallen (n-2)! Kopien an für die Permutationen der anderen Indizes außer j,k. Nur gleich permutierte Terme beider Slater-Determinanten lassen dabei wegen Orthonormalität den Faktor 1 durch. Alles zusammengefasst mit Integralvariablen r,s und einem Faktor 2*(n-2)!:

$\langle \Psi |P|\Psi \rangle =2/n/(n-1)\sum _{lm;\;l\neq m}\int d^{3}r\;d^{3}s\;[\phi _{l}^{*}(r)\phi _{l}(r)P(r,s)\phi _{m}^{*}(s)\phi _{m}(s)-\phi _{l}^{*}(r)\phi _{m}(r)P(r,s)\phi _{m}^{*}(s)\phi _{l}(s)]$

Das Ergebnis hängt nicht von den Indizes (j,k) ab. Die Summe geht über die Menge der (ungeordneten) Paare und ist nicht als Doppelsumme über zwei Laufindizes zu lesen. Vertauschte Ordnung eines Paares $(l\leftrightarrow m)$ verändert keines der zwei Integrale, weil jeweils zu kompensieren mit $r\leftrightarrow s.$ Die Symmetrie P(r,s)=P(s,r) wird ja vorausgesetzt.

Der Coulomb-Operator ist nun die Summe von n(n-1)/2 genau solchen Zweipunkt-Termen. Alle erzeugen sie Gleiches, also gilt da der letzte Ausdruck. Der Faktor vor der Summe hebt sich weg, und $P(r,s)=1/(|{\vec {r}}-{\vec {s}}|).$ Schreibweise des Coulomb-Funktionals mit Blick auf baldige Funktionalableitung:

$\langle C\rangle =\langle \Psi |C|\Psi \rangle =\sum \nolimits _{lm;\;l\neq m}\int d^{3}r\;\phi _{l}^{*}({\vec {r}}){\Big [}\int d^{3}s\;{\big (}\phi _{m}^{*}({\vec {s}})\phi _{m}({\vec {s}})\cdot \phi _{l}({\vec {r}})-(\phi _{m}^{*}(s)\phi _{m}(r))\cdot \phi _{l}(s){\big )}/|{\vec {r}}-{\vec {s}}|{\Big ]},$

$\langle C\rangle =:\langle D-X\rangle .$

Hier unterscheiden wir den direkten Anteil D und den Austausch-Teil X, der all die Terme mit Minuszeichen aufsaugt.

Zusammengefasst. Der Hamilton-Operator-Erwartungswert ist also ein Funktional auf einem Raum der Wellen-Tupel

$\phi =\{\phi _{1},...\phi _{n}\}\;\mapsto \;E[\phi ]=\langle \Psi |T+V+C|\Psi \rangle$

mit der Bedingung, dass Wellentupel orthonormal sind. Gesucht wird ein Minimum des Funktionals.

Das Energie-Funktional $E[\phi ]$ wird nun mit Lagrange-Faktoren ergänzt zu

$F[\phi ]:=E[\phi ]-\sum \nolimits _{j,k}\epsilon _{jk}\langle \phi _{j}|\phi _{k}\rangle .$

$F=F^{*}$ muss reell sein wie E, bei unbeschränkter Variation von $\phi .$ Daraus folgt eine hermitesche Multiplikator-Matrix $\epsilon _{jk}=\epsilon _{kj}^{*}.$ Außerdem ist F invariant unter den unitären Matrixtransformationen:

${\phi '}_{i}=\sum \nolimits _{j}U_{ij}\phi _{j};\;{\phi '}_{i}^{*}=\sum \nolimits _{j}U_{ij}^{*}\phi _{j}^{*};\;{\epsilon '}_{jk}=\sum \nolimits _{hl}U_{jh}\epsilon _{hl}U_{lk}^{\dagger };\;UU^{\dagger }=U^{\dagger }U=\mathbf {1}$

Denn die Summen $\sum \nolimits _{j,k}\epsilon _{jk}\phi _{j}^{*}(r)\phi _{k}(r)$ bleiben gleich und die Slater-Determinanten

$\det S\propto \det\{S_{ik}\}=\det[\phi _{i}({\vec {r}}_{k})]$

transformieren sich zu $\det[\sum \nolimits _{j}U_{ij}\phi _{j}({\vec {r}}_{k})]=(\det U)(\det S).$ Wegen Unitarität ist $(\det U)$ ein konstanter Faktor vom Betrag 1, der die Erwartungswerte von hermiteschen Operatoren nicht ändert.
Es kann nun jede hermitesche Matrix $\epsilon _{jk}$ unitär auf Diagonalform gebracht werden. Nach einer solchen Symmetrietransformation hat das Funktional F eine Darstellung mit einem physikalisch gleichen Zustand, der wesentlich einfacher aussieht:

$F[\phi ]=E[\phi ]-\sum \nolimits _{j}\epsilon _{j}\langle \phi _{j}|\phi _{j}\rangle ;\quad (\langle \phi _{j}|\phi _{j}\rangle =\int d^{3}r\;\phi _{j}^{*}({\vec {r}})\phi _{j}({\vec {r}}))$

Wird so die unitäre Redundanz rausgeworfen, reduziert sich das Variationsproblem damit auf n reelle Lagrange-Faktoren $\epsilon _{i}$ und n Gleichungen für die Einteilchenfunktionen:

$\delta F/\delta \phi _{i}^{*}=0\;\Rightarrow \;\delta E/\delta \phi _{i}^{*}({\vec {r}})=\epsilon _{i}\phi _{i}({\vec {r}}).$

Die Funktionalableitungen von $E[\phi ]=\langle T+V+D-X\rangle$ wurden schon ziemlich vorgekaut, weil oben bewusst die Teile von E als Integrale von $\phi _{i}^{*}$ mal irgendwas vorkommen. Man liest direkt die Anteile (T+V) und Q=(D-X) ab.

$\delta (T+V)/\delta \phi _{i}^{*}(r)=(-{\vec {\nabla }}^{2}/2+v({\vec {r}}))\phi _{i}({\vec {r}})$

Beim Q-Anteil wird wegen der besprochenen Paarsymmetrie aus der Summe $\sum \nolimits _{lm}$ für die Ableitung Nummer i nur l=i ausgewählt, der andere Index durchläuft alle Werte ungleich i. Daher mit $P({\vec {r}},{\vec {s}})=1/|{\vec {r}}-{\vec {s}}|$

$\delta Q/\delta \phi _{i}^{*}(r)=\int d^{3}s\sum \nolimits _{m;\;m\neq i}\phi _{m}^{*}({\vec {s}})\phi _{m}({\vec {s}})P({\vec {r}},{\vec {s}})\phi _{i}({\vec {r}})$

$-\int d^{3}s\sum \nolimits _{m;\;m\neq i}\phi _{m}^{*}({\vec {s}})\phi _{m}({\vec {r}})P({\vec {r}},{\vec {s}})\phi _{i}({\vec {s}})$

Noch ein Trick: Nun wird formal auch m=i in den Summen zugelassen, weil sich diese Terme aufheben. Die zwei Teile des effektiven Potenzials saugen jeweils einen Selbstenergie-Term auf.

Die Dichte der Elektronenwolke ist folgende Summe von Betragsquadraten

$\rho ({\vec {r}})=\sum \nolimits _{m}\phi _{m}^{*}({\vec {r}})\phi _{m}({\vec {r}})=\sum \nolimits _{m}|\phi _{m}({\vec {r}})|^{2}.$

Es gilt $\int d^{3}r\;\rho ({\vec {r}})=n.$

Eine Zweipunkt-Korrelation sei definiert als

$\zeta ({\vec {r}},{\vec {s}})=\sum \nolimits _{m}\phi _{m}^{*}({\vec {s}})\phi _{m}({\vec {r}});\quad (\zeta ({\vec {s}},{\vec {r}})=\zeta ^{*}({\vec {r}},{\vec {s}})).$

Dann liest sich die Ableitung von Q als ein klassisches lokales Potenzial mit der Ladungsdichte $\rho$ minus ein nichtlokales Potenzial mit der Korrelation.

$\delta Q/\delta \phi _{i}^{*}(r)=u_{1}({\vec {r}})\phi _{i}({\vec {r}})-\int d^{3}s\;u_{2}({\vec {r}},{\vec {s}})\phi _{i}({\vec {s}})$

$u_{1}({\vec {r}})=\int d^{3}s\;\rho ({\vec {s}})/|{\vec {r}}-{\vec {s}}|$

$u_{2}({\vec {r}},{\vec {s}})=\zeta ({\vec {r}},{\vec {s}})/|{\vec {r}}-{\vec {s}}|$

Der Teil $u_{1}$ als Coulomb-Potenzial überrascht nicht. Der andere Teil heißt Austausch-Wechselwirkung -- ein Quanteneffekt für gleiche Fermionen.

Der Erwartungswert $\langle C\rangle$ kann als Funktion von $\rho$ ausgedrückt werden, wenn man Anteile l=m dazu zählt und mit der halben Doppelsumme hantiert.

$\langle C\rangle =\langle \Psi |C|\Psi \rangle =\sum \nolimits _{lm;\;l\neq m}\int d^{3}r\;\int d^{3}s\;{\Big [}\phi _{l}^{*}({\vec {r}})\phi _{l}({\vec {r}})\phi _{m}^{*}({\vec {s}})\phi _{m}({\vec {s}})-\phi _{l}^{*}({\vec {r}})(\phi _{m}^{*}(s)\phi _{m}(r)\phi _{l}(s){\Big ]}/|{\vec {r}}-{\vec {s}}|.$

Weil für l=m die Termpaare gleich sind, funktioniert die Transformation

$\sum \nolimits _{lm;\;l\neq m}[*]\;\to \;{\tfrac {1}{2}}\sum \nolimits _{l}\sum \nolimits _{m}[*]\;\Rightarrow \;\langle C\rangle =\langle D'-X'\rangle$

$\langle D'\rangle ={\frac {1}{2}}\sum \nolimits _{l}\sum \nolimits _{m}\int d^{3}r\;\int d^{3}s\;(1/|{\vec {r}}-{\vec {s}}|)\phi _{l}^{*}({\vec {r}})\phi _{l}({\vec {r}})\phi _{m}^{*}({\vec {s}})\phi _{m}({\vec {s}})$

$={\frac {1}{2}}\int d^{3}r\;d^{3}s\;\rho ({\vec {r}})\rho ({\vec {s}})/|{\vec {r}}-{\vec {s}}|.$

$\langle X'\rangle ={\frac {1}{2}}\sum \nolimits _{l}\sum \nolimits _{m}\int d^{3}r\;\int d^{3}s\;(1/|{\vec {r}}-{\vec {s}}|)\phi _{l}^{*}({\vec {r}})\phi _{m}({\vec {s}})\phi _{m}^{*}({\vec {s}})\phi _{l}({\vec {r}})$

$={\frac {1}{2}}\int d^{3}r\;d^{3}s\;|\zeta ({\vec {r}},{\vec {s}})|^{2}/|{\vec {r}}-{\vec {s}}|.$

Die Teilsumme $\langle D'\rangle$ ist die klassische Coulomb-Energie einer Ladungswolke, die andere Summe $\langle X'\rangle$ zieht davon einen manifest positiven Anteil ab. Wurde aber nicht gesagt, dass gleiche Fermionen sich mehr abstoßen sollen als unterschiedliche Teilchen? Tun sie. Diese Abstoßung ist hier über das Pauli-Prinzip in den gezwungen orthogonalen Einteilchenwellen verpackt; der Austauschterm korrigiert offenbar etwas von dessen Übertreibung.

Das n-Tupel von Variationsgleichungen sind gekoppelte Gleichungen vom zeitunabhängigen Schrödinger-Typ, in denen die Wellen nichtlinear auf das Potenzial rückwirken -- darüber hinaus taucht ein nichtlokaler Integralkern auf. Der effektive Einteilchen-Hamilton-Operator hängt von den Lösungen ab.

${\bar {H}}\phi ({\vec {r}})=(-{\vec {\nabla }}^{2}/2+v+u_{1})\phi ({\vec {r}})-\int d^{3}s\;u_{2}({\vec {r}},{\vec {s}})\phi ({\vec {s}})$

Damit sieht der Hartree-Fock-Algorithmus für den Grundzustand etwa so aus: Zuerst errate man zu vorgegebenen Kern-Positionen eines Moleküls eine vorläufige Ladungsdichte $\rho$ und versuche sie anzupassen mit einem n-Tupel von orthonormalen Wellenfunktionen. Eine Möglichkeit, den Motor der Iteration anzukurbeln, ist es, die vereinfachten Hartree-Gleichungen zu lösen, worin der Austauschteil entfällt und Schrödingergleichungen dableiben:

$(-{\vec {\nabla }}^{2}/2+v({\vec {r}})+u_{1}({\vec {r}}))\phi _{i}({\vec {r}})=\epsilon _{i}\phi _{i}({\vec {r}});\quad u_{1}({\vec {r}})=\int d^{3}s\;\rho ({\vec {s}})/|{\vec {r}}-{\vec {s}}|$

Selbstkonsistenz-Bedingung: $\rho ({\vec {r}})=\sum \nolimits _{i}|\phi _{i}({\vec {r}})|^{2}.$

Dann berechne man mit den n niedrigsten Wellenfunktionen die Ladungsdichte neu und die Zweipunkt-Korrelation. Danach finde man n verbesserte Zustände der modifizierten Einteilchen-Schrödinger-Gleichung mit nichtlokalem Potenzial. Mit dem neuen n-Tupel wiederhole man die Prozedur Wellenfunktionen $\leftrightarrow$ Ladungsdichten, bis hoffentlich Konvergenz erfolgt. Die selbstkonsistene Iteration ist natürlich sehr rechenintensiv und garantiert keine große Präzision. Hartree-Fock ist eine Referenzmethode, auf der zahlreiche verfeinerte Algorithmen für Computereinsatz aufbauen.

Zu beachten. Die Summe der $\epsilon _{i}$ ist nicht die Energie des Grundzustands.

$\sum \nolimits _{i}\epsilon _{i}=\sum \nolimits _{i}\langle \phi _{i}|{\bar {H}}|\phi _{i}\rangle =\langle T+V\rangle +\int d^{3}r\;d^{3}s\;{\big [}\rho ({\vec {r}})\rho ({\vec {s}})-|\zeta (r,s)|^{2}{\big ]}/|{\vec {r}}-{\vec {s}}|$

$=\langle T+V+2(D'-X')\rangle =\langle E\rangle +\langle D'-X'\rangle$

Hier kommen Doppelsummen in den Elektronpaar-Wechselwirkungen vor. Sie werden in der Epsilon-Summe doppelt gezählt im Vergleich zu den einfachen Paarsummen im Energie-Funktional.
Zumindest messen die $\epsilon _{i}$ nicht ganz unrealistisch die Ionisierungs- Energie für Elektronen. Die Energie, die man aufbringen muss, um das eine oder andere Elektron vom Molekül abzutrennen.

Der Spin wurde bis hierhin nicht erwähnt. Er ändert nicht viel an der Theorie. Statt $\phi _{i}({\vec {r}})$ gibt es $\phi _{i}(\sigma ,{\vec {r}})$ mit einem diskreten Index und jedes Integral wird zu $\sum \nolimits _{\sigma }\int d^{3}r.$ Die Coulomb-Wechselwirkung hat in nullter Näherung keine Spin-Abhängigkeit. Slater-Determinanten sind nicht Null, wenn zum Beispiel zwei $\phi _{i}$ die gleiche Ortsfunktion haben, jedoch linear unabhängige Tupel (mit $\sigma$ indexiert) als Faktoren. Jedes Orbital-Niveau kann von zwei Elektronen mit verschiedenen Spinrichtungen besetzt werden.