Global Solutions of Restricted Open-shell Hartree-Fock Theory from Semidefinite Programming with Applications to Strongly Correlated Quantum Systems

# Global solutions of restricted open-shell Hartree-Fock theory from semidefinite programming with applications to strongly correlated quantum systems

## Abstract

We present a density matrix approach for computing global solutions of restricted open-shell Hartree-Fock (ROHF) theory, based on semidefinite programming (SDP), that gives upper and lower bounds on the Hartree-Fock energy of quantum systems. While wave function approaches to Hartree-Fock theory yield an upper bound to the Hartree-Fock energy, we derive a semidefinite relaxation of Hartree-Fock theory that yields a rigorous lower bound on the Hartree-Fock energy. We also develop an upper-bound algorithm in which Hartree-Fock theory is cast as an SDP with a nonconvex constraint on the rank of the matrix variable. Equality of the upper- and lower-bound energies guarantees that the computed solution is the globally optimal solution of Hartree-Fock theory. The work extends a previously presented method for closed-shell systems [S. Veeraraghavan and D. A. Mazziotti, Phys. Rev. A 89, 010502(R) (2014))]. For strongly correlated systems the SDP approach provides an alternative to the locally optimized Hartree-Fock energies and densities with a certificate of global optimality. Applications are made to the potential energy curves of C, CN, Cr and NO.

## I Introduction

The most widely used approach to obtain the ground-state Hartree-Fock (HF) energy and density matrix has been to solve the Euler-Lagrange equations associated with the minimization of the Hartree-Fock energy. In 1951 Roothaan Roothaan (1951, 1960) and Hall Hall (1951) proposed the first self-consistent-field (SCF) method to solve the Hartree-Fock equations. However, the method was soon discovered to converge only for well-behaved cases. Since then numerous algorithms have been proposed to modify the SCF method to improve its convergence properties, including level-shifting Saunders and Hillier (1973); Carbo, Hernandez, and Sanz (1977), damping Hartree (1957); Zerner and Hehenberger (1979) and direct inversion of the iterative subspace (DIIS) Pulay (1980, 1982). DIIS is currently the most popular SCF algorithm because of its computational efficiency in most cases. Nevertheless, it is not globally convergent, and in many cases it is known to fail even with a good initial guess. A method is globally convergent when it converges to a local minimum from any initial guess. Level-shifting is globally convergent for a large enough shift parameter Cancès (2000), but its speed of convergence decreases as the shift parameter increases.

Since SCF methods do not ensure an energy decrease at each iteration, a potentially more natural approach to solving the Hartree-Fock problem is to minimize the Hartree-Fock energy directly as a function of the density matrix using gradient or Hessian-based methods. In 1956 McWeeny McWeeny (1956) proposed direct-minimization methods Fletcher and Reeves (1964); Hillier and Saunders (1970); Igawa and Fukutome (1975); Seeger and Pople (1976); Camp and King (1981); Stanton (1981); Bacskay (1981), but they have not yet found wide applicability either due to their slow convergence or prohibitive cost. Recently, a combination of the monotonic energy decrease property of direct-minimization methods and the speed of SCF methods was achieved in the relaxed-constraints algorithms Cances and Le Bris (2000); Kudin, Scuseria, and Cances (2002) and trust-region methods Thogersen et al. (2004); Francisco, Martinez, and Martinez (2004, 2006) which are both globally convergent and efficient. However, none of these methods can certify that the computed solution is the global minimum.

Convex minimization problems possess the attractive property that the existence of a local minimum implies that it is the global minimum. Semidefinite programs (SDP) are a class of convex optimization problems in which a linear function of a positive semidefinite matrix is optimized subject to linear constraints. For an -electron system the minimization of the ground-state energy as a functional of the two-electron reduced density matrix (2-RDM) subject to -representability constraints Mazziotti (2007); Coleman and Yukalov (2000); Mazziotti (2012a) has been expressed as an SDP Mazziotti (2007); Erdahl (1979); Mazziotti and Erdahl (2001); Nakata et al. (2001); Mazziotti (2002), and the resulting variational 2-RDM method Mazziotti (2007, 2012b); Garrod, Mihailović, and Rosina (1975); Erdahl (1979); Mazziotti and Erdahl (2001); Nakata et al. (2001); Mazziotti (2002); Zhao et al. (2004); Mazziotti (2004); ?; Cancés, Stoltz, and Lewin (2006); Shenvi and Izmaylov (2010); Verstichel et al. (2012); Mazziotti (2012a) in conjunction with large-scale SDP solvers Mazziotti (2004, 2011) has been applied to computing directly the 2-RDMs of strongly correlated systems including molecules like polyaromatic hydrocarbons Pelzer et al. (2011) and firefly luciferin Greenman and Mazziotti (2010) as well as quantum dots Rothman and Mazziotti (2008), quantum phase transitions Gidofalvi and Mazziotti (2006), and one- and two-dimensional spin models Hammond and Mazziotti (2006); Verstichel et al. (2012). We recently presented two SDP algorithms that yield upper and lower bounds on the ground-state energy from Hartree-Fock theory. Here we extend these algorithms to treat restricted open-shell Hartree-Fock (ROHF) theory. While wave function approaches to Hartree-Fock theory yield an upper bound to the Hartree-Fock energy, we derive a semidefinite relaxation of Hartree-Fock theory that yields a rigorous lower bound on the Hartree-Fock energy. In the lower-bound SDP algorithm the idempotency constraint on the one-electron density matrix is relaxed. We also develop an upper-bound algorithm in which Hartree-Fock theory is cast as an SDP with a nonconvex constraint on the rank of the matrix variable. Whenever the upper and lower bounds are equal to each other, they provide a certificate of global optimality to the obtained solution.

