The Hartree–Fock problem provides the conceptual and mathematical underpinning of a large portion of quantum chemistry. As efforts in quantum technology aim to enhance computational chemistry algorithms, the Hartree–Fock method, central to many other numerical approaches, is a natural target for quantum enhanced algorithms. While quantum computers and quantum simulation offer many prospects for the future of modern chemistry, the non-deterministic polynomial-complete Hartree–Fock problem is not a likely candidate. We highlight this fact from a number of perspectives including computational complexity, practical examples, and the full characterization of energy landscapes for simple systems.

The present study aims to highlight the difficulty in conducting optimization problems in the context of electronic structure, with and without quantum resources. Specifically, we will focus on the Hartree–Fock problem, an optimization problem using the mean-field approximation. The Hartree–Fock problem1,2 provides the mathematical setting for molecular orbitals widely used in chemistry and beyond. The ubiquitous self-consistent field (SCF) methodology used to solve the Hartree–Fock problem is also applied to most implementations of density functional theory based on Kohn–Sham theory.3 While the solution to an instance of the Hartree–Fock problem is often insufficient for many applications, it often serves as the reference state for post-Hartree–Fock methods. Widely used post-Hartree–Fock methods include the coupled cluster ansatz,4–9 Møller–Plesset perturbation theory,10–13 equation of motion,14–20 multi-reference configuration interaction,21–23 and many more. Here, we forgo improving the Hartree–Fock ansatz and instead ask how difficult it is to find the true Hartree–Fock global minimum and its importance.

Many instances of the Hartree–Fock problem can be solved quickly using heuristic approaches. In practice, conventional algorithms for Hartree–Fock scale cubically with the number of basis functions. However, this cost only reflects the correct scaling if the number of iterations is bounded by a constant or if local minima are acceptable in place of a global solution. Linear scaling methods24,25 avoid the diagonalization of Hessian entirely and rely on localized properties of the system. This assumption may not be true generally for every case. Iterative procedures, regardless of cost per step, are prone to convergence issues. This article highlights the above-mentioned properties in the context of standard numerical and hybrid-quantum approaches to the Hartree–Fock problem.

Typical approaches to improve SCF convergence use direct inversion of the iterative subspace26 (DIIS), level shifting,27 quadratically convergent Newton–Raphson techniques,28 variational coupled cluster methods,29 or varying fractional occupation numbers,30 among many other approaches. The success of any above-mentioned methods varies at each instance and depends on the initial parameters chosen, e.g., the size of the iterative subspace, and often work well in combination, leading to attempts to build black-box SCF procedures.31,32

Unfortunately, these methods cannot work in all cases without violating fundamental assumptions on the theory of computation. In previous work, arbitrary spin-glasses have been mapped to instances of the Hartree–Fock problem.33 In addition to other works,34,35 this shows that the Hartree–Fock problem is difficult in the worst case setting and non-deterministic polynomial (NP)-complete.36 The complexity class of non-deterministic polynomial (NP) time problems is the set of problems37 that can be solved efficiently with a hint.38 It is possible that every problem that can be solved efficiently with a hint can also be efficiently solved without the hint. However, all signs of practical experience suggest that problems in the NP class do not admit efficient black-box solutions. Thus, due to the NP-completeness of the Hartree–Fock problem, it is unlikely that any classical algorithm can solve all instances in a time proportional to a polynomial of the input size. Note that because it is in the NP complexity class, the Hartree–Fock problem can be solved efficiently in polynomial time with a sufficiently strong hint, e.g., as gleamed from experience or luck. It has been known since the late 1990s that quantum computers can promise no more than a quadratic speedup over their classical counterparts on such NP-complete problems.39 

It has long been known that any classical algorithm can be simulated using quantum hardware. With the heavy reliance of the variational quantum eigensolvers on classical optimization routines, it is marginal how the quantum approach differs from the conventional approaches. In the recent variational quantum eigensolver (VQE)40,41 study, the optimization strategy used was an augmented Hessian approach42 rather than the typical DIIS.26 However, both of these optimization strategies can be employed by conventional computers to optimize the orbitals.

In this article, we consider how the use of quantum hardware enhances the ability of chemists to solve instances of the Hartree–Fock problem. Quantum hardware has many prospects for applications to physical and chemical simulations;43 however, the Hartree–Fock problem is not likely to admit drastic advances using quantum computers.

The present article is inspired by a recent study of the Google group and collaborators40 using a variational quantum–classical hybrid formulation of the Hartree–Fock problem. There, the authors pointed out that the primary purpose of considering the Hartree–Fock problem was to benchmark their quantum device. Here, we highlight obstructions to its use as a general purpose replacement for standard SCF solvers. We do so by first introducing the Hartree–Fock problem and characterizations of the solution landscape for an instance of the Hartree–Fock problem. We then give examples that are (1) simple, (2) small, and (3) well-motivated instances of the Hartree–Fock problem that present convergence problems for black-box approaches.

The Hartree–Fock ansatz is important for conventional quantum chemistry and has been thoroughly developed over the past 100 years.44 

