A key development in quantum chemistry has been the computation of accurate self-consistent-field wave functions for many diatomic and polyatomic molecules. The principles of molecular SCF MO calculations are essentially the same as those for atomic SCF calculations (Section 11.1). We shall restrict ourselves to closed-shell configurations. For open shells, the formulas are more complicated.
The molecular Hartree-Fock wave function is written as an antisymmetrized product (Slater determinant) of spin-orbitals, each spin-orbital being a product of a spatial orbital $\phi_{i}$ and a spin function (either $\alpha$ or $\beta$ ).
The expression for the Hartree-Fock molecular electronic energy $E{\mathrm{HF}}$ is given by the variation theorem as $E{\mathrm{HF}}=\langle D| \hat{H}{\mathrm{el}}+V{N N}|D\rangle$, where $D$ is the Slater-determinant Hartree-Fock wave function, and the purely electronic Hamiltonian $\hat{H}{\mathrm{el}}$ and the internuclear repulsion $V{N N}$ are given by (13.5) and (13.6). Since $V{N N}$ doesn't involve electronic coordinates and $D$ is normalized, we have $\langle D| V{N N}|D\rangle=V{N N}\langle D \mid D\rangle=V{N N}$. The operator $\hat{H}{\text {el }}$ is the sum of one-electron operators $\hat{f}{i}$ and two electron operators $\hat{g}{i j}$; we have $\hat{H}{\mathrm{el}}=\sum{i} \hat{f}{i}+\sum{j} \sum{i>j} \hat{g}_{i j}$, where (in atomic units)
\(
\hat{f}{i}=-\frac{1}{2} \nabla{i}^{2}-\sum{\alpha} \frac{Z{\alpha}}{r{i \alpha}} \text { and } \hat{g}{i j}=\frac{1}{r_{i j}}
\)
The Hamiltonian $\hat{H}{\text {el }}$ is the same as the Hamiltonian $\hat{H}$ for an atom except that $\sum{\alpha} Z{\alpha} / r{i \alpha}$ replaces $Z / r{i}$ in $\hat{f}{i}$. Hence Eq. (11.83) can be used to give $\langle D| \hat{H}_{\text {el }}|D\rangle$. Therefore, the Hartree-Fock energy of a diatomic or polyatomic molecule with only closed shells is
\(
\begin{gather}
E{\mathrm{HF}}=2 \sum{i=1}^{n / 2} H{i i}^{\mathrm{core}}+\sum{i=1}^{n / 2} \sum{j=1}^{n / 2}\left(2 J{i j}-K{i j}\right)+V{N N} \tag{14.22}\
H{i i}^{\mathrm{core}} \equiv\left\langle\phi{i}(1)\right| \hat{H}^{\mathrm{core}}(1)\left|\phi{i}(1)\right\rangle \equiv\left\langle\phi{i}(1)\right|-\frac{1}{2} \nabla{1}^{2}-\sum{\alpha} \frac{Z{\alpha}}{r{1 \alpha}}\left|\phi{i}(1)\right\rangle \tag{14.23}\
J{i j} \equiv\left\langle\phi{i}(1) \phi{j}(2)\right| 1 / r{12}\left|\phi{i}(1) \phi{j}(2)\right\rangle, \quad K{i j} \equiv\left\langle\phi{i}(1) \phi{j}(2)\right| 1 / r{12}\left|\phi{j}(1) \phi_{i}(2)\right\rangle \tag{14.24}
\end{gather}
\)
where the one-electron-operator symbol was changed from $\hat{f}_{1}$ to $\hat{H}^{\text {core }}(1)$. The oneelectron core Hamiltonian
\(
\hat{H}^{\mathrm{core}}(1) \equiv-\frac{1}{2} \nabla{1}^{2}-\sum{\alpha} \frac{Z{\alpha}}{r{1 \alpha}}
\)
is the sum of the kinetic-energy operator for electron 1 and the potential-energy operators for the attractions between electron 1 and the nuclei. $\hat{H}^{\text {core }}(1)$ omits the interactions of electron 1 with the other electrons. The sums over $i$ and $j$ are over the $n / 2$ occupied spatial
orbitals $\phi{i}$ of the $n$-electron molecule. In the Coulomb integrals $J{i j}$ and the exchange integrals $K_{i j}$, the integration goes over the spatial coordinates of electrons 1 and 2.
The Hartree-Fock method looks for those orbitals $\phi{i}$ that minimize the variational integral $E{\mathrm{HF}}$. Each MO is taken to be normalized: $\left\langle\phi{i}(1) \mid \phi{i}(1)\right\rangle=1$. Moreover, for computational convenience one takes the MOs to be orthogonal: $\left\langle\phi{i}(1) \mid \phi{j}(1)\right\rangle=0$ for $i \neq j$. It might be thought that a lower energy could be obtained if the orthogonality restriction were omitted, but this is not so. A closed-shell antisymmetric wave function is a Slater determinant, and one can use the properties of determinants to show that a Slater determinant of nonorthogonal orbitals is equal to a Slater determinant in which the orbitals have been orthogonalized by the Schmidt or some other procedure; see Section 15.9 and F. W. Bobrowicz and W. A. Goddard, Chapter 4, Section 3.1 of Schaefer, Methods of Electronic Structure Theory. In effect, the Pauli antisymmetry requirement removes nonorthogonalities from the orbitals.
The derivation of the equation that determines the orthonormal $\phi{i}$ 's that minimize $E{\mathrm{HF}}$ is complicated and is omitted. (For the derivation, see Lowe and Peterson, Appendix 7; Szabo and Ostlund, Sections 3.1 and 3.2; Parr, pages 21-23.) One finds that the closedshell orthogonal Hartree-Fock MOs satisfy
\(
\begin{equation}
\hat{F}(1) \phi{i}(1)=\varepsilon{i} \phi_{i}(1) \tag{14.25}
\end{equation}
\)
where $\varepsilon_{i}$ is the orbital energy and where the (Hartree-) Fock operator $\hat{F}$ is (in atomic units)
\(
\begin{gather}
\hat{F}(1)=\hat{H}^{\text {core }}(1)+\sum{j=1}^{n / 2}\left[2 \hat{J}{j}(1)-\hat{K}{j}(1)\right] \tag{14.26}\
\hat{H}^{\text {core }}(1) \equiv-\frac{1}{2} \nabla{1}^{2}-\sum{\alpha} \frac{Z{\alpha}}{r_{1 \alpha}} \tag{14.27}
\end{gather}
\)
where the Coulomb operator $\hat{J}{j}$ and the exchange operator $\hat{K}{j}$ are defined by
\(
\begin{align}
\hat{J}{j}(1) f(1) & =f(1) \int\left|\phi{j}(2)\right|^{2} \frac{1}{r{12}} d v{2} \tag{14.28}\
\hat{K}{j}(1) f(1) & =\phi{j}(1) \int \frac{\phi_{j}^{}(2) f(2)}{r{12}} d v{2} \tag{14.29}
\end{align*}
\)
where $f$ is an arbitrary function and the integrals are definite integrals over all space.
The first term on the right of (14.27) is the operator for the kinetic energy of one electron. The second term is the sum of the potential-energy operators for the attractions between one electron and the nuclei. The Coulomb operator $\widehat{J}{j}(1)$ is the potential energy of interaction between electron 1 and a smeared-out electron with electronic density $-\left|\phi{j}(2)\right|^{2}$. The factor 2 in (14.26) occurs because there are two electrons in each spatial orbital. The exchange operator has no simple physical interpretation but arises from the requirement that the wave function be antisymmetric with respect to electron exchange. The exchange operators are absent from the Hartree equations (11.9). The Hartree-Fock MOs $\phi{i}$ in (14.25) are eigenfunctions of the same operator $\hat{F}$, the eigenvalues being the orbital energies $\varepsilon{i}$.
The orthogonality of the MOs greatly simplifies MO calculations, causing many integrals to vanish. In contrast, the VB method uses atomic orbitals, and AOs centered on different atoms are not orthogonal. MO calculations are much simpler than VB calculations, and the MO method is used far more often than the VB method.
The true Hamiltonian operator and wave function involve the coordinates of all $n$ electrons. The Hartree-Fock Hamiltonian operator $\hat{F}$ is a one-electron operator (that is, it involves the coordinates of only one electron), and (14.25) is a one-electron differential equation. This has been indicated in (14.25) by writing $\hat{F}$ and $\phi_{i}$ as functions of the coordinates of electron 1 . Of course, the coordinates of any electron could have been used. The operator $\hat{F}$ is peculiar in that it depends on its own eigenfunctions [see Eqs. (14.26) to (14.29)], which are not known initially. Hence the HartreeFock equations must be solved by an iterative process, starting with an initial guess for the MOs.
To obtain the expression for the orbital energies $\varepsilon{i}$, we multiply (14.25) by $\phi{i}^{}(1)$ and integrate over all space. Using the fact that $\phi{i}$ is normalized and using the result of Prob. 14.8, we get $\varepsilon{i}=\int \phi_{i}^{}(1) \hat{F}(1) \phi{i}(1) d v{1}$ and
\(
\begin{gather}
\varepsilon{i}=\left\langle\phi{i}(1)\right| \hat{H}^{\mathrm{core}}(1)\left|\phi{i}(1)\right\rangle+\sum{j}\left[2\left\langle\phi{i}(1)\right| \hat{J}{j}(1)\left|\phi{i}(1)\right\rangle-\left\langle\phi{i}(1)\right| \hat{K}{j}(1)\left|\phi{i}(1)\right\rangle\right] \
\varepsilon{i}=H{i i}^{\text {core }}+\sum{j=1}^{n / 2}\left(2 J{i j}-K_{i j}\right) \tag{14.30}
\end{gather}
\)
where $H{i i}^{\text {core }}, J{i j}$, and $K_{i j}$ are defined by (14.23) and (14.24).
Summation of (14.30) over the $n / 2$ occupied orbitals gives
\(
\begin{equation}
\sum{i=1}^{n / 2} \varepsilon{i}=\sum{i=1}^{n / 2} H{i i}^{\text {core }}+\sum{i=1}^{n / 2} \sum{j=1}^{n / 2}\left(2 J{i j}-K{i j}\right) \tag{14.31}
\end{equation}
\)
Solving this equation for $\sum{i} H{i i}^{\text {core }}$ and substituting the result into (14.22), we obtain the Hartree-Fock energy as
\(
\begin{equation}
E{\mathrm{HF}}=2 \sum{i=1}^{n / 2} \varepsilon{i}-\sum{i=1}^{n / 2} \sum{j=1}^{n / 2}\left(2 J{i j}-K{i j}\right)+V{N N} \tag{14.32}
\end{equation}
\)
Since there are two electrons per MO, the quantity $2 \sum{i} \varepsilon{i}$ is the sum of the orbital energies. Subtraction of the double sum in (14.32) avoids counting each interelectronic repulsion twice, as discussed in Section 11.1.
A key development that made feasible the calculation of accurate molecular SCF wave functions was Roothaan's 1951 proposal to expand the spatial orbitals $\phi{i}$ as linear combinations of a set of one-electron basis functions $\chi{s}$ :
\(
\begin{equation}
\phi{i}=\sum{s=1}^{b} c{s i} \chi{s} \tag{14.33}
\end{equation}
\)
To exactly represent the MOs $\phi{i}$, the basis functions $\chi{s}$ should form a complete set. This requires an infinite number of basis functions. In practice, one must use a finite number $b$ of basis functions. If $b$ is large enough and the functions $\chi_{s}$ are well chosen, one can represent the MOs with negligible error.
To avoid confusion, we shall use the letters $r, s, t, u$ to label the basis functions $\chi$, and the letters $i, j, k, l$ to label the MOs $\phi$. (Often the Greek letters $\mu, \nu, \lambda, \sigma$ are used to label the basis functions.)
Substitution of the expansion (14.33) into the Hartree-Fock equations (14.25) gives
\(
\sum{s} c{s i} \hat{F} \chi{s}=\varepsilon{i} \sum{s} c{s i} \chi_{s}
\)
Multiplication by $\chi_{r}^{*}$ and integration gives
\(
\begin{gather}
\sum{s=1}^{b} c{s i}\left(F{r s}-\varepsilon{i} S{r s}\right)=0, \quad r=1,2, \ldots, b \tag{14.34}\
F{r s} \equiv\left\langle\chi{r}\right| \hat{F}\left|\chi{s}\right\rangle, \quad S{r s} \equiv\left\langle\chi{r} \mid \chi_{s}\right\rangle \tag{14.35}
\end{gather}
\)
The equations (14.34) form a set of $b$ simultaneous linear homogeneous equations in the $b$ unknowns $c{s i}, s=1,2, \ldots, b$, that describe the $\mathrm{MO} \phi{i}$ in (14.33). For a nontrivial solution, we must have
\(
\begin{equation}
\operatorname{det}\left(F{r s}-\varepsilon{i} S_{r s}\right)=0 \tag{14.36}
\end{equation}
\)
This is a secular equation whose roots give the orbital energies $\varepsilon{i}$. The (Hartree-Fock-) Roothaan equations (14.34) must be solved by an iterative process, since the $F{r s}$ integrals depend on the orbitals $\phi{i}$ (through the dependence of $\hat{F}$ on the $\phi{i}$ 's), which in turn depend on the unknown coefficients $c_{s i}$.
One starts with guesses for the occupied-MO expressions as linear combinations of the basis functions, as in (14.33). This initial set of MOs is used to compute the Fock operator $\hat{F}$ from (14.26) to (14.29). The matrix elements (14.35) are computed, and the secular equation (14.36) is solved to give an initial set of $\varepsilon{i}$ 's. These $\varepsilon{i}$ 's are used to solve (14.34) for an improved set of coefficients, giving an improved set of MOs, which are then used to compute an improved $\hat{F}$, and so on. One continues until no further improvement in MO coefficients and energies occurs from one cycle to the next. The calculations are done using a computer. (The most efficient way to solve the Roothaan equations is to use matrix-algebra methods; see the last part of this section.)
We have used the terms SCF wave function and Hartree-Fock wave function interchangeably. In practice, the term SCF wave function is applied to any wave function obtained by iterative solution of the Roothaan equations, whether or not the basis set is large enough to give a really accurate approximation to the Hartree-Fock SCF wave function. There is only one true Hartree-Fock wave function, which is the best possible wave function that can be written as a Slater determinant of spin-orbitals. With current computer power, one can use very large basis sets for small molecules and obtain wave functions that differ negligibly from the true Hartree-Fock wave functions. Because of deficiencies in properties calculated from Hartree-Fock wave functions, several methods that go beyond the Hartree-Fock method are widely used (see Chapter 16).
The Fock Matrix Elements
To solve the Roothaan equations (14.34), we first must express the Fock matrix elements (integrals) $F_{r s}$ in terms of the basis functions $\chi$. The Fock operator $\hat{F}$ is given by (14.26), so
\(
\begin{align}
& F{r s}=\left\langle\chi{r}(1)\right| \hat{F}(1)\left|\chi{s}(1)\right\rangle \
& F{r s}=\left\langle\chi{r}(1)\right| \hat{H}^{\mathrm{core}}(1)\left|\chi{s}(1)\right\rangle+\sum{j=1}^{n / 2}\left[2\left\langle\chi{r}(1) \mid \hat{J}{j}(1) \chi{s}(1)\right\rangle-\left\langle\chi{r}(1) \mid \hat{K}{j}(1) \chi_{s}(1)\right\rangle\right] \tag{14.37}
\end{align}
\)
Replacement of $f$ by $\chi_{s}$ in (14.28), followed by use of the expansion (14.33), gives
\(
\hat{J}{j}(1) \chi{s}(1)=\chi{s}(1) \int \frac{\phi{j}^{}(2) \phi{j}(2)}{r{12}} d v{2}=\chi{s}(1) \sum{t} \sum{u} c_{t j}^{} c{u j} \int \frac{\chi{t}^{*}(2) \chi{u}(2)}{r{12}} d v_{2}
\)
Multiplication by $\chi_{r}^{*}(1)$ and integration over the coordinates of electron 1 gives
\(
\left\langle\chi{r}(1) \mid \hat{J}{j}(1) \chi{s}(1)\right\rangle=\sum{t} \sum{u} c{t j}^{} c{u j} \iint \frac{\chi{r}^{}(1) \chi{s}(1) \chi{t}^{*}(2) \chi{u}(2)}{r{12}} d v{1} d v{2}
\)
\(
\begin{equation}
\left\langle\chi{r}(1) \mid \hat{J}{j}(1) \chi{s}(1)\right\rangle=\sum{t=1}^{b} \sum{u=1}^{b} c{t j}^{} c_{u j}(r s \mid t u) \tag{14.38}
\end{}
\)
where the two-electron repulsion integral $(r s \mid t u)$ is defined as
\(
\begin{equation}
(r s \mid t u) \equiv \iint \frac{\chi_{r}^{}(1) \chi{s}(1) \chi{t}^{}(2) \chi{u}(2)}{r{12}} d v{1} d v{2} \tag{14.39}
\end{equation}
\)
The widely used notation of (14.39) should not be misinterpreted as an overlap integral. Other notations, some of which are mutually contradictory, are used for electron repulsion integrals, so it is always wise to check an author's definition.
Similarly, replacement of $f$ by $\chi_{s}$ in (14.29) leads to (Prob. 14.9)
\(
\begin{equation}
\left\langle\chi{r}(1) \mid \hat{K}{j}(1) \chi{s}(1)\right\rangle=\sum{t=1}^{b} \sum{u=1}^{b} c{t j}^{} c_{u j}(r u \mid t s) \tag{14.40}
\end{}
\)
Substituting (14.40) and (14.38) into (14.37) and changing the order of summation, we get the desired expression for $F_{r s}$ in terms of integrals over the basis functions $\chi$ :
\(
\begin{gather}
F{r s}=H{r s}^{\mathrm{core}}+\sum{t=1}^{b} \sum{u=1}^{b} \sum{j=1}^{n / 2} c{t j}^{} c{u j}[2(r s \mid t u)-(r u \mid t s)] \
F{r s}=H{r s}^{\mathrm{core}}+\sum{t=1}^{b} \sum{u=1}^{b} P{t u}\left[(r s \mid t u)-\frac{1}{2}(r u \mid t s)\right] \tag{14.41}\
P{t u} \equiv 2 \sum{j=1}^{n / 2} c{t j}^{*} c{u j}, \quad t=1,2, \ldots, b, \quad u=1,2, \ldots, b \tag{14.42}\
H{r s}^{\text {core }} \equiv\left\langle\chi{r}(1)\right| \hat{H}^{\text {core }}(1)\left|\chi_{s}(1)\right\rangle
\end{gather*}
\)
The quantities $P{t u}$ are called density matrix elements or charge, bond-order matrix elements. [Some workers use the definition $P{u t} \equiv 2 \sum{j} c{t j}^{*} c_{u j}$.] Substitution of the expansion (14.33) into (14.7) for the electron probability density $\rho$ gives for a closed-shell molecule:
\(
\begin{equation}
\rho=2 \sum{j=1}^{n / 2} \phi{j}^{} \phi{j}=2 \sum{r=1}^{b} \sum{s=1}^{b} \sum{j=1}^{n / 2} c{r j}^{*} c{s j} \chi{r}^{*} \chi{s}=\sum{r=1}^{b} \sum{s=1}^{b} P{r s} \chi{r}^{} \chi_{s} \tag{14.43}
\end{equation}
\)
To express the Hartree-Fock energy in terms of integrals over the basis functions $\chi$, we first solve (14.31) for $\sum{i} \sum{j}\left(2 J{i j}-K{i j}\right)$ and substitute the result into (14.32) to get
\(
E{\mathrm{HF}}=\sum{i=1}^{n / 2} \varepsilon{i}+\sum{i=1}^{n / 2} H{i i}^{\mathrm{core}}+V{N N}
\)
We have, using the expansion (14.33),
\(
\begin{gather}
H{i i}^{\text {core }}=\left\langle\phi{i}\right| \hat{H}^{\mathrm{core}}\left|\phi{i}\right\rangle=\sum{r} \sum{s} c{r i}^{} c{s i}\left\langle\chi{r}\right| \hat{H}^{\mathrm{core}}\left|\chi{s}\right\rangle=\sum{r} \sum{s} c{r i}^{} c{s i} H{r s}^{\text {core }} \
E{\mathrm{HF}}=\sum{i=1}^{n / 2} \varepsilon{i}+\sum{r} \sum{s} \sum{i=1}^{n / 2} c_{r i}^{} c{s i} H{r s}^{\mathrm{core}}+V{N N} \
E{\mathrm{HF}}=\sum{i=1}^{n / 2} \varepsilon{i}+\frac{1}{2} \sum{r=1}^{b} \sum{s=1}^{b} P{r s} H{r s}^{\mathrm{core}}+V_{N N} \tag{14.44}
\end{gather*}
\)
An alternative expression for $E{\mathrm{HF}}$ is useful. Multiplication of $\hat{F} \phi{i}=\varepsilon{i} \phi{i}$ [Eq. (14.25)] by $\phi{i}^{*}$ and integration gives $\varepsilon{i}=\left\langle\phi{i}\right| \hat{F}\left|\phi{i}\right\rangle$. Substitution of $\phi{i}=\sum{s=1}^{b} c{s i} \chi{s}$ [Eq. (14.33)] gives $\varepsilon{i}=\sum{r} \sum{s} c{r i}^{} c{s i}\left\langle\chi{r}\right| \hat{F}\left|\chi{s}\right\rangle=\sum{r} \sum{s} c{r i}^{} c{s i} F{r s}$. The first sum in (14.44) becomes $\sum{i} \varepsilon{i}=\sum{r} \sum{s} \sum{i} c{r i}^{*} c{s i} F{r s}=\frac{1}{2} \sum{r} \sum{s} P{r s} F{r s}$, where the definition (14.42) of $P_{r s}$ was used. Equation (14.44) becomes
\(
\begin{equation}
E{\mathrm{HF}}=\frac{1}{2} \sum{r=1}^{b} \sum{s=1}^{b} P{r s}\left(F{r s}+H{r s}^{\mathrm{core}}\right)+V_{N N} \tag{14.45}
\end{equation}
\)
which expresses $E{\mathrm{HF}}$ of a closed-shell molecule in terms of the density, Fock, and coreHamiltonian matrix elements calculated with the basis functions $\chi{r}$.
EXAMPLE
Do an SCF calculation for the helium-atom ground state using a basis set of two $1 s$ STOs with orbital exponents $\zeta{1}=1.45$ and $\zeta{2}=2.91$. [By trial and error, these have been found to be the optimum $\zeta$ 's to use for this basis set; see C. Roetti and E. Clementi, J. Chem. Phys., 60, 4725 (1974).]
From (11.14), the normalized basis functions are (in atomic units)
\(
\begin{equation}
\chi{1}=2 \zeta{1}^{3 / 2} e^{-\zeta{1} r} Y{0}^{0}, \quad \chi{2}=2 \zeta{2}^{3 / 2} e^{-\zeta{2} r} Y{0}^{0}, \quad \zeta{1}=1.45, \quad \zeta{2}=2.91 \tag{14.46}
\end{equation}
\)
To solve the Roothaan equations (14.34), we need the integrals $F{r s}$ and $S{r s}$. The overlap integrals $S_{r s}$ are
\(
\begin{gathered}
S{11}=\left\langle\chi{1} \mid \chi{1}\right\rangle=1, \quad S{22}=\left\langle\chi{2} \mid \chi{2}\right\rangle=1 \
S{12}=S{21}=\left\langle\chi{1} \mid \chi{2}\right\rangle=4 \zeta{1}^{3 / 2} \zeta{2}^{3 / 2} \int{0}^{\infty} e^{-\left(\zeta{1}+\zeta{2}\right) r^{2}} d r=\frac{8 \zeta{1}^{3 / 2} \zeta{2}^{3 / 2}}{\left(\zeta{1}+\zeta_{2}\right)^{3}}=0.8366
\end{gathered}
\)
where the Appendix integral (A.8) was used.
The integrals $F{r s}$ are given by (14.41) and depend on $H{r s}^{\text {core }}, P{t u}$, and $(r s \mid t u)$. From (14.27), $\hat{H}^{\text {core }}=-\frac{1}{2} \nabla^{2}-2 / r=-\frac{1}{2} \nabla^{2}-\zeta / r+(\zeta-2) / r$. The integrals $H{r s}^{\text {core }}$ are evaluated the same way that similar integrals were evaluated in the variation treatment of He in Section 9.4. We find (Prob. 14.12)
\(
\begin{gathered}
H{11}^{\text {core }}=\left\langle\chi{1}\right| \hat{H}^{\text {core }}\left|\chi{1}\right\rangle=-\frac{1}{2} \zeta{1}^{2}+\left(\zeta{1}-2\right) \zeta{1}=\frac{1}{2} \zeta{1}^{2}-2 \zeta{1}=-1.8488 \
H{22}^{\text {core }}=\frac{1}{2} \zeta{2}^{2}-2 \zeta{2}=-1.5860 \
H{12}^{\text {core }}=H{21}^{\text {core }}=\left\langle\chi{1}\right| \hat{H}^{\text {core }}\left|\chi{2}\right\rangle=-\frac{1}{2} \zeta{2}^{2} S{12}+\frac{4\left(\zeta{2}-2\right) \zeta{1}^{3 / 2} \zeta{2}^{3 / 2}}{\left(\zeta{1}+\zeta{2}\right)^{2}} \
H{12}^{\text {core }}=H{21}^{\text {core }}=\frac{\zeta{1}^{3 / 2} \zeta{2}^{3 / 2}\left(4 \zeta{1} \zeta{2}-8 \zeta{1}-8 \zeta{2}\right)}{\left(\zeta{1}+\zeta{2}\right)^{3}}=-1.8826
\end{gathered}
\)
Many of the electron-repulsion integrals ( $r s \mid t u$ ) are equal to one another. For real basis functions, one can show that (Prob. 14.13)
\(
\begin{equation}
(r s \mid t u)=(s r \mid t u)=(r s \mid u t)=(s r \mid u t)=(t u \mid r s)=(u t \mid r s)=(t u \mid s r)=(u t \mid s r) \tag{14.47}
\end{equation}
\)
The electron-repulsion integrals are evaluated using the $1 / r_{12}$ expansion (9.124) in Prob. 9.14. One finds [see Eq. (9.53) and Prob. 14.14]
\(
\begin{aligned}
&(11 \mid 11)= \frac{5}{8} \zeta{1}=0.9062, \quad(22 \mid 22)=\frac{5}{8} \zeta{2}=1.8188 \
&(11 \mid 22)=(22 \mid 11)=\left(\zeta{1}^{4} \zeta{2}+4 \zeta{1}^{3} \zeta{2}^{2}+\zeta{1} \zeta{2}^{4}+4 \zeta{1}^{2} \zeta{2}^{3}\right) /\left(\zeta{1}+\zeta{2}\right)^{4}=1.1826 \
&(12 \mid 12)=(21 \mid 12)=(12 \mid 21)=(21 \mid 21)=20 \zeta{1}^{3} \zeta{2}^{3} /\left(\zeta{1}+\zeta{2}\right)^{5}=0.9536 \
&(11 \mid 12)=(11 \mid 21)=(12 \mid 11)=(21 \mid 11)=\frac{16 \zeta{1}^{9 / 2} \zeta{2}^{3 / 2}}{\left(3 \zeta{1}+\zeta{2}\right)^{4}}\left[\frac{12 \zeta{1}+8 \zeta{2}}{\left(\zeta{1}+\zeta{2}\right)^{2}}+\frac{9 \zeta{1}+\zeta{2}}{2 \zeta_{1}^{2}}\right]=0.9033 \
&(12 \mid 22)=(22 \mid 12)=(21 \mid 22)=(22 \mid 21) \
& \quad=\text { the }(11 \mid 12) \text { expression with } 1 \text { and } 2 \text { interchanged }=1.2980
\end{aligned}
\)
To start the calculation, we need an initial guess for the ground-state AO expansion coefficients $c{s i}$ in (14.33) so that we can get an initial estimate of the density matrix elements $P{t u}$ in (14.41). We saw in Section 9.4 that the optimum orbital exponent for a helium AO that consists of one $1 s$ STO is $\frac{27}{16}=1.6875$. Since the orbital exponent $\zeta{1}$ is much closer to 1.6875 than is $\zeta{2}$, we expect that the coefficient of $\chi{1}$ in $\phi{1}=c{11} \chi{1}+c{21} \chi{2}$ will be substantially larger than the coefficient of $\chi{2}$. Let us take as an initial guess $c{11} / c{21} \approx 2$. [A more general method to get an initial guess for the $c{s i}$ coefficients is to neglect the electron-repulsion integrals in (14.41) and approximate $F{r s}$ in the secular equation (14.36) by $F{r s} \approx H{r s}^{\text {core }}$; we then solve (14.36) and (14.34). This would give $c{11} / c{21} \approx 1.5$ (Prob 14.15).] The normalization condition $\int\left|\phi{1}\right|^{2} d \tau=1$ gives for real coefficients (Prob. 14.17)
\(
\begin{equation}
c{21}=\left(1+k^{2}+2 k S{12}\right)^{-1 / 2}, \quad \text { where } k \equiv c{11} / c{21} \tag{14.48}
\end{equation}
\)
Substitution of $k=2$ and $S{12}=0.8366$ gives $c{21} \approx 0.3461$ and $c{11} \approx 2 c{21}=0.6922$.
With $n=2$ and $b=2$, Eq. (14.42) gives
\(
\begin{equation}
P{11}=2 c{11}^{} c{11}, \quad P{12}=2 c{11}^{*} c{21}, \quad P{21}=P{12}^{}, \quad P{22}=2 c{21}^{} c_{21} \tag{14.49}
\end{}
\)
The initial guess $c{11} \approx 0.6922, c{21} \approx 0.3461$ gives as the initial density matrix elements
\(
P{11} \approx 0.9583, \quad P{12}=P{21} \approx 0.4791, \quad P{22} \approx 0.2396
\)
The Fock matrix elements are found from (14.41) with $b=2$. Using (14.47) and $P{12}=P{21}$ for real functions, we get (Prob. 14.16a)
\(
\begin{gathered}
F{11}=H{11}^{\text {core }}+\frac{1}{2} P{11}(11 \mid 11)+P{12}(11 \mid 12)+P{22}\left[(11 \mid 22)-\frac{1}{2}(12 \mid 21)\right] \
F{12}=F{21}=H{12}^{\text {core }}+\frac{1}{2} P{11}(12 \mid 11)+P{12}\left[\frac{3}{2}(12 \mid 12)-\frac{1}{2}(11 \mid 22)\right]+\frac{1}{2} P{22}(12 \mid 22) \
F{22}=H{22}^{\text {core }}+P{11}\left[(22 \mid 11)-\frac{1}{2}(21 \mid 12)\right]+P{12}(22 \mid 12)+\frac{1}{2} P{22}(22 \mid 22)
\end{gathered}
\)
Substitution of the values of the $H_{r s}^{\text {core }}$ and $(r s \mid t u)$ integrals listed previously gives (Prob. 14.16b)
\(
\begin{gather}
F{11}=-1.8488+0.4531 P{11}+0.9033 P{12}+0.7058 P{22} \tag{14.50}\
F{12}=F{21}=-1.8826+0.4516{5} P{11}+0.8391 P{12}+0.6490 P{22} \tag{14.51}\
F{22}=-1.5860+0.7058 P{11}+1.2980 P{12}+0.9094 P{22} \tag{14.52}
\end{gather}
\)
Substitution of the initial guess for the $P{t u}$ 's into (14.50) to (14.52) gives as the initial estimates of the $F{r s}$ matrix elements:
\(
F{11} \approx-0.813, \quad F{12}=F{21} \approx-0.892, \quad F{22} \approx-0.070
\)
The initial estimate of the secular equation $\operatorname{det}\left(F{r s}-S{r s} \varepsilon_{i}\right)=0$ is
\(
\begin{gathered}
\left|\begin{array}{cc}
-0.813-\varepsilon{i} & -0.892-0.8366 \varepsilon{i} \
-0.892-0.8366 \varepsilon{i} & -0.070-\varepsilon{i}
\end{array}\right| \approx 0 \
0.3001 \varepsilon{i}^{2}-0.609{5} \varepsilon{i}-0.739 \approx 0 \
\varepsilon{1} \approx-0.854, \quad \varepsilon_{2} \approx 2.885
\end{gathered}
\)
Substitution of the lower root $\varepsilon_{1}$ into the Roothaan equation (14.34) with $r=2$ gives
\(
\begin{gathered}
c{11}\left(F{21}-\varepsilon{1} S{21}\right)+c{21}\left(F{22}-\varepsilon{1} S{22}\right) \approx 0 \
-0.177{5} c{11}+0.784 c{21} \approx 0 \
c{11} / c_{21} \approx 4.42
\end{gathered}
\)
Substitution of $k=4.42$ and $S_{12}=0.8366$ in the normalization condition (14.48) gives
\(
c{21} \approx 0.189, \quad c{11}=k c_{21} \approx 0.836
\)
Substitution of these improved coefficients into (14.49) gives as the improved density matrix elements
\(
P{11} \approx 1.398, \quad P{12}=P{21} \approx 0.316, \quad P{22} \approx 0.071
\)
Substitution of these improved $P{t u}$ 's into (14.50) to (14.52) gives as the improved $F{r s}$ values
\(
F{11} \approx-0.880, \quad F{12}=F{21} \approx-0.940, \quad F{22} \approx-0.124_{6}
\)
The improved secular equation is
\(
\begin{gathered}
\left|\begin{array}{cc}
-0.880-\varepsilon{i} & -0.940-0.8366 \varepsilon{i} \
-0.940-0.8366 \varepsilon{i} & -0.124{6}-\varepsilon{i}
\end{array}\right| \approx 0 \
\varepsilon{1} \approx-0.918, \quad \varepsilon_{2} \approx 2.810
\end{gathered}
\)
The improved $\varepsilon{1}$ value gives $c{11} / c_{21} \approx 4.61$ and
\(
c{11} \approx 0.842, \quad c{21} \approx 0.183
\)
Another cycle of calculation yields (Prob. 14.18)
\(
\begin{gather}
P{11}=1.418, \quad P{12}=P{21}=0.308, \quad P{22}=0.067 \
F{11}=-0.881, \quad F{12}=F{21}=-0.940, \quad F{22}=-0.124{5} \tag{14.53}\
\varepsilon{1}=-0.918, \quad \varepsilon{2}=2.809 \tag{14.54}\
c{11}=0.842, \quad c_{21}=0.183
\end{gather}
\)
These last $c$ 's are the same as those for the previous cycle, so the calculation has converged and we are finished. The He ground-state SCF AO for this basis set is
\(
\phi{1}=0.842 \chi{1}+0.183 \chi_{2}
\)
The SCF energy is found from (14.44) with $n=2$ and $b=2$ as
\(
\begin{aligned}
E_{\mathrm{HF}} & =-0.918+\frac{1}{2}[1.418(-1.8488)+2(0.308)(-1.8826)+0.067(-1.5860)]+0 \
& =-2.862 \text { hartrees }=-77.9 \mathrm{eV}
\end{aligned}
\)
A more precise calculation with $\zeta{1}=1.45363$ and $\zeta{2}=2.91093$ gives an SCF energy of -2.8616726 hartrees, as compared with the limiting Hartree-Fock energy -2.8616799 hartrees found with five basis functions [C. Roetti and E. Clementi, J. Chem. Phys., 60, 4725 (1974)].
Matrix Form of the Roothaan Equations
The Roothaan equations are most efficiently solved using matrix methods. The Roothaan equations (14.34) read
\(
\sum{s=1}^{b} F{r s} c{s i}=\sum{s=1}^{b} S{r s} c{s i} \varepsilon_{i}, \quad r=1,2, \ldots, b
\)
The coefficients $c{s i}$ relate the MOs $\phi{i}$ to the basis functions $\chi{s}$ according to $\phi{i}=\sum{s} c{s i} \chi{s}$. Let $\mathbf{C}$ be the square matrix of order $b$ whose elements are the coefficients $c{s i}$. Let $\mathbf{F}$ be the square matrix of order $b$ whose elements are $F{r s}=\left\langle\chi{r}\right| \hat{F}\left|\chi{s}\right\rangle$. Let $\mathbf{S}$ be the square matrix whose elements are $S{r s}=\left\langle\chi{r} \mid \chi{s}\right\rangle$. Let $\boldsymbol{\varepsilon}$ be the diagonal square matrix whose diagonal elements are the orbital energies $\varepsilon{1}, \varepsilon{2}, \ldots, \varepsilon{b}$ so that the elements of $\boldsymbol{\varepsilon}$ are $\varepsilon{m i}=\delta{m i} \varepsilon{i}$, where $\delta_{m i}$ is the Kronecker delta.
Use of the matrix multiplication rule (7.107) gives the $(s, i)$ th element of the matrix product $\mathbf{C} \boldsymbol{\varepsilon}$ as $(\mathbf{C} \boldsymbol{\varepsilon}){s i}=\sum{m} c{s m} \varepsilon{m i}=\sum{m} c{s m} \delta{m i} \varepsilon{i}=c{s i} \varepsilon{i}$. Hence the Roothaan equations read
\(
\begin{equation}
\sum{s=1}^{b} F{r s} c{s i}=\sum{s=1}^{b} S{r s}(\mathbf{C} \boldsymbol{\varepsilon}){s i} \tag{14.55}
\end{equation}
\)
From the matrix multiplication rule, the left side of (14.55) is the $(r, i)$ th element of $\mathbf{F C}$, and the right side is the $(r, i)$ th element of $\mathbf{S}(\mathbf{C} \boldsymbol{\varepsilon})$. Since the general element of $\mathbf{F C}$ equals the general element of $\mathbf{S C \varepsilon}$, these matrices are equal:
\(
\begin{equation}
\mathrm{FC}=\mathbf{S C} \boldsymbol{\varepsilon} \tag{14.56}
\end{equation}
\)
This is the matrix form of the Roothaan equations.
The set of basis functions $\chi{s}$ used to expand the MOs is not an orthogonal set. However, one can use the Schmidt or some other procedure to form orthogonal linear combinations of the basis functions to give a new set of basis functions $\chi{s}^{\prime}$ that is an orthonormal set: $\chi{s}^{\prime}=\sum{t} a{t s} \chi{t}$ and $S{r s}^{\prime}=\left\langle\chi{r}^{\prime} \mid \chi{s}^{\prime}\right\rangle=\delta{r s}$. (See Probs. 8.57 and 8.58 and Szabo and Ostlund, Section 3.4.5, for details of the orthogonalization procedure.) With this orthonormal basis set, the overlap matrix is a unit matrix, and the Roothaan equations (14.56) have the simpler form
\(
\begin{equation}
\mathbf{F}^{\prime} \mathbf{C}^{\prime}=\mathbf{C}^{\prime} \boldsymbol{\varepsilon} \tag{14.57}
\end{equation}
\)
where $F{r s}^{\prime}=\left\langle\chi{r}^{\prime}\right| \hat{F}\left|\chi{s}^{\prime}\right\rangle$ and $\mathbf{C}^{\prime}$ is the matrix of the coefficients that relate the MOs $\phi{i}$ to the orthonormal basis functions: $\phi{i}=\sum{s} c{s i}^{\prime} \chi{s}^{\prime}$. It was shown in Prob. 8.57c that the $\mathbf{F}$ and $\mathbf{F}^{\prime}$ matrices and the $\mathbf{C}$ and $\mathbf{C}^{\prime}$ matrices are related by
\(
\mathbf{F}^{\prime}=\mathbf{A}^{\dagger} \mathbf{F A} \text { and } \mathbf{C}=\mathbf{A C}^{\prime}
\)
where $\mathbf{A}$ is the matrix of coefficients $a{t s}$ in $\chi{s}^{\prime}=\sum{s} a{t s} \chi_{t}$, so we can readily calculate $\mathbf{F}^{\prime}$ from $\mathbf{F}$ and $\mathbf{C}$ from $\mathbf{C}^{\prime}$. [ $\mathbf{H}$ in Prob. 8.57 corresponds to $\mathbf{F}$ in (14.56).]
The matrix equation (14.57) has the same form as Eq. (8.87), which is $\mathbf{H C}=\mathbf{C W}$, where $\mathbf{C}$ and $\mathbf{W}$ [defined by (8.86)] are the eigenvector matrix and eigenvalue matrix, respectively, of $\mathbf{H}$. Thus, the orbital energies $\varepsilon{i}$ are the eigenvalues of the Fock matrix $\mathbf{F}^{\prime}$ and each column of $\mathbf{C}^{\prime}$ is an eigenvector of $\mathbf{F}^{\prime}$. Because the Fock operator $\hat{F}$ is Hermitian, the Fock matrix $\mathbf{F}^{\prime}$ is a Hermitian matrix. As noted in the paragraph preceding Eq. (8.94), the eigenvector matrix $\mathbf{C}^{\prime}$ of the Hermitian matrix $\mathbf{F}^{\prime}$ can be chosen to be unitary, meaning that its inverse equals its conjugate transpose [Eq. (8.92)] $\mathbf{C}^{\prime-1}=\mathbf{C}^{\prime \dagger}$. (With a unitary coefficient matrix $\mathbf{C}^{\prime}$, the MOs $\phi{i}$ are orthonormal; see Prob. 14.22.) Multiplication of (14.57) on the left by $\mathbf{C}^{\prime-1}=\mathbf{C}^{\prime \dagger}$ gives [see Eqs. (8.88) and (8.94)]
\(
\begin{equation}
\mathbf{C}^{\prime \dagger} \mathbf{F}^{\prime} \mathbf{C}^{\prime}=\boldsymbol{\varepsilon} \tag{14.58}
\end{equation}
\)
which has the same form as Eq. (8.94).
The following procedure is commonly used to do an SCF MO calculation at a specified molecular geometry.
- Choose a basis set $\chi_{s}$.
- Evaluate the $H{r s}^{\text {core }}, S{r s}$, and ( $r s \mid t u$ ) integrals.
- Use the overlap integrals $S{r s}$ and an orthogonalization procedure to calculate the $\mathbf{A}$ matrix of coefficients $a{t s}$ that will produce orthonormal basis functions $\chi{s}^{\prime}=\sum{t} a{t s} \chi{t}$.
- Make an initial guess for the coefficients $c{s i}$ in the MOs $\phi{i}=\sum{s} c{s i} \chi_{s}$. From the initial guess of coefficients, calculate the density matrix $\mathbf{P}$ in (14.42).
- Use (14.41) to calculate an estimate of the Fock matrix elements $F{r s}$ from $\mathbf{P}$ and the $(r s \mid t u)$ and $H{r s}^{\text {core }}$ integrals.
- Calculate the matrix $\mathbf{F}^{\prime}$ using $\mathbf{F}^{\prime}=\mathbf{A}^{\dagger} \mathbf{F A}$.
- Use a matrix-diagonalization method (Section 8.6) to find the eigenvalue and eigenvector matrices $\boldsymbol{\varepsilon}$ and $\mathbf{C}^{\prime}$ of $\mathbf{F}^{\prime}$.
- Calculate the coefficient matrix $\mathbf{C}=\mathbf{A C}^{\prime}$.
- Calculate an improved estimate of the density matrix from $\mathbf{C}$ using $\mathbf{P}^{*}=2 \mathbf{C} \mathbf{C}^{\dagger}$, which is the matrix form of (14.42) (Prob. 14.10c).
- Compare the improved $\mathbf{P}$ with the preceding estimate of $\mathbf{P}$. If all corresponding matrix elements differ by negligible amounts from each other, the calculation has converged and one uses the converged SCF wave function to calculate molecular properties. If the calculation has not converged, go back to step (5) to calculate an improved $\mathbf{F}$ matrix from the current $\mathbf{P}$ matrix and then do the succeeding steps.
One way to begin an SCF calculation is to initially estimate the Fock matrix elements by $F{r s} \approx H{r s}^{\text {core }}$, which amounts to neglecting the double sum in (14.41). This gives a very crude estimate. More commonly, SCF calculations get the initial estimate of the density matrix by doing a semiempirical calculation (Section 17.4) on the molecule. Semiempirical calculations are very fast. Still another possibility is to construct a guess for the $\mathbf{P}$ matrix by using the density matrices for the atoms composing the molecule. To find the equilibrium geometry of a molecule, one does a series of SCF calculations at many successive geometries (see Section 15.10). For the second and later SCF calculations of the series, one takes the initial guess of $\mathbf{P}$ as $\mathbf{P}$ for the SCF wave function of a nearby geometry.