To illustrate these algorithms, we apply them to computing the symmetry-broken spin-restricted Hartree-Fock potential energy curves for C, CN, Cr and NO. This problem is challenging because there are multiple local minima of different spatial symmetries on the potential energy surfaces. Traditional methods that solve the Euler-Lagrange equations often converge to minima higher in energy than the global minimum. In principle, Čížek-Paldus stability analysis can be applied to locate multiple solutions to the Euler-Lagrange equations, but it can be computationally expensive and it cannot determine whether a local solution is also a global solution. We find that the upper-bound SDP algorithm consistently converges to the lowest energy solutions and that the lower-bound SDP algorithm generates a tight lower bound. Neither the upper- or lower-bound SDP algorithm relies on the quality of the initial guess for the density matrix, and in all SDP calculations presented here the initial density matrix is equated to a matrix whose elements are formed by a random number generator. When the upper and lower bounds are equal, the SDP algorithms provide a certificate of global optimality for the Hartree-Fock solution. The energetically degenerate symmetry-broken solutions are important because they can be combined convexly into an ensemble density matrix that not only has the desired molecular symmetry but yields a size-consistent energy. All of the Hartree-Fock solutions that we obtain using the SDP approach are restricted, i.e. has exactly the correct expectation value. In many strongly correlated cases, as shown in Section III, employing the symmetry-broken Hartree-Fock wave function as a reference in single-reference correlation methods like coupled cluster singles-doubles improves the correlated solution.

## Ii Theory

### ii.1 Canonical Hartree-Fock Theory

The quantity that we have been discussing as the density matrix can be more precisely called the one-particle reduced density matrix (1-RDM), which we denote as . The Hartree-Fock problem for an -electron system in an orthonormal basis of rank is typically expressed as the following minimization problem over the set of Hermitian matrices ()

 minimize1D∈Hr EHF(1D) (1) subjectto Tr(1D)=N (2) 1D2=1D (3)

The self-consistent-field Hartree-Fock method for an -electron system iteratively solves a system of Euler-Lagrange equations for a stationary point. The stationary point yields a ground-state Hartree-Fock energy and a set of occupied orbitals. The computed Hartree-Fock energy is not guaranteed to be the global-energy minimum. From the perspective of reduced density matrices (RDMs) Mazziotti (2007); Coleman and Yukalov (2000), we can understand the self-consistent-field method as iteratively checking extreme points of the set of 1-RDMs for satisfaction of the Euler-Lagrange equations where each extreme point corresponds to a 1-RDM with a Slater-determinant preimage Coleman (1963); Harriman (1978); Mazziotti (2012a). The set of extreme 1-RDMs (those with an -electron Slater determinant as a preimage) can be characterized by the idempotency constraint in Eq. (3).

### ii.2 SDP Hartree-Fock Theory

#### Convex relaxation

The optimization of the Hartree-Fock energy over the set of extreme 1-RDMs can be replaced without approximation by an optimization over the larger (and convex) set of -representable 1-RDMs (those with any -electron wave function as a preimage) Lieb (1981); Cancès (2000):

 minimize1D,1Q∈Hr+ EHF(1D) (4) subjectto Tr(1D)=N (5) 1D+1Q=I (6)

where is the following quadratic function of the 1-RDM:

 EHF(1D) =r∑ij1Kij1Dij+r∑ijkl1Dik2Vikjl1Djl (7) 1Kij =⟨i|^h|j⟩ (8) 2Vikjl =12(⟨ij|kl⟩−⟨ij|lk⟩). (9)

The one-electron Hamiltonian operator contains the kinetic energy operator and electron-nuclei potential, represents the electron-electron repulsion integrals, and the indices , , , and denote the orbitals in the one-electron basis set of rank . The notation , equivalent to and , indicates that both the 1-particle RDM and the 1-hole RDM are contained in the set of Hermitian positive semidefinite matrices.

The reduced-density-matrix formulation of Hartree-Fock theory can be recast as a convex semidefinite program by embedding the quadratic product of 1-RDMs in in a higher dimensional (two-electron) matrix . Rewriting as a linear functional of

 E(1D,2M) =Tr(1K1D)+Tr(2V2M), (10)