In this subsection, we wish to succinctly sketch the key ideas of Hartree–Fock theory. First, we give a standard presentation in terms of orbital rotations, and then, we reformulate the problem using bond density matrices as a useful representation of the problem.

The Hartree–Fock optimization problem refers to the collection of instances each specific by a given electronic Hamiltonian characterized by one- and two-electron integrals ({hij}, {hijkl}). Once an instance is specified, the Hartree–Fock problem can be stated succinctly as minimizing the electronic energy in the space of single Fock states,

EHF=minΨF1Ψ|ijM2hijaiaj+ijklM4hijklaiajakal|Ψ.
(1)

The space of single Fock states F1 is the set of all rank 1 [vs rank MN in the general case] N-electron Fock states. Here and throughout, N is the number of electrons and M is the number of spatial basis functions. Each single Fock state is of the form

Ψ=b1b2bN|Ω,
(2)

where we have [bi,bj]+=bibj+bjbi=δij and [bj,bk]+=0. The vacuum state with no particles is denoted by |Ω⟩ and is characterized by aj|Ω⟩ = 0 for all j. Alternative formulations of the Hartree–Fock problem can be done using Slater determinants.45 

After the restriction to this ansatz, the variational space is characterized by rotations of the orbitals from an initial set of orbitals (e.g., the atomic orbitals). If the initial set corresponds to operators {aj:[ai,aj]+=dxϕi(x)ϕj*(x)=Sij,[ai,aj]+=0}, then

bi=jWjiaj.
(3)

Here, we require that WW = 1 when S = 1 and W̃SW̃=1 in the general case. For applications to molecular physics, single-electron spin-orbitals ϕi(x) corresponding to the fermionic operators aj and aj are often contracted Gaussians, e.g., STO-3G.46 If the orbitals in Eq. (3) yield the minimal energy, then they are called the molecular orbitals.

The optimization problem stated in Eq. (1) can be solved using any of a variety of techniques. The SCF approach to instances of the Hartree–Fock problem is a direct application of the steepest descent approach, but other methods such as semi-definite programing can be used.47 In practice, the SCF method to solve Eq. (1) uses an effective potential term that takes into account the averaged two-body interaction (mean-field). In this approach, the N-body problem is reduced to a non-linear single particle problem. At each iteration of the simplest implementation, the Fock matrix, F, is formed as a function of a previous bond density matrix (Pprev), also known as one-particle reduced density matrix or the two-point correlation function,

Pprev=Ψ|(ajak+ajak)|ΨjkM=WprevηWprev,
(4)

with η as the orbital occupancies written as a diagonal matrix. The new transformation matrix, W, is determined using the gradient of Eq. (1) with respect to the bond density matrix P. The new coefficient matrix is used to form a new bond density matrix. At each iteration, the (real-valued) bond density matrix satisfies the following three properties:48 

P=PT,
(5)
Tr(P)=N/2,
(6)
P=P2.
(7)

Next, we consider transformation between bond density matrices. A convenient parameterization of the set of bond density matrices that avoids redundancy is

P(A)=WAP0WA=eAblockP0eAblock,
(8)

with Ablock = P0A(1P0) + (1P0)AP0 for arbitrary skew-symmetric A = −AT. In Eq. (8), WA = exp(−Ablock) is a change of basis as in Eq. (3).

The objective for the optimization problem in Eq. (1) is also given by the energy matrix functional49 

EHF=minPE[P]=minP2Tr[hP]+Tr[GP].
(9)

Here, the mean-field, G, is a function of the bond density matrix

G=G[P]μν=κ,λMhκμνλAPκλ,
(10)

with the antisymmetrized two-electron integral over spatial degrees of freedom defined by hμκλνA=hμκλνhμκνλ. Note that Eq. (9) follows directly from evaluating the energy’s expectation value in Eq. (1) using, e.g., Slater–Condon rules.50 The matrix derivative of Eq. (9) gives the Fock operator

Fij=2(hij+Gij).
(11)

By rotating the charge density matrix with all possible rotations, we are able to do brute force exploration of the whole space for some small examples. Below, we expand Eq. (8) about P0 = Pcore where the N lowest eigenmodes of Hcore=hijaiaj are occupied.

It is important to note that the global minimum of the functional equation (9) has significance. If one wishes to obtain approximations to the excited states of the Hamiltonian, Eq. (9) should not be used; instead, excited state functionals51 can be derived using excited state mean-field theory.52 

Before turning to examples, we will introduce the quantum ansatz for the Hartree–Fock method.

The quantum circuit for creating the Hartree–Fock ansatz state53 was applied in the context of VQE.40 The description of the ansatz circuit will be aided through the use of QR decomposition.54 

The Givens rotations provide a useful canonical characterization of an arbitrary orthogonal matrix, W. The QR decomposition of a real n × n orthogonal matrix W can be done using T = n(n − 1)/2 Givens rotations such that

W=G1G2G4GTD.
(12)

When W has a determinant of 1, D is just the identity matrix. Each Givens rotation, Gi, is of the form Gi = g(a, b, θ) with gkk = 1, unless k is either a or b when instead gkk = cos(θ). All off-diagonal elements are zero, except gab = −gba = − sin(θ).

