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 H2 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.

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

\begin{equation}\hat{H} \Psi _k = E_k \Psi _k,\end{equation}
ĤΨk=EkΨk,
(1.1)

where Ψk and Ek are the eigenfunctions (wave functions) and eigenvalues (energies), respectively, of the Hamiltonian operator

$\hat{H}$
Ĥ⁠. However, the high-dimensionality and nature of
$\hat{H}$
Ĥ
make it difficult to solve (1.1) accurately.

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

\begin{equation}\hat{F} = \sum _{i=1}^n \hat{f}(\mathbf {r}_i).\end{equation}
F̂=i=1nf̂(ri).
(1.2)

The one-electron eigenfunctions, ψi, of

$\hat{f}$
f̂ are molecular orbitals (MOs) and a determinant of these is an antisymmetric eigenfunction of
$\hat{F}$
F̂
. Introducing a one-electron basis of functions, {ϕμ}, allows us to expand the MOs

\begin{equation}\psi _i = \sum _{\mu =1}^N C_{\mu i} \phi _{\mu }\end{equation}
ψi=μ=1NCμiϕμ
(1.3)

and convert the problem of finding eigenfunctions of

$\hat{f}$
f̂ into the generalized matrix eigenvalue problem

\begin{equation}\mathbf {F C} = \mathbf {S C \varepsilon },\end{equation}
FC=SCɛ,
(1.4)

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

\begin{equation}\mathbf {O} = (\mathbf {C}^\mathrm{old})^{\dagger } \mathbf {S} \mathbf {C}^\mathrm{new}\end{equation}
O=(C old )SC new
(1.5)

and calculating the projection pj = (∑iOij)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 H2 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 H2, 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) theory12 throughout. Most of the states that we will consider dissociate into dissimilar fragments, e.g., H(1s) + H(2pσ), and we therefore face the “symmetry dilemma” discussed by Löwdin13 many years ago. However, because our primary goal is to model the energetic characteristics of the various states of H2, 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.

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(1s) = –0.4999996, H(2s) = –0.1249831, H(2p) = –0.1249875, and He(1s2) = –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 H2 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

${\rm H}_2^+$
H2+ cation were satisfactory in all cases. Solutions for both short and long bond lengths were obtained in this way, and these were propagated for other bond lengths by moving in steps of 0.1 a.u., reading in the orbitals from the previous step, and using the MOM to maintain the desired state. By propagating in both directions, we were able to obtain both symmetric and broken-symmetry solutions for several of the states. The Allard–Kielkopf correlation diagram17 was used to check the dissociation products.

Fig. 1 shows the MO diagram for H2 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, re, dissociation energies, De, and harmonic frequencies, ωe, from Huber and Herzberg18 are listed in Table I, along with HF errors, Δ = HF – Exact. The potential energy curve for each state is shown in Fig. 2.

FIG. 1.

MO diagram for the H2 molecule.

FIG. 1.

MO diagram for the H2 molecule.

Close modal
Table I.

Exact values18 and HF deviations, Δ, for equilibrium bond lengths, re, dissociation energies, De, and harmonic vibrational frequencies, ωe, in the H2 molecule.

  re (Å)De (eV)ωe (cm−1)
HOMOStateExactΔExactΔExactΔ
g X |$^1\Sigma _g^+$|Σg+1 0.741 −0.008 4.748 −1.111 4401 90 
u b |$^3\Sigma _u^+$|Σu+3 … … … … … … 
  B |$^1\Sigma _u^+$|Σu+1 1.293 0.529 3.582 1.860 1358 −344 
g a |$^3\Sigma _g^+$|Σg+3 0.989 −0.004 3.057 −0.065 2665 38 
  EF |$^1\Sigma _g^+$|Σg+1 1.011 −0.015 2.543 0.059 2589 33 
u EF |$^1\Sigma _g^+$|Σg+1 2.315 −0.043 2.440 −0.967 1359 −159 
u e |$^3\Sigma _u^+$|Σu+3 1.107 −0.011 1.589 −0.042 2196 46 
  B|$^1\Sigma _u^+$|Σu+1 1.119 0.000 1.109 0.163 2039 58 
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 
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 
  re (Å)De (eV)ωe (cm−1)
HOMOStateExactΔExactΔExactΔ
g X |$^1\Sigma _g^+$|Σg+1 0.741 −0.008 4.748 −1.111 4401 90 
u b |$^3\Sigma _u^+$|Σu+3 … … … … … … 
  B |$^1\Sigma _u^+$|Σu+1 1.293 0.529 3.582 1.860 1358 −344 
g a |$^3\Sigma _g^+$|Σg+3 0.989 −0.004 3.057 −0.065 2665 38 
  EF |$^1\Sigma _g^+$|Σg+1 1.011 −0.015 2.543 0.059 2589 33 