we can relax the non-convex Hartree-Fock optimization to a convex semidefinite program:

 minimize1D,1Q∈Hr+,2M∈Hr2+ E(1D,2M) (11) subjectto Tr(1D)=N (12) Tr(2M)≤N (13) 1D+1Q=I (14) r∑j=12Mikjj=N1Dik. (15)

The solution of this SDP relaxation yields a lower bound to the Hartree-Fock energy. Because the constraints on the matrix are minimal, this convex SDP formulation will typically yield energies that are significantly below the Hartree-Fock energy. To reproduce Hartree-Fock, further constraints on are required.

#### Upper-bound SDP algorithm

Two separate sets of additional conditions on the matrix that yield upper and lower bounds on the Hartree-Fock energy, respectively, will be considered. The first set of constraints, yielding the upper bound, consists of a single rank constraint

 rank(2M)=1. (16)

The matrix with its rank-one constraint and the contraction constraint in Eq. (15), we can show, is a tensor product of two identical 1-RDMs

 2Mikjl =1Dik1Djl. (17)

It follows that the solution of the optimization program in Eqs. (11-15) with the rank constraint in Eq. (16) is equivalent to the solution of the RDM formulation of Hartree-Fock theory in Eqs. (4-9). We have mapped Hartree-Fock theory exactly onto a rank-constrained semidefinite program (rc-SDP HF) Dattorro (2013). The rank-constrained semidefinite program is convex except for the rank constraint; the nonconvexity of the Hartree-Fock energy functional in the RDM formulation has been transferred to the rank restriction in the SDP formulation. Because of the rank constraint, the solution of rc-SDP HF is not necessarily a global solution, meaning that the solution can be a local minimum in the Hartree-Fock energy and hence, an upper bound on the global energy minimum. Unlike traditional formulations of Hartree-Fock theory, however, rc-SDP HF optimizes the 1-RDM over the convex set of -representable 1-RDMs, and in practice, we find that this difference makes it much more robust than traditional formulations in locating the global solution.

#### Lower-bound SDP algorithm

The second set of conditions, yielding a lower bound, consists of four constraints including

 r∑j=12Mijjk=1Dik (18)

and three additional constraints from permuting the indices and and/or and symmetrically. These convex conditions are a relaxation of the idempotency of the 1-RDM. They are necessary but not sufficient for the idempotency of the 1-RDM at the Hartree-Fock solution, and hence, optimization of the SDP program in Eqs. (11-15) with these additional constraints (lb-SDP) is an SDP relaxation of the reduced-density-matrix formulation of Hartree-Fock theory in Eqs. (4-9). The lb-SDP method yields a lower bound on the energy from the global Hartree-Fock solution. In practice, this lower bound is found to be quite tight, and in some cases it agrees exactly with the global Hartree-Fock solution. If lb-SDP produces a 1-RDM solution that is idempotent, then that solution is the global Hartree-Fock solution. Furthermore, when the upper and lower bounds from rc-SDP and lb-SDP agree, we have a guaranteed certificate that these computed bounds correspond to the global energy minimum of Hartree-Fock theory.

#### Closed- and open-shell spin restriction

We have formulated rc-SDP HF and lb-SDP in the spin-orbital basis set. To perform RHF and ROHF calculations in a spatial-orbital basis set, one needs to take into account the spin structure of the Hamiltonian and density matrices. For any RHF or ROHF calculation on an -electron system, , , and will have the following block structures:

 1D =[1Dα001Dβ] 2M =[2Mαα2Mαβ2Mtαβ2Mββ] (19) 1K =[1Kα001Kα] 2V =[2Vαα2Vαβ2Vαβ2Vαα]. (20)

For RHF the two blocks and the four blocks of are identical. Therefore the only necessary modifications in rc-SDP HF and lb-SDP are replacing by , by and rewriting as follows:

 E(1Dα,2Mαβ)=2Tr(1Kα1Dα)+Tr(2V2Mαβ) (21)
 2V =2(2Vαα+2Vαβ) (22) 2Vikjl =2⟨ij|kl⟩−⟨ij|lk⟩. (23)

For ROHF the and blocks of are not identical but (assuming ) because the spatial orbitals of paired electrons are required to be the same, the row space of is a subset of the row space of . In order to enforce this relation, is divided into closed-shell and open-shell blocks which are and respectively. With these blocks and SDP ROHF can be rewritten as follows:

 E(1Dc,1Do,2M)=2Tr(1Kα1Dc)+Tr(1Kα1Do)+Tr(2V2M) (24)
 minimize1Dc,1Do,1Q∈Hr+,2M∈Hr2+ E(1Dc,1Do,2M) subjectto Tr(1Dc)=Nβ Tr(1Do)=Nα−Nβ 1Dc+1Do+1Q=I r∑j=12Miαkαjαjα=Nα(1Dcik+1Doik) r∑j=12Miβkβjβjβ=Nβ1Dcik r∑j=12Miαkαjβjβ=Nβ(1Dcik+1Doik) r∑j=12Miβkβjαjα=Nα1Dcik .