Applications of the Givens decomposition to fermionic orbital rotations has been worked out elsewhere,53,55 resulting in a quantum circuit that is able to prepare arbitrary Slater determinants following the parameters of the QR decomposition. By ordering the QR decomposition appropriately, a fermionic swap network can be used to rotate each pair of orbitals using the appropriate Givens rotation parameters. This results in an efficient state preparation circuit of the form depicted in Fig. 1. The full compilation down to gates including hardware optimization is given elsewhere.40,41

FIG. 1.

Quantum circuit for the two fermion and three molecular orbital (for example, H3+ in STO-3G basis).

FIG. 1.

Quantum circuit for the two fermion and three molecular orbital (for example, H3+ in STO-3G basis).

Close modal

Our characterization of the fermionic space in Eq. (8) gives us a set of parameters, Θ, that also characterizes the mixing between pairs of orbitals. The resulting orthogonal transformation W(Θ) is then given to the QR decomposition and forwarded to the quantum circuit construction.

All calculations of the molecular system are done in the STO-3G basis.46 Energies are reported in Hartrees, angles of rotation are in radians, and bond lengths are in angstroms.

The Hartree–Fock energy surfaces (HESs) were computed using PySCF.56,57 In this paper, we only consider Restricted Hartree–Fock (RHF) solutions where the alpha and beta spatial orbitals are restricted to be identical. The quantum optimization routines were that of OpenFermion-Cirq,41 and we only modify the initial state routines and the input molecular data.58 Internally, the optimization of the quantum circuit is performed using a conjugate gradient as implemented in SciPy.59 

We consider H2, H3+ as minimal basis model systems whose Hartree–Fock instances we can completely characterize. We begin with the H2 example.

When considering H2 in the minimal basis, there is only a single orbital mixing parameter. In Fig. 2, we have plotted the 1D HES surface as a function of bond length for H2. The number of minima in HES(θ) changes with bond length. Before a bond length of ∼1.2 Å, there is only a single minimum. However, at larger bond lengths, an additional minimum begins to appear.

FIG. 2.

H2 HES(θ, rHH). Each fixed value of the nuclear separation, rHH, generates a Hartree–Fock instance characterized by a single orbital rotation parameter, θ. Note the appearance of a second HES minimum at a higher energy around r ≈ 1.2 Å.

FIG. 2.

H2 HES(θ, rHH). Each fixed value of the nuclear separation, rHH, generates a Hartree–Fock instance characterized by a single orbital rotation parameter, θ. Note the appearance of a second HES minimum at a higher energy around r ≈ 1.2 Å.

Close modal

We continue with our two electron examples with the iso-electronic H3+. Now, instead of a 1D HES, we now have two parameters that mix the one occupied orbital with the two virtual orbitals. We plot the HES in Fig. 3 for a linear configuration with hydrogen atoms separated by 2.5 Å. There are three minima for HES(θ1, θ2). In Fig. 4, we give the HES of H3+ at 4.36 Å where there are several minima with the same globally optimal value.

FIG. 3.

H3+ HES(θ1, θ2), showing three different minima (global minimum—white circles, second minimum—orange triangles, and third minimum—green crosses) at rHH = 2.5 Å. Here, the surface is expanded about P0 = Pcore.

FIG. 3.

H3+ HES(θ1, θ2), showing three different minima (global minimum—white circles, second minimum—orange triangles, and third minimum—green crosses) at rHH = 2.5 Å. Here, the surface is expanded about P0 = Pcore.

Close modal
FIG. 4.

The H3+ HES(θ1, θ2) surface for fixed bond length at rHH = 4.36 Å. Here, the surface is expanded about P0 = Pcore.

FIG. 4.

The H3+ HES(θ1, θ2) surface for fixed bond length at rHH = 4.36 Å. Here, the surface is expanded about P0 = Pcore.

Close modal

In both the cases of H3+ and H2, there are single occupied (spatial) orbital and m virtual orbitals. For m = 2, this leads to an Ablock generator of the form

Ablock=0θ1θ2θ100θ200.
(13)

The eigenvalues of this matrix are λ±={0,±iθ12+θ22}. Since the matrix exponential of Ablock merely exponentiates the eigenvalues, when λ± = , the rotation acts trivially on the density matrix. This underlies periodicity to the plots seen in Figs. 3 and 4.

We can explain the periodicity in terms of this invariant by converting to polar coordinates where θ1 = R cos ϕ and θ2 = R sin ϕ. Now, the nontrivial eigenvalue is λ± = ±iR, and we can express the periodicity of the plots as HES(θ1, θ2) = HES(R, ϕ) = HES(R + , ϕ), with n being an integer.

There is a nice generalization of this fact. For a single spatial orbital that is doubly occupied with two electrons and m virtual orbitals, the generalization of Eq. (13) is

Ablock=0θ1θ2θ3θmθ10000θ20000θ30000θm0000.
(14)

