D.54 De­riva­tion of the Hartree-Fock equa­tions

This note de­rives the canon­i­cal Hartree-Fock equa­tions. The de­riva­tion be­low will be per­formed un­der the nor­mally stated rules of en­gage­ment that the or­bitals are of the form $\pe{n}//u//$ or $\pe{n}//d//$, so that only the spa­tial or­bitals $\pe{n}////$ are con­tin­u­ously vari­able. The de­riva­tions al­low for the fact that some spa­tial spin states may be con­strained to be equal.

First, you can make things a lot less messy by a pri­ori spec­i­fy­ing the or­der­ing of the or­bitals. The or­der­ing makes no dif­fer­ence phys­i­cally, and it sim­pli­fies the math­e­mat­ics. In par­tic­u­lar, in re­stricted Hartree-Fock some spa­tial or­bitals ap­pear in pairs, but you can only count each spa­tial or­bital as one un­known func­tion. The eas­i­est way to han­dle that is to push the spin-down ver­sions of the du­pli­cated spa­tial or­bitals to the end of the Slater de­ter­mi­nant. Then the start of the de­ter­mi­nant com­prises a list of unique spa­tial or­bitals.

So, it will be as­sumed that the or­bitals are or­dered as fol­lows:

1.
the paired spa­tial states in their spin-up ver­sion; as­sume there are $N_p$ $\raisebox{-.5pt}{$\geqslant$}$ 0 of them;
2.
un­paired spin-up states; as­sume there are $N_u$ of them;
3.
un­paired spin-down states; as­sume there are $N_d$ of them;
4.
and fi­nally, the paired spa­tial states in their spin-down ver­sion.
That means that the Slater-de­ter­mi­nant wave func­tion looks like:

