Hartree-Fock (HF) theory is most often applied to study the electronic ground states of molecular systems. However, with the advent of numerical techniques for locating higher solutions of the self-consistent field equations, it is now possible to examine the extent to which such mean-field solutions are useful approximations to electronic excited states. In this Communication, we use the maximum overlap method to locate 11 low-energy solutions of the HF equation for the H_{2} molecule and we find that, with only one exception, these yield surprisingly accurate models for the low-lying excited states of this molecule. This finding suggests that the HF solutions could be useful first-order approximations for correlated excited state wavefunctions.

## I. INTRODUCTION

Quantum chemistry is concerned with the calculation of eigenstates of the Schrödinger equation

where Ψ_{k} and *E*_{k} are the eigenfunctions (wave functions) and eigenvalues (energies), respectively, of the Hamiltonian operator

One successful strategy is the Hartree-Fock (HF) method,^{1} which approximates Ψ_{k} by a single Slater determinant with optimized spin orbitals, leading to the approximation of the many-body Hamiltonian by a sum of one-electron operators

The one-electron eigenfunctions, ψ_{i}, of

_{μ}}, allows us to expand the MOs

and convert the problem of finding eigenfunctions of

where **C** is the MO coefficient matrix, **F** is the Fock matrix, **S** is the overlap matrix, and **ε** is the (diagonal) orbital energy matrix.

It is easy to show that the HF energy is a quartic polynomial in **C** and that solving (1.4) is equivalent to making this energy stationary with respect to **C**. According to the variational principle, the lowest solution is an upper bound to the true energy of the ground state of the system.^{2} However, a quartic polynomial in many variables usually has many stationary points and one may ask whether it is useful to interpret the higher solutions as mean-field approximations to excited states.

There are two objections to the proposal that excited states be modeled in this way. First, because it can be difficult to find the higher solutions of (1.4) and, second, because conventional wisdom holds that most excited states are intrinsically multi-configurational in nature.^{3}

The first objection has been partly addressed by the recent introduction of a number of interesting numerical techniques.^{4–9} In one of these, the Maximum Overlap Method (MOM),^{4} the occupation choice on each cycle of the self-consistent field (SCF) calculation is determined not by the aufbau principle but by maximizing overlap with the span of the occupied orbitals from the preceding cycle. This is implemented by computing the matrix

and calculating the projection *p*_{j} = (∑_{i}*O*_{ij})^{1/2} of each new orbital onto the space spanned by the old occupied orbitals. We then occupy the orbitals with the largest projections. Using the MOM, a SCF calculation can converge to higher energy solutions (stationary points) of the HF equations. This approach is general in that it can be applied to arbitrary molecules and is not limited to the lowest energy solutions of a particular symmetry.

The second objection requires further study and, in this Communication, we examine the extent to which single-determinant HF solutions are reasonable approximations to the ground state and ten singly excited states of the H_{2} molecule. We should emphasise from the outset that, within a single determinant model, the open-shell singlets are necessarily contaminated by the associated triplet. Attempts to remedy this have been made,^{10} but these approaches move beyond the single-determinant HF framework. Notwithstanding this spin contamination, we will refer to these contaminated solutions as “singlets” to differentiate them from the triplets.

Because the electrons no longer share an orbital in the excited states of H_{2}, dynamic correlation is less important than in the ground state and it is plausible that the excited states may, in fact, be modeled more accurately than the ground state. By a similar argument based on the Fermi hole which separates same-spin electrons in the HF approximation, it is possible that triplet states will be more accurately modeled than singlets.^{11} However, in states where static correlation plays a prominent role, it is possible that neither of these expectations may be met.

In order to describe homolytic bond-breaking, we have used spin-unrestricted Hartree-Fock (UHF) theory^{12} throughout. Most of the states that we will consider dissociate into dissimilar fragments, e.g., H(1*s*) + H(2*p*_{σ}), and we therefore face the “symmetry dilemma” discussed by Löwdin^{13} many years ago. However, because our primary goal is to model the energetic characteristics of the various states of H_{2}, whenever a solution splits into symmetric and broken-symmetry solutions, we choose the broken-symmetry branch, and whenever the solution splits into restricted and unrestricted solutions, we choose the unrestricted branch.