It is straightforward to calculate that the eigenvalues of this matrix are zero except λ±=±iθ12+θ22+θ32++θm2. Following the same argument as in the m = 2 case above, we can say that

HES(Θ)=HES(R,Φ)=HES(R+nπ,Φ),
(15)

where R=θ12++θm2.

The range of the search space for each θj is [0, 2π). Straightforward combination of parameter results is a parameter space for bond density matrices that is a hyper-cube in m dimensions with equal side lengths of 2π. However, due to the constraint of Eq. (15), the minimal search space is now confined to the interior of a hyper-sphere of radius π (also in m dimension). Therefore, the ratio of the minimal required search space to the full search goes to zero as m tends to infinity. This is a well-known consequence of the vanishing ratio of the volume of a hyper-sphere to the volume of the corresponding hyper-cube.60,61

We used the quantum algorithm outlined in Ref. 40 for obtaining RHF solutions for four examples. Depending on the initial guess, it may converge to local rather global solutions.

The different initial guesses were generated using the Givens rotations corresponding to different minima in Figs. 2 and 3, respectively. The results for H2, converging to two different minima on the quantum simulator, are shown in Fig. 5. Similarly, results for H3+, converging to two different minima on the quantum simulator, are shown in Fig. 6. In Fig. 5, the values of λ = 1 and λ = 1.2 were chosen as states far enough in parameter space to have energy sufficiently large. If we select states with small λ, the convergence to the minimum is highly likely so long as the system does not climb uphill in energy since the initial state has energy less than all minima except the global minimum.

FIG. 5.

Convergence of the HF quantum optimization to different minima for the Hartree–Fock instance of H2 at rHH = 1.5 Å in a minimal basis. Depicted are two initial conditions obtained by adding offsets of differing strengths, λ, to the global ground state. The value used in the OpenFermion implementation was λ = 0.1.

FIG. 5.

Convergence of the HF quantum optimization to different minima for the Hartree–Fock instance of H2 at rHH = 1.5 Å in a minimal basis. Depicted are two initial conditions obtained by adding offsets of differing strengths, λ, to the global ground state. The value used in the OpenFermion implementation was λ = 0.1.

Close modal
FIG. 6.

Convergence of the HF quantum optimization to different minima for the Hartree–Fock instance of H3+ at rHH = 2.5 Å in a minimal basis. Depicted are two initial conditions obtained by adding offsets of differing strengths, λ, to the global ground state.

FIG. 6.

Convergence of the HF quantum optimization to different minima for the Hartree–Fock instance of H3+ at rHH = 2.5 Å in a minimal basis. Depicted are two initial conditions obtained by adding offsets of differing strengths, λ, to the global ground state.

Close modal

As final examples, we choose diatomic carbon and its cation. We also consider C2 and C22+ as instances that are commonly known to confound solvers due to the appearance of saddle points in the optimization landscape. To illustrate the complications of convergence, we modified the initial parameters following the method used in Ref. 40. Namely, we begin with the solution provided by the classical SCF solver, perturb from those solutions, and observe if the quantum algorithm still converges to the correct minimum.

Since the Hartree–Fock problem is an optimization problem, quantum computing via a modified Grover search39,62 can be used to find the global solution quadratically faster than the classical brute force approaches. However, in both conventional and quantum solvers, local searches are employed whereby no guarantees on finding the globally optimal solutions are given.

This example allows us to highlight the importance of using the information of Hessian to avoid saddle points during SCF optimization in either quantum or classical methods. In Figs. 7 and 8, convergence results are shown for C2 and C22+, respectively. In each of the plots, we have plotted the performance of the quantum circuit optimization routine and the classical optimization routines as implemented in PySCF. At an inter-nuclear separation of 1.5 Å, there is a saddle point appearing in the potential energy landscape.

FIG. 7.

SCF convergence for different optimizers for the minimal basis C2 at rCC = 1.5 Å. Plotted are the optimization trajectories for the quantum circuit optimizer and PySCF with and without an extra step. The energy of the saddle point and of the global minimum of the HES is also shown. The PySCF guess is Pcore.

FIG. 7.

SCF convergence for different optimizers for the minimal basis C2 at rCC = 1.5 Å. Plotted are the optimization trajectories for the quantum circuit optimizer and PySCF with and without an extra step. The energy of the saddle point and of the global minimum of the HES is also shown. The PySCF guess is Pcore.

Close modal
FIG. 8.

Optimization to local minima for C22+ on the quantum simulator while avoiding convergence to the saddle point at rCC = 1.0 Å. The PySCF guess is Pcore.

FIG. 8.

Optimization to local minima for C22+ on the quantum simulator while avoiding convergence to the saddle point at rCC = 1.0 Å. The PySCF guess is Pcore.

Close modal

The use of the orbital Hessian helps solvers avoid saddle points where the gradient may vanish at non-optimal values. The solver uses Hessian, which provides a notion of the curvature of the landscape. This allows the solver to avoid saddle points. The augmented Hessian conjugate-gradient method42 was used for the quantum optimization of the Hartree–Fock circuit. This allows the solver to avoid convergence to saddle points, as illustrated in Figs. 7 and 8. Many of the quantum chemistry do not check the RHF solutions by default. This makes them prone to failure by convergence to a saddle point rather than the RHF solution.