\begin{displaymath}
{\left\vert{\rm det}\;
\pe 1//u//,\ldots,\pe N_p//u//,\pe ...
...\ldots,\pe N//d//,\pe N+1//d//,\ldots,\pe I//d//\right\rangle}
\end{displaymath}


\begin{displaymath}
N_1=N_p+N_u \qquad N=N_p+N_u+N_d
\end{displaymath}


\begin{displaymath}
\pe N+1////=\pe 1////,\quad \pe N+2////=\pe 2////,\quad \ldots,\quad
\pe I////=\pe N_p////
\end{displaymath}

The to­tal num­ber of un­known spa­tial or­bitals is $N$, and you need a cor­re­spond­ing $N$ equa­tions for them.

The vari­a­tional method dis­cussed in chap­ter 9.1 says that the ex­pec­ta­tion en­ergy must be un­changed un­der small changes in the or­bitals, pro­vided that penalty terms are added for changes that vi­o­late the or­tho­nor­mal­ity re­quire­ments on the or­bitals.

The ex­pec­ta­tion value of en­ergy was in chap­ter 9.3.3 found to be:

\begin{eqnarray*}
\left\langle{E}\right\rangle & = &\sum_{n=1}^I \langle\pe n//...
...angle{\updownarrow}_n\vert{\updownarrow}_{\underline n}\rangle^2
\end{eqnarray*}

(From here on, the ar­gu­ment of the first or­bital of a pair in ei­ther side of an in­ner prod­uct is taken to be the first in­ner prod­uct in­te­gra­tion vari­able ${\skew0\vec r}$ and the ar­gu­ment of the sec­ond or­bital is the sec­ond in­te­gra­tion vari­able ${\underline{\skew0\vec r}}$)

The penalty terms re­quire penalty fac­tors called La­grangian vari­ables. The penalty fac­tor for vi­o­la­tions of the nor­mal­iza­tion re­quire­ment

\begin{displaymath}
\langle\pe n////\vert\pe n////\rangle=1
\end{displaymath}

will be called $\epsilon_{nn}$ for rea­sons ev­i­dent in a sec­ond. Or­thog­o­nal­ity be­tween any two spa­tial or­bitals $\pe{n}////$ and $\pe{\underline n}////$ re­quires

\begin{displaymath}
{\textstyle\frac{1}{2}}\Big(\langle\pe n////\vert\pe{\under...
...e -
\langle\pe{\underline n}////\vert\pe n////\rangle\Big)=0.
\end{displaymath}

where the first con­straint says that the real part of $\langle\pe{n}////\vert\pe{\underline n}////\rangle$ must be zero and the sec­ond that its imag­i­nary part must be zero too. (Re­mem­ber that if you switch sides in an in­ner prod­uct, you turn it into its com­plex con­ju­gate.) To avoid in­clud­ing the same or­thog­o­nal­ity con­di­tion twice, the con­straint will be writ­ten only for ${\underline n}$ $\raisebox{.3pt}{$>$}$ $n$. The penalty fac­tor for the first con­straint will be called $2\epsilon_{n{\underline n},r}$, and the one for the sec­ond con­straint $2\epsilon_{n{\underline n},i}$.

In those terms, the pe­nal­ized change in ex­pec­ta­tion en­ergy be­comes, in the re­stricted case that all unique spa­tial or­bitals are mu­tu­ally or­thog­o­nal,

\begin{displaymath}
\delta \left\langle{E}\right\rangle -
\sum_{n=1}^N \sum_{{...
...} \delta \langle\pe n////\vert\pe{\underline n}////\rangle = 0
\end{displaymath}

where $\epsilon_{n{\underline n}}$ is an Her­mit­ian ma­trix, with $\epsilon_{n{\underline n}}$ $\vphantom0\raisebox{1.5pt}{$=$}$ $\epsilon_{n{\underline n},r}+{\rm i}\epsilon_{n{\underline n},i}$. The no­ta­tions for the La­grangian vari­ables were cho­sen above to achieve this fi­nal re­sult.

But for un­re­stricted Hartree-Fock, spa­tial or­bitals are not re­quired to be or­thog­o­nal if they have op­po­site spin, be­cause the spins will take care of or­thog­o­nal­ity. You can re­move the er­ro­neously added con­straints by sim­ply spec­i­fy­ing that the cor­re­spond­ing La­grangian vari­ables are zero:

\begin{displaymath}
\epsilon_{n{\underline n}}=0\mbox{ if unrestricted Hartree-...
...e{\updownarrow}_n\vert{\updownarrow}_{\underline n}\rangle=0$}
\end{displaymath}

or equiv­a­lently, if $n$ $\raisebox{-.3pt}{$\leqslant$}$ $N_u$, ${\underline n}$ $\raisebox{.3pt}{$>$}$ $N_u$ or $n$ $\raisebox{.3pt}{$>$}$ $N_u$, ${\underline n}$ $\raisebox{-.3pt}{$\leqslant$}$ $N_u$.

Now work out the pe­nal­ized change in ex­pec­ta­tion en­ergy due to a change in the val­ues of a se­lected spa­tial or­bital $\pe{m}////$ with $m$ $\raisebox{-.3pt}{$\leqslant$}$ $N$. It is

\begin{eqnarray*}
&&
\sum_{\pe n////=\pe m////}
\Big(
\langle\delta\pe m////...
...N \epsilon_{nm} \langle\pe n////\vert\delta\pe m////\rangle
= 0
\end{eqnarray*}

OK, OK it is a mess. Sums like $\pe{n}////$ $\vphantom0\raisebox{1.5pt}{$=$}$ $\pe{m}////$ are for the re­stricted Hartree-Fock case, in which spa­tial or­bital $\pe{m}////$ may ap­pear twice in the Slater de­ter­mi­nant. From now on, just write them as $[2]$, mean­ing, put in a fac­tor 2 if or­bital $\pe{m}////$ ap­pears twice. The ex­cep­tion is for the ex­change in­te­grals which pro­duce ex­actly one nonzero spin prod­uct; write that as $[\langle{\updownarrow}_n\vert{\updownarrow}_m\rangle^2]$, mean­ing take out that prod­uct if the or­bital ap­pears twice.

Next, note that the sec­ond term in each row is just the com­plex con­ju­gate of the first. Con­sid­er­ing ${\rm i}\delta\pe{m}////$ as a sec­ond pos­si­ble change in or­bital, as was done in the ex­am­ple in chap­ter 9.1, it is seen that the first terms by them­selves must be zero, so you can just ig­nore the sec­ond term in each row. And the in­te­grals with the fac­tors ${\textstyle\frac{1}{2}}$ are pair­wise the same; the dif­fer­ence is just a name swap of the first and sec­ond sum­ma­tion and in­te­gra­tion vari­ables. So all that you re­ally have left is

\begin{eqnarray*}
&& [2] \langle\delta\pe m////\vert h^{\rm e}\vert\pe m////\ra...
...{n=1}^N \epsilon_{mn}\langle\delta\pe m////\vert\pe n////\rangle
\end{eqnarray*}

Now note that if you write out the in­ner prod­uct over the first po­si­tion co­or­di­nate, you will get an in­te­gral of the gen­eral form

\begin{eqnarray*}
\lefteqn{\int_{{\rm all}\;{\skew0\vec r}}\delta\pe m////\stru...
..._{n=1}^N\epsilon_{mn}\pe n////
\bigg) {\,\rm d}^3{\skew0\vec r}
\end{eqnarray*}

If this in­te­gral is to be zero for what­ever you take $\delta\pe{m}////$, then the terms within the paren­the­ses must be zero. (Just take $\delta\pe{m}////$ pro­por­tional to the par­en­thet­i­cal ex­pres­sion; you would get the in­te­gral of an ab­solute square, only zero if the square is.) Un­avoid­ably, you must have that

\begin{displaymath}[2]h^{\rm e}\pe m////
+ \sum_{n=1}^I [2] \langle\pe n////\ve...
... m////\rangle\pe n////
=
\sum_{n=1}^N \epsilon_{mn}\pe n////
\end{displaymath}

You can di­vide by [2]:

\begin{displaymath}
\fbox{$\displaystyle
h^{\rm e}\pe m////
+ \sum_{n=1}^I \l...
...end{array} \hspace{-4pt}
\right\}
\epsilon_{mn}\pe n////
$}
\end{displaymath} (D.35)

where you use the lower of the choices be­tween the braces in the case that spa­tial or­bital $\pe{m}////$ ap­pears twice in the Slater de­ter­mi­nant, or equiv­a­lently, if $m$ $\raisebox{-.3pt}{$\leqslant$}$ $N_p$. There you have the Hartree-Fock equa­tions, one for each $m$ $\raisebox{-.3pt}{$\leqslant$}$ $N$. Re­call that they ap­ply as­sum­ing the or­der­ing of the Slater de­ter­mi­nant given in the be­gin­ning, and that for un­re­stricted Hartree-Fock, $\epsilon_{mn}$ is zero if $\langle{\updownarrow}_m\vert{\updownarrow}_n\rangle$ $\vphantom0\raisebox{1.5pt}{$=$}$ 0 is.

How about those $\epsilon_{mn}$, you say? Shouldn’t the right hand side just be $\epsilon_m\pe{m}////$? Ah, you want the canon­i­cal Hartree-Fock equa­tions, not just the plain vanilla ver­sion.

OK, let’s do the re­stricted closed-shell Hartree-Fock case first, then, since it is the eas­i­est one. Every state is paired, so the lower choice in the curly brack­ets al­ways ap­plies, and the num­ber of unique un­known spa­tial states is $N$ $\vphantom0\raisebox{1.5pt}{$=$}$ $I$$\raisebox{.5pt}{$/$}$​2. Also, you can re­duce the sum­ma­tion up­per lim­its $n$ $\vphantom0\raisebox{1.5pt}{$=$}$ $I$ to $n$ $\vphantom0\raisebox{1.5pt}{$=$}$ $N$ $\vphantom0\raisebox{1.5pt}{$=$}$ $I$$\raisebox{.5pt}{$/$}$​2 if you add a fac­tor 2, since the sec­ond half of the spa­tial or­bitals are the same as the first half. So you get

\begin{displaymath}
h^{\rm e}\pe m////
+ 2 \sum_{n=1}^N \langle\pe n////\vert ...
... =
\sum_{n=1}^N {\textstyle\frac{1}{2}}\epsilon_{mn}\pe n////
\end{displaymath}

Now, sup­pose that you de­fine a new set of or­bitals, each a lin­ear com­bi­na­tion of the cur­rent ones:

\begin{displaymath}
\overline\pe\nu//// \equiv \sum_{n=1}^N u_{\nu n} \pe n////
\qquad\mbox{for $\nu=1,2,\ldots,N$}
\end{displaymath}

where the $u_{{\nu}n}$ are the mul­ti­ples of the orig­i­nal or­bitals. Will the new ones still be an or­tho­nor­mal set? Well, they will be if

\begin{displaymath}
\langle\overline\pe\mu////\vert\overline\pe\nu////\rangle=\delta_{\mu\nu}
\end{displaymath}

where $\delta_{\mu\nu}$ is the Kro­necker delta, one if $\mu$ $\vphantom0\raisebox{1.5pt}{$=$}$ $\nu$, zero oth­er­wise. Sub­sti­tut­ing in the de­f­i­n­i­tion of the new or­bitals, mak­ing sure not to use the same name for two dif­fer­ent in­dices,

\begin{displaymath}
\sum_{m=1}^N \sum_{n=1}^N u^*_{\mu m} u_{\nu n}
\langle\pe m////\vert\pe n////\rangle=\delta_{\mu\nu}.
\end{displaymath}

Now note that the $\pe{n}////$ are or­tho­nor­mal, so to get a nonzero value, $m$ must be $n$, and you get

\begin{displaymath}
\sum_{n=1}^N u^*_{\mu n} u_{\nu n} = \delta_{\mu\nu}.
\end{displaymath}

Con­sider $n$ to be the com­po­nent in­dex of a vec­tor. Then this re­ally says that vec­tors $\vec{u}_\mu$ and $\vec{u}_\nu$ must be or­tho­nor­mal. So the ma­trix of co­ef­fi­cients must con­sist of or­tho­nor­mal vec­tors. Math­e­mati­cians call such ma­tri­ces “uni­tary,” rather than or­tho­nor­mal, since it is eas­ily con­fused with unit, and that keeps math­e­mati­cians in busi­ness ex­plain­ing all the con­fu­sion.

Call the com­plete ma­trix $U$. Then, ac­cord­ing to the rules of ma­trix mul­ti­pli­ca­tion and Her­mit­ian ad­joint, the or­tho­nor­mal­ity con­di­tion above is equiv­a­lent to $UU^{\rm {H}}$ $\vphantom0\raisebox{1.5pt}{$=$}$ $I$ where $I$ is the unit ma­trix. That means $U^{\rm {H}}$ is the in­verse ma­trix to $U$, $U^{\rm {H}}$ $\vphantom0\raisebox{1.5pt}{$=$}$ $U^{-1}$ and then you also have $U^{\rm {H}}U$ $\vphantom0\raisebox{1.5pt}{$=$}$ $I$:

\begin{displaymath}
\sum_{\nu=1}^N u^*_{\nu m} u_{\nu n} = \delta_{mn}.
\end{displaymath}

Now pre­mul­ti­ply the de­f­i­n­i­tion of the new or­bitals above by $U^{\rm {H}}$; you get

\begin{displaymath}
\sum_{\nu=1}^N u^*_{\nu m} \overline\pe\nu//// =
\sum_{\nu=1}^N \sum_{n=1}^N u^*_{\nu m} u_{\nu n} \pe n////
\end{displaymath}

but the sum over $\nu$ in the right hand side is just $\delta_{mn}$ and you get

\begin{displaymath}
\sum_{\nu=1}^N u^*_{\nu m} \overline\pe\nu//// =
\sum_{n=1}^N \delta_{mn} \pe n//// = \pe m////.
\end{displaymath}

That gives you an ex­pres­sion for the orig­i­nal or­bitals in terms of the new ones. For aes­thetic rea­sons, you might just as well reno­tate $\nu$ to $\mu$, the Greek equiv­a­lent of $m$, to get

\begin{displaymath}
\pe m//// = \sum_{\mu=1}^N u^*_{\mu m} \overline\pe\mu////.
\end{displaymath}

Now plug that into the non­canon­i­cal re­stricted closed-shell Hartree-Fock equa­tions, with equiv­a­lent ex­pres­sions for $\pe{n}////$ us­ing what­ever sum­ma­tion vari­able is still avail­able,

\begin{displaymath}
\pe n//// = \sum_{\mu=1}^N u^*_{\mu n} \overline\pe\mu////
...
.../ = \sum_{\lambda=1}^N u^*_{\lambda n} \overline\pe\lambda////
\end{displaymath}

and use the re­duc­tion for­mula $UU^{\rm {H}}$ $\vphantom0\raisebox{1.5pt}{$=$}$ $I$,

\begin{displaymath}
\sum_{m=1}^N u_{\nu m} u^*_{\mu m} = \delta_{\mu\nu}
\quad...
..._{n=1}^N u_{\kappa n} u^*_{\lambda n} = \delta_{\kappa\lambda}
\end{displaymath}

pre­mul­ti­ply­ing all by $U$, i.e. put $\sum_{m=1}^{N}u_{\nu{m}}$ be­fore each term. You get

\begin{displaymath}
h^{\rm e}\overline\pe\nu////
+ 2 \sum_{\lambda=1}^N
\lang...
...}{2}}
u_{\nu m} \epsilon_{mn} u^*_{\mu n} \overline\pe\mu////
\end{displaymath}

Note that the only thing that has changed more than just by sym­bol names is the ma­trix in the right hand side. Now for each value of $\mu$, take $u^*_{{\mu}n}$ as the $\mu$-​th or­tho­nor­mal eigen­vec­tor of Her­mit­ian ma­trix $\epsilon_{mn}$, call­ing the eigen­value $2\epsilon_\mu$. Then the right hand side be­comes

\begin{displaymath}
\sum_{m=1}^{I} \sum_{\mu=1}^N
u_{\nu m} \epsilon_{\mu} u^*...
...mu} \overline\pe\mu////
=
\epsilon_{\nu} \overline\pe\nu////
\end{displaymath}

So, in terms of the new or­bitals de­fined by the re­quire­ment that $u^*_{{\mu}n}$ gives the eigen­vec­tors of $\epsilon_{mn}$, the right hand side sim­pli­fies to the canon­i­cal one.

Since you no longer care about the old or­bitals, you can drop the over­lines on the new ones, and re­vert to sen­si­ble ro­man in­dices $n$ and ${\underline n}$ in­stead of the Greek ones $\nu$ and $\lambda$. You then have the canon­i­cal re­stricted closed-shell Hartree-Fock equa­tions

\begin{displaymath}
\fbox{$\displaystyle
h^{\rm e}\pe n////
+ 2 \sum_{{\under...
...n////\rangle\pe{\underline n}////
=
\epsilon_n \pe n////
$}
\end{displaymath} (D.36)

that, as noted, as­sume that the Slater de­ter­mi­nant is or­dered so that the $I$$\raisebox{.5pt}{$/$}$​2 spin-up or­bitals are at the start of it. Note that the left hand side di­rectly pro­vides a Her­mit­ian Fock op­er­a­tor if you iden­tify it as ${\cal F}\pe{n}////$; there is no need to in­volve spin here.

In the un­re­stricted case, the non­canon­i­cal equa­tions are

\begin{displaymath}
h^{\rm e}\pe m////
+ \sum_{n=1}^I \langle\pe n////\vert v^...
... m////\rangle\pe n////
=
\sum_{n=1}^I \epsilon_{mn}\pe n////
\end{displaymath}

In this case the spin-up and spin-down spa­tial states are not mu­tu­ally or­tho­nor­mal, and you want to re­de­fine the group of spin up states and the group of spin down states sep­a­rately.

The term in lin­ear al­ge­bra is that you want to par­ti­tion your $U$ ma­trix. What that means is sim­ply that you sep­a­rate the or­bital num­bers into two sets. The set of num­bers $n$ $\raisebox{-.3pt}{$\leqslant$}$ $N_u$ of spin-up or­bitals will be in­di­cated as U, and the set of val­ues $n$ $\raisebox{.3pt}{$>$}$ $N_u$ of spin-down ones by D. So you can par­ti­tion (sep­a­rate) the non­canon­i­cal equa­tions above into equa­tions for $m\in{\rm {U}}$ (mean­ing $m$ is one of the val­ues in set U):

\begin{displaymath}
h^{\rm e}\pe m////
+ \sum_{n\in{\rm U}} \langle\pe n////\v...
...\rangle\pe n////
=
\sum_{n\in{\rm U}} \epsilon_{mn}\pe n////
\end{displaymath}

and equa­tions for $m\in{\rm {D}}$

\begin{displaymath}
h^{\rm e}\pe m////
+ \sum_{n\in{\rm U}} \langle\pe n////\v...
...\rangle\pe n////
=
\sum_{n\in{\rm D}} \epsilon_{mn}\pe n////
\end{displaymath}

In these two equa­tions, the fact that the up and down spin states are or­thog­o­nal was used to get rid of one pair of sums, and an­other pair was elim­i­nated by the fact that there are no La­grangian vari­ables $\epsilon_{mn}$ link­ing the sets, since the spa­tial or­bitals in the sets are al­lowed to be mu­tu­ally nonorthog­o­nal.

Now sep­a­rately re­place the or­bitals of the up and down states by a mod­i­fied set just like for the re­stricted closed-shell case above, for each us­ing the uni­tary ma­trix of eigen­vec­tors of the $\epsilon_{mn}$ co­ef­fi­cients ap­pear­ing in the right hand side of the equa­tions for that set. It leaves the equa­tions in­tact ex­cept for changes in names, but gets rid of the equiv­a­lent of the $\epsilon_{mn}$ for $m$ $\raisebox{.2pt}{$\ne$}$ $n$, leav­ing only $\epsilon_{mm}$-equiv­a­lent val­ues. Then com­bine the spin-up and spin-down equa­tions again into a sin­gle ex­pres­sion. You get, in terms of re­vised names,

\begin{displaymath}
\fbox{$\displaystyle
h^{\rm e}\pe n////
+ \sum_{{\underli...
...////\rangle\pe {\underline n}////
=
\epsilon_n \pe n////
$}
\end{displaymath} (D.37)

In the re­stricted open-shell Hartree-Fock method, the par­ti­tion­ing also needs to in­clude the set P of or­bitals $n$ $\raisebox{-.3pt}{$\leqslant$}$ $N_p$ whose spa­tial or­bitals ap­pear both with spin-up and spin-down in the Slater de­ter­mi­nant. In that case, the pro­ce­dure above to elim­i­nate the $\epsilon_{mn}$ val­ues for $m$ $\raisebox{.2pt}{$\ne$}$ $n$ no longer works, since there are co­ef­fi­cients re­lat­ing the sets. This (even more) elab­o­rate case will be left to the ref­er­ences that you can find in [45].

Woof.