## II. METHOD AND RESULTS

First, we constructed a large basis set by adding two diffuse *p*-shells with exponents α_{1} = 0.0513930 and α_{2} = 0.0185657 to the aug-pc4 basis.^{14} These exponents were obtained by fitting the existing seven *p*-shell exponents to an exponential function, and extrapolating to obtain exponents for the eighth and ninth *p*-shells. Addition of a tenth *p*-shell had no effect on the results. The quality of the resulting basis set is indicated by the following HF energies (in atomic units): H(1*s*) = –0.4999996, H(2*s*) = –0.1249831, H(2*p*) = –0.1249875, and He(1*s*^{2}) = –2.8615949. The complete-basis limit for the last of these, which is the united-atom limit, is –2.8616799.^{15}

Using a modified version of the Q-Chem 4.0 package,^{16} we then used the MOM to compute the UHF energy curves of 11 states of the H_{2} molecule. The success of the MOM depends critically on the initial guesses for the occupied MOs and we found that the appropriate ground-state orbitals of the

^{17}was used to check the dissociation products.

Fig. 1 shows the MO diagram for H_{2} and each of the excited states considered has a single electron in the 1σ_{g} orbital and another in a higher orbital, as indicated in Table I. All of these states have been well studied, both experimentally and theoretically.^{18–30} Exact equilibrium bond lengths, *r*_{e}, dissociation energies, *D*_{e}, and harmonic frequencies, ω_{e}, from Huber and Herzberg^{18} are listed in Table I, along with HF errors, Δ = HF – Exact. The potential energy curve for each state is shown in Fig. 2.

. | . | r_{e} (Å)
. | D_{e} (eV)
. | ω_{e} (cm^{−1})
. | |||
---|---|---|---|---|---|---|---|

HOMO . | State . | Exact . | Δ . | Exact . | Δ . | Exact . | Δ . |

1σ_{g} | X |$^1\Sigma _g^+$|$\Sigma g+1$ | 0.741 | −0.008 | 4.748 | −1.111 | 4401 | 90 |

1σ_{u} | b |$^3\Sigma _u^+$|$\Sigma u+3$ | … | … | … | … | … | … |

B |$^1\Sigma _u^+$|$\Sigma u+1$ | 1.293 | 0.529 | 3.582 | 1.860 | 1358 | −344 | |

2σ_{g} | a |$^3\Sigma _g^+$|$\Sigma g+3$ | 0.989 | −0.004 | 3.057 | −0.065 | 2665 | 38 |

EF |$^1\Sigma _g^+$|$\Sigma g+1$ | 1.011 | −0.015 | 2.543 | 0.059 | 2589 | 33 | |

1σ_{u} | EF |$^1\Sigma _g^+$|$\Sigma g+1$ | 2.315 | −0.043 | 2.440 | −0.967 | 1359 | −159 |

2σ_{u} | e |$^3\Sigma _u^+$|$\Sigma u+3$ | 1.107 | −0.011 | 1.589 | −0.042 | 2196 | 46 |

B^{′} |$^1\Sigma _u^+$|$\Sigma u+1$ | 1.119 | 0.000 | 1.109 | 0.163 | 2039 | 58 | |

1π_{u} | c ^{3}Π_{u} | 1.037 | −0.006 | 3.069 | −0.127 | 2467 | 39 |

C ^{1}Π_{u} | 1.033 | −0.006 | 2.542 | 0.074 | 2444 | 62 | |

1π_{g} | i ^{3}Π_{g} | 1.070 | −0.009 | 0.925 | −0.067 | 2253 | 31 |

I ^{1}Π_{g} | 1.069 | −0.009 | 0.924 | −0.072 | 2259 | 44 |