u EF |$^1\Sigma _g^+$|Σg+1 2.315 −0.043 2.440 −0.967 1359 −159 
u e |$^3\Sigma _u^+$|Σu+3 1.107 −0.011 1.589 −0.042 2196 46 
  B|$^1\Sigma _u^+$|Σu+1 1.119 0.000 1.109 0.163 2039 58 
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 
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 
FIG. 2.

UHF potential energy curves for low-lying states of the H2 molecule.

FIG. 2.

UHF potential energy curves for low-lying states of the H2 molecule.

Close modal

The strengths and weaknesses of HF theory when describing the X

$^1\Sigma _g^+$
Σg+1 ground state of H2 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 De reflects the large correlation energy (Ec = −1.11 eV22,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

$^3\Sigma _u^+$
Σu+3⁠, dissociates correctly and agrees qualitatively with the exact curve. The exact curve has a very shallow minimum at re ≈ 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

$^1\Sigma _u^+$
Σu+1 state is poorly modeled by HF. All of the open-shell singlets are heavily spin-contaminated but, in this case, the contaminating state (b
$^3\Sigma _u^+$
Σu+3
) is repulsive and grossly distorts the energy curve of the B state. Like the spin-restricted HF solution for the ground state, it has half an electron of each spin on each center and is qualitatively incorrect in the dissociation limit. The correct dissociation products for this state are H(1s) + H(2pσ) which have a total energy of –0.625 a.u.

The deviations in the re and ωe value for the a

$^3\Sigma _g^+$
Σg+3 state are about half those for the ground state and the deviation in De is much smaller. Furthermore, its dissociation into H(1s) + H(2pσ) is qualitatively correct.

The singly excited E

$^1\Sigma _g^+$
Σg+1 and doubly excited F
$^1\Sigma _g^+$
Σg+1
solutions have the electron configurations
$1\sigma _g^1 2\sigma _g^1$
1σg12σg1
and
$1\sigma _u^2$
1σu2
, respectively. These interact strongly to form an EF
$^1\Sigma _g^+$
Σg+1
state23,30 with a double-minimum potential energy curve. The HF solution for the E
$1\sigma _g^1 2\sigma _g^1$
1σg12σg1
configuration yields re, De, and ωe values which agree well with the exact values at the first EF minimum. The HF solution for the F
$1\sigma _u^2$
1σu2
configuration is in poorer agreement with the exact curve at the second EF minimum (for example, De is underestimated by 1 eV) but this is not surprising because, as in the X
$1\sigma _g^2$
1σg2
state, the electrons share an orbital and are therefore strongly correlated. In HF theory, these two solutions cannot mix. However, it is known from elementary quantum mechanics that two same-symmetry solutions form an avoided crossing. Thus, these two HF solutions combine to form upper- and lower energy curves and we choose the lower one to model the EF state. The HF energy curve for the EF state has a cusp at 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 re, De, and ωe parameters for the e

$^3\Sigma _u^+$
Σu+3 state are well modeled by HF and the correct dissociation products, H(1s) + H(2s), are also obtained.

For the B

$^1\Sigma _u^+$
Σu+1 state, HF predicts re accurately but overestimates De by about 15% (0.006 a.u.). This De 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(1s) + H(2s), 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 De 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(1s) + H(2pπ) and both HF solutions dissociate correctly.

The I 1Πg and i 3Πg states have almost identical re, De, 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(1s) + H(2pπ) 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\Sigma _u^+$
Σu+1⁠, EF
$^1\Sigma _g^+$
Σg+1
, B
$^1\Sigma _u^+$
Σu+1
, and C 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(1s) + H(2pπ) 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 Hessian33,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.

FIG. 3.

Lowest eigenvalues λ of the HF MO rotation Hessian for the symmetric C 1Πu solution in the H2 molecule.

FIG. 3.

Lowest eigenvalues λ of the HF MO rotation Hessian for the symmetric C 1Πu solution in the H2 molecule.

Close modal
FIG. 4.

Lowest eigenvalues λ of the HF MO rotation Hessian for the symmetric I 1Πg solution in the H2 molecule.

FIG. 4.

Lowest eigenvalues λ of the HF MO rotation Hessian for the symmetric I 1Πg solution in the H2 molecule.

Close modal

Hartree-Fock theory (spin-unrestricted and allowing broken symmetry) is surprisingly effective for modeling the low-lying excited states of the H2 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.

1.
A.
Szabo
and
N. S.
Ostlund
,
Modern Quantum Chemistry
(
Dover
,
New York
,
1996
).
2.
D. J.
Thouless
,
The Quantum Mechanics of Many-Body Systems
(
Academic
,
New York
,
1961
).
3.
L.
González
,
D.
Escudero
, and
L.
Serrano-Andrés
,
ChemPhysChem
13
,
28
(
2012
).
4.
A. T. B.
Gilbert
,
N. A.
Besley
, and
P. M. W.
Gill
,
J. Phys. Chem. A
112
,
13164
(
2008
).
5.
E. R.
Davidson
and
L. Z.
Stenkamp
,
J. Chem. Phys.
41
,
656
(
1964
).
6.
A. L.
Thom
and
M.
Head-Gordon
,
Phys. Rev. Lett.
101
,
193001
(
2008
).
7.
J.
Cullen
,
M.
Krykunov
, and
T.
Ziegler
,
Chem. Phys.
391
,
11
(
2011
).
8.
B.
Peng
,
B. E.
Van Kuiken
,
F.
Ding
, and
X.
Li
,
J. Chem. Theory Comput.
9
,
3933
(
2013
).
9.
M.
Tassi
,
I.
Theophilou
, and
S.
Thanos
,
Int. J. Quantum Chem.
113
,
690
(
2013
).
10.
E. R.
Davidson
and
L. Z.
Stenkamp
,
Int. J. Quantum Chem.
10
,
21
(
1976
).
11.
F.
Jensen
,
Introduction to Computational Chemistry
(
Wiley
,
New York
,
1999
).
12.
J. A.
Pople
and
R. K.
Nesbet
,
J. Chem. Phys.
22
,
571
(
1954
).
13.
P.-O.
Löwdin
,
Rev. Mod. Phys.
35
,
496
(
1963
).
14.
F.
Jensen
,
J. Chem. Phys.
117
,
9234
(
2002
).
15.
K.
Szalewicz
and
H. J.
Monkhorst
,
J. Chem. Phys.
75
,
5785
(
1981
).
16.
A. I.
Krylov
and
P. M. W.
Gill
,
WIREs Comput. Mol. Sci.
3
,
317
(
2013
).
17.
N. F.
Allard
and
J. F.
Kielkopf
,
Astron. Astrophys.
493
,
1155
(
2009
).
18.
K. P.
Huber
and
G.
Herzberg
,
Molecular Spectra and Molecular Structure. IV. Constants of Diatomic Molecules
(
Van Nostrand Reinhold
,
New York
,
1979
).
19.
W.
Kołos
and
L.
Wolniewicz
,
J. Chem. Phys.
43
,
2429
(
1965
).
20.
W.
Kołos
and
L.
Wolniewicz
,
J. Chem. Phys.
45
,
509
(
1966
).
21.
W.
Kołos
and
L.
Wolniewicz
,
J. Chem. Phys.
48
,
3672
(
1968
).
22.
W.
Kołos
and
L.
Wolniewicz
,
J. Chem. Phys.
49
,
404
(
1968
).
23.
W.
Kołos
and
L.
Wolniewicz
,
J. Chem. Phys.
50
,
3228
(
1969
).
24.
W.
Kołos
and
J.
Rychlewski
,
J. Mol. Spectrosc.
62
,
109
(
1976
).
25.
W.
Kołos
and
J.
Rychlewski
,
J. Mol. Spectrosc.
66
,
428
(
1977
).
26.
L.
Wolniewicz
and
K.
Dressler
,
J. Chem. Phys.
88
,
3861
(
1988
).
27.
W.
Kołos
,
J. Mol. Spectrosc.
62
,
429
(
1976
).
28.
G.
Staszewska
and
L.
Wolniewicz
,
J. Mol. Spectrosc.
212
,
208
(
2002
).
29.
L.
Wolniewicz
and
G.
Staszewska
,
J. Mol. Spectrosc.
220
,
45
(
2003
).
30.
G.
Dieke
,
J. Mol. Spectrosc.
2
,
494
(
1958
).
31.
V. A.
Rassolov
,
M. A.
Ratner
, and
J. A.
Pople
,
J. Chem. Phys.
112
,
4014
(
2000
).
32.
Y.
Shao
,
M.
Head-Gordon
, and
A. I.
Krylov
,
J. Chem. Phys.
118
,
4807
(
2003
).
33.
J.
Čížek
and
J.
Paldus
,
J. Chem. Phys.
47
,
3976
(
1967
).
34.
P.-O.
Löwdin
,
J.
Calais
, and
J. M.
Calazans
,
Int. J. Quantum Chem.
20
,
1201
(
1981
).
35.
N. S.
Ostlund
,
J. Chem. Phys.
57
,
2994
(
1972
).
36.
W. D.
Allen
,
D. A.
Horner
,
R. L.
Dekock
,
R. B.
Remington
, and
H. F.
Schaefer
 III
,
Chem. Phys.
133
,
11
(
1989
).