In this study, we showed that the convergence to a local minimum or global minimum is a function of the initial guess. The proposed algorithm on the quantum simulator carries this feature from the classical algorithm. There is no a priori guarantee that it will find the global solution.

Since the Hartree–Fock problem is an optimization problem, quantum computing via a modified Grover search39,62 can be used to find the global solution quadratically faster than the classical brute force approaches. However, in both conventional and quantum solvers, local searches are employed whereby no guarantees on finding the globally optimal solutions are given.

In analyzing the Hartree–Fock functional for H2, H3+, C2, and C22+, we note that the number of critical points in the solution space63 changes as the nuclear separation changes. This also implies that there are additional difficulties in applying Hartree–Fock for nuclear dynamics or other non-equilibrium configurations. Most solvers for the Hartree–Fock problem do not explore the entire space and usually only choose a single starting point (rather than multiple starting points). While this is not often an issue, it can cause serious pathology when used in post-Hartree–Fock methods, e.g., the coupled cluster method. This will be true for both on conventional computers and in its unitary formulation for quantum computers.

The recent application of quantum technology to the Hartree–Fock problem may serve as a hardware benchmark but is unlikely to have a dramatic impact on the practical approaches to this problem.

While the application to the Hartree–Fock problem is not likely to change the workhorse routines used on conventional computers, there are still interesting use cases for the results from Ref. 40. The projection methods and purity extrapolation presented in Ref. 40 will still be useful.

We should point out that the obstructions to efficiently solving the Hartree–Fock problem using quantum resources is not directly related to the concentration of measure phenomena, leading to barren plateaus discussed in the context of VQE ansatz initialization.64 It would be interesting to combine the insights of computational complexity with similar arguments for the Hartree–Fock ansatz.

When considering application areas of quantum computers, it is far more likely to make major breakthroughs when considering time-dependent phenomena. For example, the quantum–classical hybrid algorithm for obtaining the Kohn–Sham potential of the time-dependent density functional65–67 also requires measuring the bond density matrix.

This article highlights the wide body of knowledge on the SCF method and its difficulty and shows that the lack of verification of the solution is material in both the quantum and conventional computing domains.

The authors thank Z. Zimboras and N. Anand for many useful discussions. J.D.W. acknowledges support from the Department of Energy under Grant No. DE-SC0019374 “Quantum Chemistry for Quantum Computers.” Additional support for J.D.W. comes from the NSF (Grant Nos. PHYS-1820747 and EPSCoR-1921199) and from the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research under programs Quantum Computing Application Teams and Accelerated Research for Quantum Computing. S.G. would like to thank USC for the Dornsife Graduate School Fellowship.

The data that support the findings of this study are available from the authors upon reasonable request.