. | . | r_{e} (Å)
. | D_{e} (eV)
. | ω_{e} (cm^{−1})
. | |||
---|---|---|---|---|---|---|---|

HOMO . | State . | Exact . | Δ . | Exact . | Δ . | Exact . | Δ . |

1σ_{g} | X |$^1\Sigma _g^+$|$\Sigma g+1$ | 0.741 | −0.008 | 4.748 | −1.111 | 4401 | 90 |

1σ_{u} | b |$^3\Sigma _u^+$|$\Sigma u+3$ | … | … | … | … | … | … |

B |$^1\Sigma _u^+$|$\Sigma u+1$ | 1.293 | 0.529 | 3.582 | 1.860 | 1358 | −344 | |

2σ_{g} | a |$^3\Sigma _g^+$|$\Sigma g+3$ | 0.989 | −0.004 | 3.057 | −0.065 | 2665 | 38 |

EF |$^1\Sigma _g^+$|$\Sigma g+1$ | 1.011 | −0.015 | 2.543 | 0.059 | 2589 | 33 | |

1σ_{u} | EF |$^1\Sigma _g^+$|$\Sigma g+1$ | 2.315 | −0.043 | 2.440 | −0.967 | 1359 | −159 |

2σ_{u} | e |$^3\Sigma _u^+$|$\Sigma u+3$ | 1.107 | −0.011 | 1.589 | −0.042 | 2196 | 46 |

B^{′} |$^1\Sigma _u^+$|$\Sigma u+1$ | 1.119 | 0.000 | 1.109 | 0.163 | 2039 | 58 | |

1π_{u} | c ^{3}Π_{u} | 1.037 | −0.006 | 3.069 | −0.127 | 2467 | 39 |

C ^{1}Π_{u} | 1.033 | −0.006 | 2.542 | 0.074 | 2444 | 62 | |

1π_{g} | i ^{3}Π_{g} | 1.070 | −0.009 | 0.925 | −0.067 | 2253 | 31 |

I ^{1}Π_{g} | 1.069 | −0.009 | 0.924 | −0.072 | 2259 | 44 |

## III. DISCUSSION

The strengths and weaknesses of HF theory when describing the X

_{2}are well-known. HF underestimates the equilibrium bond length by 0.008 Å, underestimates the dissociation energy by 1.111 eV, and overestimates the harmonic frequency by 90 cm

^{−1}. The large error in

*D*

_{e}reflects the large correlation energy (E

_{c}= −1.11 eV

^{22,31}) at the equilibrium geometry. (The correlation energy of the dissociated atoms is exactly zero.) These results provide a benchmark for determining how well HF models the excited states.

The lowest energy triplet state, b

*r*

_{e}≈ 7.8 a.u.

^{19}due to dispersion, but HF is incapable of modeling this and the HF solution is therefore purely repulsive.

The B

*s*) + H(2

*p*

_{σ}) which have a total energy of –0.625 a.u.

The deviations in the *r*_{e} and ω_{e} value for the a

*D*

_{e}is much smaller. Furthermore, its dissociation into H(1

*s*) + H(2

*p*

_{σ}) is qualitatively correct.

The singly excited E

^{23,30}with a double-minimum potential energy curve. The HF solution for the E

*r*

_{e},

*D*

_{e}, and ω

_{e}values which agree well with the exact values at the first EF minimum. The HF solution for the F

*D*

_{e}is underestimated by 1 eV) but this is not surprising because, as in the X

*R*= 3.52 a.u. (the maximum on the true EF curve lies at

*R*= 3.13 a.u.). An analogous cusp arises in the HF energy curve for the torsional rotation in ethene.

^{32}

The *r*_{e}, *D*_{e}, and ω_{e} parameters for the e

*s*) + H(2

*s*), are also obtained.

For the B^{′}

*r*

_{e}accurately but overestimates

*D*

_{e}by about 15% (0.006 a.u.). This

*D*