For rc-SDP ROHF, the only additional constraint is on the rank, as in the spin-orbital formulation

 rank(2M)=1. (25)

For lb-SDP, the set of four constraints described earlier in Eq. (18) have to be enforced for all the four blocks of as follows:

 r∑j=12Miαjαjαkα =(1Dcik+1Doik) (26) r∑j=12Miβjβjβkβ =1Dcik (27) r∑j=12Miαjαjβkβ =1Dcik (28) r∑j=12Miβjβjαkα =1Dcik. (29)

Although formally there are 16 constraints in all, 6 of them are redundant due to the Hermiticity of resulting in 10 linearly independent constraints.

## Iii Applications

To illustrate the rc-SDP and lb-SDP methods, we apply them to computing the dissociation curves for C, CN, Cr and NO.

### iii.1 Methodology

The GAMESS electronic structure package is used to perform self-consistent-field Hartree-Fock calculations (SCF HF with DIIS) and coupled cluster singles doubles (CCSD) Purvis III and Bartlett (1982); Bartlett and Musiał (2007); Piecuch et al. (2002) calculations. The rc-SDP and lb-SDP are solved using the SDP solver RRSDP Mazziotti (2004). Since DIIS is the standard accelerator for SCF HF calculations, we compare rc-SDP HF results with DIIS results. Both rc-SDP HF and DIIS methods are performed without enforcing a specific spatial symmetry. The DIIS solution at the internuclear distance where is differentially larger than the distance is obtained by using the DIIS solution at as an initial guess.

The SDP solver RRSDP imposes the semidefinite constraint on each matrix through the factorization . For rc-SDP HF, the rank-one constraint on is readily enforced by defining to be a rectangular matrix. Scaling of RRSDP Mazziotti (2004) is determined by the matrix multiplication for the largest matrix block, which is for both rc-SDP and lb-SDP. For rc-SDP the rank of is one, and hence, the matrix multiplication scales approximately as . For lb-SDP the rank of scales as after applying the bound on the maximum rank from Pataki Pataki (1998) and Barvinok Barvinok (1995), and hence, the matrix multiplication scales approximately as .

### iii.2 C2 stretch

Because the C molecule has many low-lying excited states, it is a significantly multireferenced system even at equilibrium, which makes it a challenging system for both Hartree-Fock and correlation energy calculations Gidofalvi and Mazziotti (2005). Figure 1 shows various RHF energy curves and the lower bound in the 6-31G* basis set as a function of the C-C bond distance. The D and D curves Abrams and Sherrill (2004), generated by seeding the scan of the potential energy surface at two different values of were obtained using DIIS. The rc-SDP curve corresponds to D for  Å, C for  Å  Å, D for  Å  Å, and to C for  Å. The rc-SDP curve bifurcates from the D curve and joins the D curve without ever being non-differentiable. For  Å and  Å  the lb-SDP solution is lower than the rc-SDP solution by less than 0.005 a.u. thereby certifying it to be the global minimum within that threshold. For the intermediate region, rc-SDP likely continues to give the globally optimal curve although we do not have a formal mathematical guarantee.

### iii.3 Cr2 stretch

Cr is known to be an extremely challenging molecule to describe correctly by ab initio electronic structure theory Goodgame and Goddard (1985); Roos (2003); Celani et al. (2004); Muller (2009); Andersson et al. (1994); Bauschlicher Jr. and Partridge (1994); Roos and Andersson (1995); Zgid et al. (2009); Angeli et al. (2006). Figure 2 shows the HF energy in the valence triple-zeta (TZV) SchÃ¤fer, Huber, and Ahlrichs (1994) basis set as a function of the Cr-Cr distance. The large number of HF solutions that are energetically close to each other, shown in Fig. 3, provides a novel characterization of the substantial multireference correlation in Cr. The number of energetically close HF solutions is comparable in the STO-6G basis, indicating that this feature is not significantly dependent upon the basis set. The solution found by DIIS has D symmetry for all whereas the solution found by rc-SDP has D symmetry for  Å and C symmetry for  Å. Although the rc-SDP solution is symmetry-broken for  Å, it is globally optimal within the bound provided by lb-SDP. Further verification of the rc-SDP solutions being HF minima is provided by the fact that DIIS is able to obtain them when they are employed as initial guesses.