1.
D. R.
Hartree
, “
The wave mechanics of an atom with a non-Coulomb central field. Part I: Theory and methods
,”
Proc. Cambridge Philos. Soc.
24
,
89
312
(
1928
).
2.
V.
Fock
, “
Näherungsmethode zur lösung des quantenmechanischen Mehrkörperproblems
,”
Z. Phys.
61
,
723
734
(
1930
).
3.
W.
Kohn
and
L. J.
Sham
, “
Self-consistent equations including exchange and correlation effects
,”
Phys. Rev.
140
,
1133
(
1965
).
4.
F.
Coester
, “
Bound states of a many-particle system
,”
Nucl. Phys.
7
,
421
424
(
1958
).
5.
F.
Coester
and
H.
Kümmel
, “
Short-range correlations in nuclear wave functions
,”
Nucl. Phys.
17
,
477
485
(
1960
).
6.
J.
Čížek
, “
On the correlation problem in atomic and molecular systems. Calculation of wavefunction components in Ursell-type expansion using quantum-field theoretical methods
,”
J. Chem. Phys.
45
,
4256
4266
(
1966
).
7.
R. J.
Bartlett
, “
Many-body perturbation theory and coupled cluster theory for electron correlation in molecules
,”
Annu. Rev. Phys. Chem.
32
,
359
401
(
1981
).
8.
J.
Paldus
and
X.
Li
, “
A critical assessment of coupled cluster method in quantum chemistry
,”
Adv. Chem. Phys.
110
,
1
175
(
1999
).
9.
I.
Shavitt
and
R. J.
Bartlett
,
Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory
(
Cambridge University Press
,
2009
).
10.
W.
Kutzelnigg
, “
How many-body perturbation theory (MBPT) has changed quantum chemistry
,”
Int. J. Quantum Chem.
109
,
3858
3884
(
2009
).
11.
I.
Lindgren
and
J.
Morrison
,
Atomic Many-Body Theory
, Springer Series in Chemical Physics Vol. 13 (
Springer Berlin
,
1982
).
12.
J. S.
Binkley
and
J. A.
Pople
, “
Møller–Plesset theory for atomic ground state energies
,”
Int. J. Quantum Chem.
9
,
229
236
(
1975
).
13.
D.
Cremer
, “
Møller–Plesset perturbation theory: From small molecule methods to methods for thousands of atoms
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
1
,
509
530
(
2011
).
14.
D. J.
Rowe
, “
Equations-of-motion method and the extended shell model
,”
Rev. Mod. Phys.
40
,
153
(
1968
).
15.
K.
Emrich
, “
An extension of the coupled cluster formalism to excited states (I)
,”
Nucl. Phys. A
351
,
379
396
(
1981
).
16.
J.
Geertsen
,
M.
Rittby
, and
R. J.
Bartlett
, “
The equation-of-motion coupled-cluster method: Excitation energies of Be and Co
,”
Chem. Phys. Lett.
164
,
57
62
(
1989
).
17.
J. F.
Stanton
and
R. J.
Bartlett
, “
The equation of motion coupled-cluster method: A systematic biorthogonal approach to molecular excitation energies, transition probabilities, and excited state properties
,”
J. Chem. Phys.
98
,
7029
7039
(
1993
).
18.
S. V.
Levchenko
and
A. I.
Krylov
, “
Equation-of-motion spin-flip coupled-cluster model with single and double substitutions: Theory and application to cyclobutadiene
,”
J. Chem. Phys.
120
,
175
185
(
2004
).
19.
D.
Sinha
,
S.
Mukhopadhyay
, and
D.
Mukherjee
, “
A note on the direct calculation of excitation energies by quasi-degenerate mbpt and coupled-cluster theory
,”
Chem. Phys. Lett.
129
,
369
374
(
1986
).
20.
A. I.
Krylov
, “
Equation-of-motion coupled-cluster methods for open-shell and electronically excited species: The Hitchhiker’s guide to Fock space
,”
Annu. Rev. Phys. Chem.
59
,
433
462
(
2008
).
21.
H. J.
Werner
and
P. J.
Knowles
, “
An efficient internally contracted multiconfiguration–reference configuration interaction method
,”
J. Chem. Phys.
89
,
5803
5814
(
1988
).
22.
P. J.
Knowles
and
N. C.
Handy
, “
A determinant based full configuration interaction program
,”
Comput. Phys. Commun.
54
,
75
83
(
1989
).
23.
P. J.
Knowles
and
H.-J.
Werner
, “
Internally contracted multiconfiguration-reference configuration interaction calculations for excited states
,”
Theor. Chim. Acta
84
,
95
103
(
1992
).
24.
S.
Goedecker
, “
Linear scaling electronic structure methods
,”
Rev. Mod. Phys.
71
,
1085
1123
(
1999
).
25.
Linear-Scaling Techniques in Computational Chemistry and Physics
, edited by
R.
Zalesny
,
M. G.
Papadopoulos
,
P. G.
Mezey
, and
J.
Leszczynski
(
Springer Netherlands
,
2011
).
26.
P.
Pulay
, “
Convergence acceleration of iterative sequences. The case of SCF iteration
,”
Chem. Phys. Lett.
73
,
393
(
1980
).
27.
V. R.
Saunders
and
I. H.
Hillier
, “
A level-shifting method for converging closed shell Hartree–Fock wave functions
,”
Int. J. Quantum Chem.
7
,
699
705
(
1973
).
28.
G. B.
Bacskay
, “
A quadratically convergent Hartree-Fock (QC-SCF) method. Application to closed shell systems
,”
Chem. Phys.
61
,
385
404
(
1981
).
29.
J.
Šimunek
and
J.
Noga
, “
Hartree-Fock via variational coupled cluster theory: An alternative way to diagonalization free algorithm
,”
AIP Conf. Proc.
1504
,
143
151
(
2012
).
30.
A. D.
Rabuck
and
G. E.
Scuseria
, “
Improving self-consistent field convergence by varying occupation numbers
,”
J. Chem. Phys.
110
,
695
(
1999
).
31.
L.
Thøgersen
,
J.
Olsen
,
D.
Yeager
,
P.
Jorgensen
,
P.
Salek
, and
T.
Helgaker
, “
The trust-region self-consistent field method: Towards a black-box optimization in Hartree-Fock and Kohn-Sham theories
,”
J. Chem. Phys.
121
,
16
27
(
2004
).
32.
K. N.
Kudin
,
G. E.
Scuseria
, and
E.
Cancès
, “
A black-box self-consistent field convergence algorithm: One step closer
,”
J. Chem. Phys.
116
,
8255
(
2002
).
33.
J. D.
Whitfield
,
N.
Schuch
, and
F.
Verstraete
, “
The computational complexity of density functional theory
,” in
Many-Electron Approaches in Physics, Chemistry and Mathematics: A Multidisciplinary View
, Lecture Notes in Physics, edited by
V.
Bach
and
L. D.
Site
(
Springer
,
2014
), pp.
245
260
.
34.
N.
Schuch
and
F.
Verstraete
, “
Computational complexity of interacting electrons and fundamental limitations of density functional theory
,”
Nat. Phys.
5
,
732
(
2009
), also see the Appendix of arxiv:0712.0483.
35.
J. D.
Whitfield
and
Z.
Zimborás
,
J. Chem. Phys.
141
,
234103
(
2014
).
36.
J. D.
Whitfield
,
P. J.
Love
, and
A.
Aspuru-Guzik
, “
Computational complexity in electronic structure
,”
Phys. Chem. Chem. Phys.
15
,
397
411
(
2013
).
37.