_{e}error is among the largest in our study, but is still much smaller than that for the ground state. Although the dissociation limit for this state, H(1

*s*) + H(2

*s*), can be accurately modeled using a broken-symmetry UHF wave function, we were unable to locate a HF solution connecting these fragments to the equilibrium structure. Thus, the HF curve for this state shown in Fig. 2 dissociates to an unphysical solution with half an alpha electron and half a beta electron on each centre.

The HF descriptions of the c ^{3}Π_{u} and C ^{1}Π_{u} states are similar. The deviation in both bond lengths is 0.006 Å, which is similar to, but slightly better than, the ground state deviation of 0.008 Å. The *D*_{e} and ω_{e} errors are also smaller than those of the ground state, with the triplet state being reproduced slightly better. The dissociation products for both states are H(1*s*) + H(2*p*_{π}) and both HF solutions dissociate correctly.

The I ^{1}Π_{g} and i ^{3}Π_{g} states have almost identical *r*_{e}, *D*_{e}, and ω_{e} values. This similarity is not unexpected because the electron in the diffuse 1π_{g} orbital rarely encounters the electron in the compact 1σ_{g} orbital and their spins are therefore largely irrelevant. HF reproduces this, giving almost the same errors for the two states. However, neither solution dissociates correctly to the H(1*s*) + H(2*p*_{π}) fragments. These states have small well depths and the errors of roughly 0.07 eV, although comparable to those for other states, are relatively larger.

It is interesting to note that HF theory *overestimates* the dissociation energies of four of the singlet states that we studied, viz. B

^{′}

^{1}Π

_{u}. This implies that each of these states has a

*positive*correlation energy and, therefore, that the higher solutions of the HF equation are non-variational.

We also note that, when a HF solution breaks symmetry as the bond is stretched, this always occurs beyond the equilibrium bond length of the target solution, thereby ensuring that the data in Table I are unique. To understand the symmetry-breaking phenomenon in greater detail, we have analyzed the C ^{1}Π_{u} and I ^{1}Π_{g} states. Both states dissociate to H(1*s*) + H(2*p*_{π}) but only the former breaks symmetry and can dissociate to the correct energy. The C ^{1}Π_{u} solution bifurcates at 4.3 a.u. into a higher energy symmetric and a lower energy broken-symmetry solution. No such splitting occurs for the I ^{1}Π_{g} solution.

Fig. 3 shows the lowest eigenvalues of the MO rotation Hessian^{33,34} for the symmetric C ^{1}Π_{u} as a function of bond length. The bifurcation arises at the bond length where one of the eigenvalues becomes negative, indicating the existence of a new, lower energy solution.^{35,36} Closer inspection of the corresponding eigenvector reveals that the rotations involved include mixing of the occupied π_{u} orbital with the virtual π_{g} orbitals, thus breaking the symmetry of the solution. The eigenvalues for the symmetric I ^{1}Π_{g} solution (Fig. 4) do not exhibit the same behaviour, which explains the absence of a broken symmetry solution for this state. In each case, the zero eigenvalue indicates that each of the Π solutions is degenerate.

## IV. CONCLUSIONS

Hartree-Fock theory (spin-unrestricted and allowing broken symmetry) is surprisingly effective for modeling the low-lying excited states of the H_{2} molecule. We find that it provides useful first-order descriptions of ten of the eleven states considered and, in fact, with the single exception of the B ^{1}Σ_{u} state, it is more accurate for the excited states (near their equilibrium bond lengths) than it is for the ground state. We also find, as anticipated, that the triplet states are usually more accurately modeled than the singlets.

The accuracy of the results in Table I strongly suggests that there is a correspondence between solutions of the HF equation and solutions of the Schrödinger equation. It appears likely that the higher HF solutions may often (though not always) be a satisfactory first-order approximation for excited states, which can subsequently be improved, if desired, using standard correlation methods. We are investigating this, for both singly- and multiply excited states, and will present our findings elsewhere.