Changes in Hartree-Fock energies and densities can impact correlation energy calculations in two ways: (1) any change in the Hartree-Fock energy changes the correlation energy by its very definition and (2) any change in the Hartree-Fock density (or the Hilbert space spanned by the molecular orbitals) changes the reference wave function employed in many-electron correlation methods including coupled cluster Purvis III and Bartlett (1982); Bartlett and Musiał (2007); Piecuch et al. (2002) and parametric RDM methods Mazziotti (2008, 2010). In this and other examples considered, rc-SDP helps to identify symmetry-broken Hartree-Fock solutions that often generate improved CCSD solutions. While rc-SDP may identify a piecewise smooth potential energy surface, each piece can be analytically continued to generate a smooth Hartree-Fock surface from which a smooth CCSD surface can be computed. Cr is known to be extremely challenging for single-reference methods like coupled cluster theories Scuseria and Schaefer III (1990). Figure 4 explores the effect of using the global symmetry-broken C solution rather than the local D solution as the reference wave function in CCSD. The results in Fig. 4 show that much of the failure noted in the literature can be attributed to the D reference wave function rather than CCSD. While CCSD with the D reference diverges beyond 1.7 Å , CCSD with the C reference at least yields a physically realistic dissociation curve for the ground state.

### iii.4 CN stretch

The CN radical is of astrophysical interest due its presence in the interstellar medium Combi (1980). Although its low-lying excited states Knowles et al. (1988) enhance its utility as a optical probe  Mele and Okabe (1969); West and Berry (1974) for studying various properties, they also make it a multireferenced system. Figure 5 shows various HF energy curves and the lower bound in the cc-pVDZ basis Dunning, Jr. (1989) plotted as a function of the C-N internuclear distance . The C, C and C’ curves were obtained using DIIS. The rc-SDP curve corresponds to C for Å and to C for all other values of . The rc-SDP method manages to obtain the lowest curve among three different Hartree-Fock curves for all values of . As is evident from the figure, the ground-state Hartree-Fock curve is a piecewise defined function of the three Hartree-Fock curves with there being three points of non-smoothness where the curves intersect at = 1.2, 1.8 and 2.25 Å . If DIIS is given the Hückel guess, it converges to different curves in different regions which are not always the lowest solutions for those regions. If DIIS is used to generate the curve from left to right with the solution for a smaller bond distance being the guess for a larger bond distance, only the C curve is obtained. Consistent generation of the C curve might be the reason why the C and C’ curves have not been previously reported ThÃ¸gersen and Olsen (2004). The fact that the lb-SDP curve is never lower than rc-SDP by more than 0.006 a.u. for  Å provides a certificate of global optimality for those rc-SDP points within that threshold. After 1.6 Å rc-SDP likely continues to give the globally optimal curve although we do not have a formal mathematical guarantee. This is corroborated by the fact that the rc-SDP (and C) curve is size consistent, meaning that it is asymptotically equal to exactly the sum of Hartree-Fock energies of doublet N and singlet C. The energy of doublet N and singlet C is the same from both rc-SDP and DIIS and certified by lb-SDP to be globally optimal within 0.002 a.u. and 0.0004 a.u. respectively.

It is also worth noting that rc-SDP (and the C curve) does not dissociate CN into quadruplet N and triplet C in spite of them being lower in energy than doublet N and singlet C respectively. This is not due to convergence to a local minimum but instead is due to the inability of ‘restricted’ orbitals in ROHF (and RHF) to dissociate a molecule into fragments which have spatially separated electron pairs. This is a consequence of using the same spatial orbital to describe electron pairs. Since (doublet) CN has one unpaired electron, ROHF can describe its dissociation into fragments which have a total of one unpaired electron at most which is why it dissociates CN into doublet N and singlet C.

Figure 6 shows the CCSD curves obtained using the C and C Hartree-Fock solutions of which the latter was identified using rc-SDP. As is evident from the figure, although the CCSD with the C reference curve is lower in energy for R  1.4 Å , it rises rapidly after that point as it dissociates into charged species. The CCSD with the C reference curve is lower in energy after 1.4 Å and gives a qualitatively correct representation of dissociation into neutral species. It is also worth noting that the C and C Hartree-Fock curves cross at 1.2 Å, which is before the corresponding CCSD curves cross.

### iii.5 No2 bend

The NO radical is known to have a complicated, extensively studied photochemistry Wilkinson and J. Whitaker (2010). It has a conical intersection between the ground and first excited states Worner et al. (2011); Kraus et al. (2012). Figure 7 shows various HF energy curves (C symmetry) and the lower bound in the cc-pVDZ basis Dunning, Jr. (1989) plotted as a function of the O-N-O angle () for a symmetric configuration with a N-O bond length of 1.197 Å . Despite the existence of multiple HF minima which are energetically close, rc-SDP obtains the lowest minimum for all bond angles. For lb-SDP certifies the C curve from both DIIS and rc-SDP to be globally optimal. Furthermore, since lb-SDP is never lower than rc-SDP by more than 0.01 a.u., the entire rc-SDP curve is globally optimal within that threshold. Although the rc-SDP curve for corresponds to a saddle point on the complete NO potential energy surface (as a function of the two bond lengths in addition to the bond angle) it is indeed the global minimum (within 0.01 a.u.) for the fixed values of N-O bond lengths used in the calculation.

## Iv Discussion