In this article, we use the term problem to refer to a specific collection of problem instances. That is to say, the Hartree–Fock problem is the collection of all instances of Hartree–Fock, i.e., Eq. (1) with a fixed choice {hij, hijkl}.

38.
M.
Sipser
,
Introduction to the Theory of Computation
(
PWS Publishing Company
,
1997
).
39.
C. H.
Bennett
,
E.
Bernstein
,
G.
Brassard
, and
U.
Vazirani
, “
Strengths and weaknesses of quantum computing
,”
SIAM J. Comput.
26
,
1510
1524
(
1997
).
40.
Google AI Quantum and Collaborators, “
Hartree-fock on a superconducting qubit quantum computer
,”
Science
369
,
1084
1089
(
2020
).
42.
Q.
Sun
, “
Co-iterative augmented Hessian method for orbital optimization
,” arXiv:1610.08423 [physics.chem-ph] (
2016
).
43.
S.
Lloyd
, “
Universal quantum simulators
,”
Science
273
,
1073
1078
(
1996
).
44.
E.
Cancès
,
M.
Defranceschi
,
W.
Kutzelnigg
,
C. L.
Bris
, and
Y.
Maday
, “
Computational quantum chemistry: A primer
,” in
Special Volume, Computational Chemistry
, Handbook of Numerical Analysis Vol. 10 (
Elsevier
,
2003
), pp.
3
270
.
45.

Each Fock state corresponds to a Slater determinant ΨJ(x1,x2,xN)=Ω|ψ^(x1)ψ^(x2)ψ^(xN)|ΨJ/N! once orbitals {ϕj(x)} have been assigned to operators {aj} via field modes ψ^(x)=jϕj(x)aj.

46.
W. J.
Hehre
,
R. F.
Stewart
, and
J. A.
Pople
, “
Self-consistent molecular-orbital methods. I. Use of Gaussian expansions of slater-type atomic orbitals
,”
J. Chem. Phys.
51
,
2657
2664
(
1969
).
47.
S.
Veeraraghavan
and
D. A.
Mazziotti
, “
Global solutions of Hartree-Fock theory and their consequences for strongly correlated quantum systems
,”
Phys. Rev. A
89
,
010502
(
2014
).
48.

If S1, then a bond density matrix in this space P̃ satisfies Tr(P̃S)=N/2 and P̃=P̃SP̃. We can convert to an orthogonal basis using an orthogonalizing transform X, where XSX = 1, e.g., canonical orthogonalization Xcan=USs1/2 with S=USsUS. Then, we orthogonalize as P=XP̃X.

49.

To derive Eq. (9), one must substitute the orbital expansion (3) into Eq. (2) and then evaluate Eq. (1) using the commutator relations.

