D.55 Derivation of the Hartree-Fock equations

This note derives the canonical Hartree-Fock equations. The derivation below will be performed under the normally stated rules of engagement that the orbitals are of the form $\pe{n}//u//$ or $\pe{n}//d//$, so that only the spatial orbitals $\pe{n}////$ are continuously variable. The derivations allow for the fact that some spatial spin states may be constrained to be equal.

First, you can make things a lot less messy by a priori specifying the ordering of the orbitals. The ordering makes no difference physically, and it simplifies the mathematics. In particular, in restricted Hartree-Fock some spatial orbitals appear in pairs, but you can only count each spatial orbital as one unknown function. The easiest way to handle that is to push the spin-down versions of the duplicated spatial orbitals to the end of the Slater determinant. Then the start of the determinant comprises a list of unique spatial orbitals.

So, it will be assumed that the orbitals are ordered as follows:

1. the paired spatial states in their spin-up version; assume there are $N_p\ge 0$ of them;
2. unpaired spin-up states; assume there are $N_u$ of them;
3. unpaired spin-down states; assume there are $N_d$ of them;
4. and finally, the paired spatial states in their spin-down version.

That means that the Slater-determinant wave function looks like:

$$\frac{1}{\sqrt{I!}} \Big\vert{\rm det}( \pe 1//u//,\ldot... ...,\ldots,\pe N//d//,\pe N+1//d//,\ldots,\pe I//d//)\Big\rangle$$

$$N_1=N_p+N_u \qquad N=N_p+N_u+N_d$$

$$\pe N+1////=\pe 1////,\quad \pe N+2////=\pe 2////,\quad \ldots,\quad \pe I////=\pe N_p////$$

The total number of unknown spatial orbitals is $N$, and you need a corresponding $N$ equations for them.

The variational method discussed in chapter 9.1 says that the expectation energy must be unchanged under small changes in the orbitals, provided that penalty terms are added for changes that violate the orthonormality requirements on the orbitals.

The expectation value of energy was in chapter 9.3.3 found to be:

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

(From here on, the argument of the first orbital of a pair in either side of an inner product is taken to be the first inner product integration variable ${\skew0\vec r}$ and the argument of the second orbital is the second integration variable ${\underline{\skew0\vec r}}$)

The penalty terms require penalty factors called Lagrangian variables. The penalty factor for violations of the normalization requirement

$$\langle\pe n////\vert\pe n////\rangle=1$$

will be called $\epsilon_{nn}$ for reasons evident in a second. Orthogonality between any two spatial orbitals $\pe{n}////$ and $\pe{\underline n}////$ requires

$${\textstyle\frac{1}{2}}\Big(\langle\pe n////\vert\pe{\unde... ...- \langle\pe{\underline n}////\vert\pe n////\rangle\Big)=0.$$

where the first constraint says that the real part of $\langle\pe{n}////\vert\pe{\underline n}////\rangle$ must be zero and the second that its imaginary part must be zero too. (Remember that if you switch sides in an inner product, you turn it into its complex conjugate.) To avoid including the same orthogonality condition twice, the constraint will be written only for ${\underline n}>n$. The penalty factor for the first constraint will be called $2\epsilon_{n{\underline n},r}$, and the one for the second constraint $2\epsilon_{n{\underline n},i}$.

In those terms, the penalized change in expectation energy becomes, in the restricted case that all unique spatial orbitals are mutually orthogonal,

$$\delta \langle E\rangle - \sum_{n=1}^N \sum_{{\underline... ... \delta \langle\pe n////\vert\pe{\underline n}////\rangle = 0$$

where $\epsilon_{n{\underline n}}$ is an Hermitian matrix, with $\epsilon_{n{\underline n}}=\epsilon_{n{\underline n},r}+{\rm i}\epsilon_{n{\underline n},i}$. The notations for the Lagrangian variables were chosen above to achieve this final result.

But for unrestricted Hartree-Fock, spatial orbitals are not required to be orthogonal if they have opposite spin, because the spins will take care of orthogonality. You can remove the erroneously added constraints by simply specifying that the corresponding Lagrangian variables are zero:

$$\epsilon_{n{\underline n}}=0\mbox{ if unrestricted Hartree... ...{\updownarrow}_n\vert{\updownarrow}_{\underline n}\rangle=0$}$$

or equivalently, if $n\le N_u,\; {\underline n}>N_u$ or $n>N_u,\; {\underline n}\le N_u$.

Now work out the penalized change in expectation energy due to a change in the values of a selected spatial orbital $\pe m////$ with $m\le N$. It is

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

OK, OK it is a mess. Sums like $\pe n////=\pe m////$ are for the restricted Hartree-Fock case, in which spatial orbital $\pe m////$ may appear twice in the Slater determinant. From now on, just write them as $[2]$, meaning, put in a factor 2 if orbital $\pe{m}////$ appears twice. The exception is for the exchange integrals which produce exactly one nonzero spin product; write that as $[\langle{\updownarrow}_n\vert{\updownarrow}_m\rangle^2]$, meaning take out that product if the orbital appears twice.

Next, note that the second term in each row is just the complex conjugate of the first. Considering ${\rm i}\delta\pe m////$ as a second possible change in orbital, as was done in the example in chapter 9.1, it is seen that the first terms by themselves must be zero, so you can just ignore the second term in each row. And the integrals with the factors ${\textstyle\frac{1}{2}}$ are pairwise the same; the difference is just a name swap of the first and second summation and integration variables. So all that you really have left is

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

Now note that if you write out the inner product over the first position coordinate, you will get an integral of the general form

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

If this integral is to be zero for whatever you take $\delta\pe m////$, then the terms within the parentheses must be zero. (Just take $\delta\pe m////$ proportional to the parenthetical expression; you would get the integral of an absolute square, only zero if the square is.) Unavoidably, you must have that

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

You can divide by [2]:

$$\fbox{$\displaystyle h^{\rm e}\pe m//// + \sum_{n=1}^I... ...ray} \hspace{-4pt} \right\} \epsilon_{mn}\pe n//// $}$$ (D.33)

where you use the lower of the choices between the braces in the case that spatial orbital $\pe m////$ appears twice in the Slater determinant, or equivalently, if $m\le N_p$. There you have the Hartree-Fock equations, one for each $m\le N$. Recall that they apply assuming the ordering of the Slater determinant given in the beginning, and that for unrestricted Hartree-Fock, $\epsilon_{mn}$ is zero if $\langle{\updownarrow}_m\vert{\updownarrow}_n\rangle=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 canonical Hartree-Fock equations, not just the plain vanilla version.

OK, let’s do the restricted closed-shell Hartree-Fock case first, then, since it is the easiest one. Every state is paired, so the lower choice in the curly brackets always applies, and the number of unique unknown spatial states is $N=I/2$. Also, you can reduce the summation upper limits $n=I$ to $n=N=I/2$ if you add a factor 2, since the second half of the spatial orbitals are the same as the first half. So you get

$$h^{\rm e}\pe m//// + 2 \sum_{n=1}^N \langle\pe n////\ver... ... \sum_{n=1}^N {\textstyle\frac{1}{2}}\epsilon_{mn}\pe n////$$

Now, suppose that you define a new set of orbitals, each a linear combination of the current ones:

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

where the $u_{\nu n}$ are the multiples of the original orbitals. Will the new ones still be an orthonormal set? Well, they will be if

$$\langle\overline\pe\mu////\vert\overline\pe\nu////\rangle=\delta_{\mu\nu}$$

where $\delta_{\mu\nu}$ is the Kronecker delta, one if $\mu=\nu$, zero otherwise. Substituting in the definition of the new orbitals, making sure not to use the same name for two different indices,

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

Now note that the $\pe n////$ are orthonormal, so to get a nonzero value, $m$ must be $n$, and you get

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

Consider $n$ to be the component index of a vector. Then this really says that vectors $\vec u_\mu$ and $\vec u_\nu$ must be orthonormal. So the “matrix of coefficients” must consist of orthonormal vectors. Mathematicians call such matrices “unitary,” rather than orthonormal, since it is easily confused with “unit,” and that keeps mathematicians in business explaining all the confusion.

Call the complete matrix $U$. Then, according to the rules of matrix multiplication and Hermitian adjoint, the orthonormality condition above is equivalent to $UU^{\rm {H}}=I$ where $I$ is the unit matrix. That means $U^{\rm {H}}$ is the inverse matrix to $U$, $U^{\rm {H}}=U^{-1}$ and then you also have $U^{\rm {H}}U=I$:

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

Now premultiply the definition of the new orbitals above by $U^{\rm {H}}$; you get

$$\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////$$

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

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

That gives you an expression for the original orbitals in terms of the new ones. For aesthetic reasons, you might just as well renotate $\nu$ to $\mu$, the Greek equivalent of $m$, to get

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

Now plug that into the noncanonical restricted closed-shell Hartree-Fock equations, with equivalent expressions for $\pe n////$ using whatever summation variable is still available,

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

and use the reduction formula $UU^{\rm {H}}=I$,

$$\sum_{m=1}^N u_{\nu m} u^*_{\mu m} = \delta_{\mu\nu} \qu... ...{n=1}^N u_{\kappa n} u^*_{\lambda n} = \delta_{\kappa\lambda}$$

premultiplying all by $U$, i.e. put $\sum_{m=1}^{N}u_{\nu{m}}$ before each term. You get

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

Note that the only thing that has changed more than just by symbol names is the matrix in the right hand side. Now for each value of $\mu$, take $u^*_{\mu n}$ as the $\mu$-th orthonormal eigenvector of Hermitian matrix $\epsilon_{mn}$, calling the eigenvalue $2\epsilon_\mu$. Then the right hand side becomes

$$\sum_{m=1}^{I} \sum_{\mu=1}^N u_{\nu m} \epsilon_{\mu} u... ... \overline\pe\mu//// = \epsilon_{\nu} \overline\pe\nu////$$

So, in terms of the new orbitals defined by the requirement that $u^*_{\mu n}$ gives the eigenvectors of $\epsilon_{mn}$, the right hand side simplifies to the canonical one.

Since you no longer care about the old orbitals, you can drop the overlines on the new ones, and revert to sensible roman indices $n$ and ${\underline n}$ instead of the Greek ones $\nu$ and $\lambda$. You then have the canonical restricted closed-shell Hartree-Fock equations

$$\fbox{$\displaystyle h^{\rm e}\pe n//// + 2 \sum_{{\un... .../\rangle\pe{\underline n}//// = \epsilon_n \pe n//// $}$$ (D.34)

that, as noted, assume that the Slater determinant is ordered so that the $I/2$ spin-up orbitals are at the start of it. Note that the left hand side directly provides a Hermitian Fock operator if you identify it as ${\cal F}\pe n////$; there is no need to involve spin here.

In the unrestricted case, the noncanonical equations are

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

In this case the spin-up and spin-down spatial states are not mutually orthonormal, and you want to redefine the group of spin up states and the group of spin down states separately.

The term in linear algebra is that you want to partition your $U$ matrix. What that means is simply that you separate the orbital numbers into two sets. The set of numbers $n\le N_u$ of spin-up orbitals will be indicated as U, and the set of values $n>N_u$ of spin-down ones by D. So you can partition (separate) the noncanonical equations above into equations for $m\in{\rm U}$ (meaning $m$ is one of the values in set U):

$$h^{\rm e}\pe m//// + \sum_{n\in{\rm U}} \langle\pe n////... ...ngle\pe n//// = \sum_{n\in{\rm U}} \epsilon_{mn}\pe n////$$

and equations for $m\in{\rm D}$

$$h^{\rm e}\pe m//// + \sum_{n\in{\rm U}} \langle\pe n////... ...ngle\pe n//// = \sum_{n\in{\rm D}} \epsilon_{mn}\pe n////$$

In these two equations, the fact that the up and down spin states are orthogonal was used to get rid of one pair of sums, and another pair was eliminated by the fact that there are no Lagrangian variables $\epsilon_{mn}$ linking the sets, since the spatial orbitals in the sets are allowed to be mutually nonorthogonal.

Now separately replace the orbitals of the up and down states by a modified set just like for the restricted closed-shell case above, for each using the unitary matrix of eigenvectors of the $\epsilon_{mn}$ coefficients appearing in the right hand side of the equations for that set. It leaves the equations intact except for changes in names, but gets rid of the equivalent of the $\epsilon_{mn}$ for $m\ne n$, leaving only $\epsilon_{mm}$-equivalent values. Then combine the spin-up and spin-down equations again into a single expression. You get, in terms of revised names,

$$\fbox{$\displaystyle h^{\rm e}\pe n//// + \sum_{{\unde... ...\rangle\pe {\underline n}//// = \epsilon_n \pe n//// $}$$ (D.35)

In the restricted open-shell Hartree-Fock method, the partitioning also needs to include the set P of orbitals $n\le N_p$ whose spatial orbitals appear both with spin-up and spin-down in the Slater determinant. In that case, the procedure above to eliminate the $\epsilon_{mn}$ values for $m\ne n$ no longer works, since there are coefficients relating the sets. This (even more) elaborate case will be left to the references that you can find in [44].

Woof.