An RDM formulation of Hartree-Fock theory, based on semidefinite programming, has been presented that yields upper and lower bounds on the Hartree-Fock solution. When these bounds are equal, they provide a certificate guaranteeing the globally optimal Hartree-Fock solution. As electrons become more strongly correlated, methods for Hartree-Fock based on the self-consistent-field approach like DIIS converge to stationary points of Hartree-Fock theory with potentially non-global energies and densities. In most instances we have been able to certify global optimality for known solutions of Hartree-Fock theory for the first time. Although there are methods to determine whether the obtained stationary point is a local minimum, maximum, or saddle point Čížek and Paldus (1967); Seeger and Pople (1977), our approach is unique in that it certifies global optimality.

Semidefinite relaxation of Hartree-Fock theory, which we derived in Section II.2.3, yields a rigorous lower bound on the Hartree-Fock energy. In contrast, wave function approaches to Hartree-Fock theory, such as the traditional optimization of a Slater determinant, yield upper bounds on the Hartree-Fock energy. Minimization of the electronic energy with respect to the orbitals of a Slater determinant generates a local stationary point which may or may not be the global minimum. Higher derivatives, such as those found in stability analysis, can be employed to search for additional local minima, each of which provides an upper bound on the energy of the global minimum. While not previously developed, the lower-bound approach enables us in many cases to certify that a solution is the global minimum of Hartree-Fock theory. If the 1-RDM obtained by the lower-bound SDP algorithm is idempotent, then the 1-RDM and its associated energy represent the global solution to the Hartree-Fock calculation. Furthermore, even if the 1-RDM is not idempotent, agreement of the lower-bound energy with the upper-bound energy from either a traditional wave function-based Hartree-Fock calculation or an SDP-based upper-bound calculation guarantees that the computed energy is the global minimum. As shown in the theory section, an upper bound to the Hartree-Fock energy can be computed through a rank-constrained SDP in which a nonconvex rank constraint is added to the optimization. Importantly, this upper-bound formulation shows that Hartree-Fock theory is convex except for the presence of the rank constraint.

Symmetry breaking and restoration can be employed to capture correlation effects at a lower computational cost Jiminez-Hoyos, Rodriguez-Guzman, and Scuseria (2013). Recently, they have been employed in variational quasi-particle theory Scuseria et al. (2011) to compute the ground-state energies from an antisymmetrized geminal power wave functionMazziotti (2000); ? at an computational cost where is the number of orbitals. In the presence of strong electron correlation the lowest energy solutions of Hartree-Fock theory can be spatially symmetry broken. We employ the SDP methods to distinguish the global solution from multiple local solutions of different spatial symmetries. The symmetry breaking generates multiple wave functions at the global minimum that are energetically degenerate. In the present case symmetry restoration can be accomplished by two methods. First, the full molecular symmetry can be reestablished by taking the ensemble of the energetically degenerate symmetry-broken solutions. The ensemble nature of the ground-state density matrix is a consequence of pursuing a mean-field description of the strongly correlated system. Second, a linear combination of the symmetry-broken solutions can be taken to generate a wave function by a non-orthogonal configuration interaction Broer and Nieuwpoort (1988); Hollauer and Nascimento (1993). In the second approach the degenerate Hartree-Fock solutions become entangled to form a pure density matrix composed of a single correlated wave function. While in the present work we pursue the first approach, the second approach provides insight into how the different symmetry-broken solutions of Hartree-Fock theory contribute information to the correlated ground-state wave function.

Direct computation of the two-electron reduced density matrix (2-RDM) has been previously accomplished by minimizing the energy as a functional of the 2-RDM subject to -representability conditions Mazziotti (2007, 2012b); Garrod, Mihailović, and Rosina (1975); Erdahl (1979); Mazziotti and Erdahl (2001); Nakata et al. (2001); Mazziotti (2002); Zhao et al. (2004); Mazziotti (2006); Cancés, Stoltz, and Lewin (2006); Shenvi and Izmaylov (2010); Verstichel et al. (2012); Mazziotti (2012a). Constrained optimization is performed by SDP. Although we have taken a different path in the derivation of the SDP algorithms for Hartree-Fock theory, they can be viewed within variational 2-RDM theory as the addition of further constraints on the 2-RDM to ensure that it represents the mean-field (or Hartree-Fock) limit. The upper-bound algorithm requires a nonconvex rank constraint, and the lower-bound algorithm requires relaxed idempotency conditions. As described above, the combination of the upper-bound and lower-bound SDP algorithms provides a mechanism for certifying the global minimum of Hartree-Fock theory. We have shown that global solutions are useful for seeding either wave function or reduced density matrix methods for describing strongly correlated quantum systems. The SDP-based restricted closed- and open-shell Hartree-Fock method, described here, is directly extendable to an unrestricted Hartree-Fock method, which will be presented elsewhere.

###### Acknowledgements.
DAM gratefully acknowledges the National Science Foundation, the Army Research Office, and Microsoft Corporation for their generous support.

### References