50.
A.
Szabo
and
N.
Ostlund
,
Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory
(
Dover Publications
,
1996
).
51.
T. S.
Hardikar
and
E.
Neuscamman
, “
A self consistent field formulation of excited state mean field theory
,”
J. Chem. Phys.
153
,
164108
(
2020
).
52.
J. A. R.
Shea
,
E.
Gwin
, and
E.
Neuscamman
, “
A generalized variational principle with applications to excited state mean field theory
,”
J. Chem. Theory Comput.
16
,
1526
1540
(
2020
).
53.
I. D.
Kivlichan
,
J.
McClean
,
N.
Wiebe
,
C.
Gidney
,
A.
Aspuru-Guzik
,
G. K.-L.
Chan
, and
R.
Babbush
, “
Quantum simulation of electronic structure with linear depth and connectivity
,”
Phys. Rev. Lett.
120
,
110501
(
2018
).
54.
R. A.
Horn
and
C. R.
Johnson
,
Matrix Analysis
(
Cambridge University Press
,
2005
).
55.
Z.
Jiang
,
K. J.
Sung
,
K.
Kechedzhi
,
V. N.
Smelyanskiy
, and
S.
Boixo
, “
Quantum algorithms to simulate many-body physics of correlated fermions
,”
Phys. Rev. Appl.
9
,
044036
(
2018
).
56.
Q.
Sun
,
X.
Zhang
,
S.
Banerjee
,
P.
Bao
,
M.
Barbry
,
N. S.
Blunt
,
N. A.
Bogdanov
,
G. H.
Booth
,
J.
Chen
,
Z.-H.
Cui
,
J. J.
Eriksen
,
Y.
Gao
,
S.
Guo
,
J.
Hermann
,
M. R.
Hermes
,
K.
Koh
,
P.
Koval
,
S.
Lehtola
,
Z.
Li
,
J.
Liu
,
N.
Mardirossian
,
J. D.
McClain
,
M.
Motta
,
B.
Mussard
,
H. Q.
Pham
,
A.
Pulkin
,
W.
Purwanto
,
P. J.
Robinson
,
E.
Ronca
,
E.
Sayfutyarova
,
M.
Scheurer
,
H. F.
Schurkus
,
J. E. T.
Smith
,
C.
Sun
,
S.-N.
Sun
,
S.
Upadhyay
,
L. K.
Wagner
,
X.
Wang
,
A.
White
,
J. D.
Whitfield
,
M. J.
Williamson
,
S.
Wouters
,
J.
Yang
,
J. M.
Yu
,
T.
Zhu
,
T. C.
Berkelbach
,
S.
Sharma
,
A.
Sokolov
, and
G. K.-L.
Chan
, “
Recent developments in the PYSCF program package
,”
J. Chem. Phys.
153
,
024109
(
2020
).
57.
Q.
Sun
,
T. C.
Berkelbach
,
N. S.
Blunt
,
G. H.
Booth
,
S.
Guo
,
Z.
Li
,
J.
Liu
,
J.
McClain
,
E. R.
Sayfutyarova
,
S.
Sharma
,
S.
Wouters
, and
G. K.-L.
Chan
, “
The python-based simulations of chemistry framework (PYSCF)
,” arXiv:1701.08223 [physics.chem-ph] (
2017
).
58.
J. R.
McClean
,
K. J.
Sung
,
I. D.
Kivlichan
,
Y.
Cao
,
C.
Dai
,
E. S.
Fried
,
C.
Gidney
,
B.
Gimby
,
T.
Häner
,
T.
Hardikar
,
V.
Havlíček
,
O.
Higgott
,
C.
Huang
,
J.
Izaac
,
Z.
Jiang
,
X.
Liu
,
S.
McArdle
,
M.
Neeley
,
T.
O’Brien
,
B.
O’Gorman
,
I.
Ozfidan
,
M. D.
Radin
,
J.
Romero
,
N.
Rubin
,
N. P. D.
Sawaya
,
K.
Setia
,
S.
Sim
,
D. S.
Steiger
,
M.
Steudtner
,
Q.
Sun
,
W.
Sun
,
D.
Wang
,
F.
Zhang
, and
R.
Babbush
, “
OpenFermion: The electronic structure package for quantum computers
,”
Quantum Science and Technology
5
,
034014
(
2020
).
59.
P.
Virtanen
,
R.
Gommers
,
T. E.
Oliphant
,
M.
Haberland
,
T.
Reddy
,
D.
Cournapeau
,
E.
Burovski
,
P.
Peterson
,
W.
Weckesser
,
J.
Bright
,
S. J.
van der Walt
,
M.
Brett
,
J.
Wilson
,
K. J.
Millman
,
N.
Mayorov
,
A. R. J.
Nelson
,
E.
Jones
,
R.
Kern
,
E.
Larson
,
C. J.
Carey
,
İ.
Polat
,
Y.
Feng
,
E. W.
Moore
,
J.
VanderPlas
,
D.
Laxalde
,
J.
Perktold
,
R.
Cimrman
,
I.
Henriksen
,
E. A.
Quintero
,
C. R.
Harris
,
A. M.
Archibald
,
A. H.
Ribeiro
,
F.
Pedregosa
,
P.
van Mulbregt
, and
SciPy 1.0 Contributors
, “
SciPy 1.0: Fundamental algorithms for scientific computing in Python
,”
Nat. Methods
17
,
261
272
(
2020
).
60.
D. W.
Scott
,
Multivariate Density Estimation: Theory, Practice, and Visualization
(
John Wiley & Sons
,
2015
).
61.
R. E.
Stanton
, “
Multiple solutions to the Hartree-Fock problem. I. General treatment of two-electron closed-shell systems
,”
J. Chem. Phys.
48
,
257
262
(
1968
).
62.
L. K.
Grover
, “
A fast quantum mechanical algorithm for database search
,” in
Proceedings of the 28th ACM Symposium on Theory of Computing (STOC’96)
(
ACM
,
1996
), p.
212
, also see arXiv:quant-ph/9605043.
63.

Note that only the global minima are constructed to have any connection to the energies of the Hamiltonian.

64.
J. R.
McClean
,
S.
Boixo
,
V. N.
Smelyanskiy
,
R.
Babbush
, and
H.
Neven
, “
Barren plateaus in quantum neural network training landscapes
,”
Nat. Commun.
9
,
4812
(
2018
).
65.
J.
Yang
,
J.
Brown
, and
J. D.
Whitfield
, “
Measurement on quantum devices with applications to time-dependent density functional theory
,” arXiv:1909.03078 (
2019
).
66.
J.
Brown
,
J.
Yang
, and
J. D.
Whitfield
, “
Solver for the electronic V-representation problem of time-dependent density functional theory
,”
J. Chem. Theory Comput.
16
,
6014
(
2020
).
67.
J. D.
Whitfield
,
M.-H.
Yung
,
D. G.
Tempel
,
S.
Boixo
, and
A.
Aspuru-Guzik
, “
Computational complexity of time-dependent density functional theory
,”
New J. Phys.
16
,
083035
(
2014
).