1. C. C. J. Roothaan, Rev. Mod. Phys. 23, 69 (1951).
2. C. C. J. Roothaan, Rev. Mod. Phys. 32, 179 (1960).
3. G. Hall, Proc. R. Soc. A 205, 541 (1951).
4. V. R. Saunders and I. H. Hillier, Int. J. Quantum Chem. 7, 699 (1973).
5. R. Carbo, J. Hernandez,  and F. Sanz, Chem. Phys. Lett. 47, 581 (1977).
6. D. R. Hartree, The calculation of atomic structures, Structure of matter series (J. Wiley, 1957).
7. M. C. Zerner and M. Hehenberger, Chem. Phys. Lett. 62, 550 (1979).
8. P. Pulay, Chem. Phys. Lett. 73, 393 (1980).
9. P. Pulay, J. Comp. Chem. 3, 556 (1982).
10. E. Cancès, Mathematical Models and Methods for Ab Initio Quantum Chemistry, 1st ed., edited by M. Defranceschi and C. L. Bris (Springer, New York, 2000) Chap. 2.
11. R. McWeeny, Proc. R. Soc. A 235, 496 (1956).
12. R. Fletcher and C. M. Reeves, The Computer Journal 7, 149 (1964).
13. I. Hillier and V. Saunders, Proc. R. Soc. A 320, 161 (1970).
14. A. Igawa and H. Fukutome, Prog. Theor. Phys. 54, 1266 (1975).
15. R. Seeger and J. A. Pople, J. Chem. Phys. 65, 265 (1976).
16. R. N. Camp and H. F. King, J. Chem. Phys. 75, 268 (1981).
17. R. E. Stanton, J. Chem. Phys. 75, 3426 (1981).
18. G. B. Bacskay, Chem. Phys. 61, 385 (1981).
19. E. Cances and C. Le Bris, Int. J. Quantum Chem. 79, 82 (2000).
20. K. Kudin, G. Scuseria,  and E. Cances, J. Chem. Phys. 116, 8255 (2002).
21. L. Thogersen, J. Olsen, D. Yeager, P. Jorgensen, P. Salek,  and T. Helgaker, J. Chem. Phys. 121, 16 (2004).
22. J. B. Francisco, J. M. Martinez,  and L. Martinez, J. Chem. Phys. 121, 10863 (2004).
23. J. Francisco, J. Martinez,  and L. Martinez, J. Math. Chem. 40, 349 (2006).
24. D. A. Mazziotti, ed., Reduced-Density-Matrix Mechanics: With Application to Many-electron Atoms and Molecules, Advances in Chemical Physics, Vol. 134 (Wiley, New York, 2007).
25. A. J. Coleman and V. I. Yukalov, Reduced Density Matrices: Coulson’s Challenge (Springer-Verlag, New York, 2000).
26. D. A. Mazziotti, Phys. Rev. Lett. 108, 263002 (2012a).
27. R. M. Erdahl, Reports Math. Phys. 15, 147 (1979).
28. D. A. Mazziotti and R. M. Erdahl, Phys. Rev. A 63, 042113 (2001).
29. M. Nakata, H. Nakatsuji, M. Ehara, M. Fukuda, K. Nakata,  and K. Fujisawa, J. Chem. Phys. 114, 8282 (2001).
30. D. A. Mazziotti, Phys. Rev. A 65, 062511 (2002).
31. D. A. Mazziotti, Chem. Rev. 12, 244 (2012b).
32. C. Garrod, V. Mihailović,  and M. Rosina, J. Math. Phys. 10, 1855 (1975).
33. Z. Zhao, B. J. Braams, H. Fukuda, M. L. Overton,  and J. K. Percus, J. Chem. Phys. 120, 2095 (2004).
34. D. A. Mazziotti, Phys. Rev. Lett. 93, 213001 (2004).
35. D. A. Mazziotti, Phys. Rev. A 74, 032501 (2006).
36. E. Cancés, G. Stoltz,  and M. Lewin, J. Chem. Phys. 125, 064101 (2006).
37. N. Shenvi and A. F. Izmaylov, Phys. Rev. Lett. 105, 213003 (2010).
38. B. Verstichel, H. van Aggelen, W. Poelmans,  and D. V. Neck, Phys. Rev. Lett. 108, 213001 (2012).
39. D. A. Mazziotti, Phys. Rev. Lett. 106, 083001 (2011).
40. K. Pelzer, L. Greenman, G. Gidofalvi,  and D. A. Mazziotti, J. Phys. Chem. A 115, 5632 (2011).
41. L. Greenman and D. A. Mazziotti, J. Chem. Phys. 133, 164110 (2010).
42. A. E. Rothman and D. A. Mazziotti, Phys. Rev. A 78, 032510 (2008).
43. G. Gidofalvi and D. A. Mazziotti, Phys. Rev. A 74, 012501 (2006).
44. J. R. Hammond and D. A. Mazziotti, Phys. Rev. A 73, 062505 (2006).
45. A. J. Coleman, Rev. Mod. Phys. 35, 668 (1963).
46. J. E. Harriman, Phys. Rev. A. 17, 1249 (1978).
47. E. H. Lieb, Phys. Rev. Lett. 46, 457 (1981).
48. J. Dattorro, Convex Optimization & Euclidean Distance Geometry (Meboo, Palo Alto, 2013) Chap. 4, pp. 308–333.
49. G. D. Purvis III and R. J. Bartlett, J. Chem. Phys. 76, 1910 (1982).
50. R. J. Bartlett and M. Musiał, Rev. Mod. Phys. 79, 291 (2007).
51. P. Piecuch, S. A. Kucharski, K. Kowalski,  and M. MusiaÅ, Comp. Phys. Comm. 149, 71 (2002).
52. G. Pataki, Math. Oper. Res. 23, 339 (1998).
53. A. Barvinok, Discrete Comput. Geom. 13, 189 (1995).
54. G. Gidofalvi and D. A. Mazziotti, J. Chem. Phys. 122, 194104 (2005).
55. M. L. Abrams and C. D. Sherrill, J. Chem. Phys. 121, 9211 (2004).
56. M. M. Goodgame and W. A. Goddard, Phys. Rev. Lett. 54, 661 (1985).
57. B. O. Roos, Collection of Czechoslovak Chemical Communications 68, 265 (2003).
58. P. Celani, H. Stoll, H.-J. Werner,  and P. Knowles, Mol. Phys. 102, 2369 (2004).
59. T. Muller, J. Phys. Chem. A 113, 12729 (2009).
60. K. Andersson, B. Roos, P.-A. Malmqvist,  and P.-O. Widmark, Chem. Phys. Lett. 230, 391 (1994).
61. C. W. Bauschlicher Jr. and H. Partridge, Chem. Phys. Lett. 231, 277 (1994).
62. B. O. Roos and K. Andersson, Chem. Phys. Lett. 245, 215 (1995).
63. D. Zgid, D. Ghosh, E. Neuscamman,  and G. K.-L. Chan, J. Chem. Phys. 130, 194107 (2009).
64. C. Angeli, B. Bories, A. Cavallini,  and R. Cimiraglia, J. Chem. Phys. 124, 054108 (2006).
65. A. SchÃ¤fer, C. Huber,  and R. Ahlrichs, J. Chem. Phys. 100 (1994).
66. D. A. Mazziotti, Phys. Rev. Lett. 101, 253002 (2008).
67. D. A. Mazziotti, Phys. Rev. A. 81, 062515 (2010).
68. G. E. Scuseria and H. F. Schaefer III, Chem. Phys. Lett. 174, 501 (1990).
69. M. R. Combi, Astrophys. J. 241, 830 (1980).
70. P. J. Knowles, H.-J. Werner, P. J. Hay,  and D. C. Cartwright, J. Chem. Phys. 89, 7334 (1988).
71. A. Mele and H. Okabe, J. Chem. Phys. 51, 4798 (1969).
72. G. A. West and M. J. Berry, J. Chem. Phys. 61, 4700 (1974).
73. T. H. Dunning, Jr., J. Chem. Phys. 90, 1007 (1989).
74. L. ThÃ¸gersen and J. Olsen, Chem. Phys. Lett. 393, 36 (2004).
75. I. Wilkinson and B. J. Whitaker, Annu. Rep. Prog. Chem., Sect. C: Phys. Chem. 106, 274 (2010).
76. H. J. Worner, J. B. Bertrand, B. Fabre, J. Higuet, H. Ruf, A. Dubrouil, S. Patchkovskii, M. Spanner, Y. Mairesse, V. Blanchet, E. MÃ©vel, E. Constant, P. B. Corkum,  and D. M. Villeneuve, Science 334, 208 (2011).
77. P. M. Kraus, Y. Arasaki, J. B. Bertrand, S. Patchkovskii, P. B. Corkum, D. M. Villeneuve, K. Takatsuka,  and H. J. Wörner, Phys. Rev. A 85, 043409 (2012).
78. J. Čížek and J. Paldus, J. Chem. Phys. 47, 3976 (1967).
79. R. Seeger and J. A. Pople, J. Chem. Phys. 66, 3045 (1977).
80. C. A. Jiminez-Hoyos, R. Rodriguez-Guzman,  and G. E. Scuseria, J. Chem. Phys. 139, 204102 (2013).
81. G. E. Scuseria, C. A. Jimenez-Hoyos, T. M. Henderson, K. Samanta,  and J. K. Ellis, J. Chem. Phys. 135, 124108 (2011).
82. D. A. Mazziotti, J. Chem. Phys. 112 (2000).
83. D. A. Mazziotti, Chem. Phys. Lett. 338, 323 (2001).
84. R. Broer and W. Nieuwpoort, Theor. chim. acta 73, 405 (1988).
85. E. Hollauer and M. A. C. Nascimento, J. Chem. Phys. 99 (1993).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters