Spatial and spin symmetry breaking in semidefinite-programming-based Hartree-Fock theory
aa r X i v : . [ phy s i c s . c h e m - ph ] A p r Spatial and spin symmetry breaking in semidefinite-programming-based Hartree-Focktheory
Daniel R. Nascimento and A. Eugene DePrince III Department of Chemistry and Biochemistry, Florida State University, Tallahassee, FL 32306-4390
The Hartree-Fock problem was recently recast as a semidefinite optimization over the space ofrank-constrained two-body reduced-density matrices (RDMs) [Phys. Rev. A , 010502(R) (2014)].This formulation of the problem transfers the non-convexity of the Hartree-Fock energy functionalto the rank constraint on the two-body RDM. We consider an equivalent optimization over thespace of positive semidefinite one-electron RDMs (1-RDMs) that retains the non-convexity of theHartree-Fock energy expression. The optimized 1-RDM satisfies ensemble N -representability condi-tions, and ensemble spin-state conditions may be imposed as well. The spin-state conditions placeadditional linear and nonlinear constraints on the 1-RDM. We apply this RDM-based approach toseveral molecular systems and explore its spatial (point group) and spin ( ˆ S and ˆ S ) symmetrybreaking properties. When imposing ˆ S and ˆ S symmetry but relaxing point group symmetry,the procedure often locates spatial-symmetry-broken solutions that are difficult to identify usingstandard algorithms. For example, the RDM-based approach yields a smooth, spatial-symmetry-broken potential energy curve for the well-known Be–H insertion pathway. We also demonstratenumerically that, upon relaxation of ˆ S and ˆ S symmetry constraints, the RDM-based approach isequivalent to real-valued generalized Hartree-Fock theory. I. INTRODUCTION
Hartree-Fock theory holds an important place in quan-tum chemistry. It seldom provides a quantitative de-scription of electronic structure, but it serves as a use-ful starting point for more sophisticated electronic struc-ture methods, such as coupled-cluster theory [1–4]. Thecanonical form of the Hartree-Fock problem is that ofRoothan [5] and Hall [6], which involves the repeated di-agonalization of the Fock matrix. When combined withconvergence acceleration procedures such as the directinversion of the iterative subspace (DIIS) [7, 8] and fasttwo-electron repulsion integral generation (using, for ex-ample, graphical processing units [9–12]), Hartree-Focktheory can be routinely applied to molecular systems con-taining thousands of atoms [11]. Nevertheless, for largeenough systems, the diagonalization of the Fock matrixcan eventually become problematic and complicates thedevelopment of linearly scaling algorithms.The direct optimization of the one-electron reduced-density matrix (1-RDM) is an attractive alternative tothe iterative solution of the Roothan-Hall equations forthe molecular orbital coefficient matrix; this idea hasbeen widely explored since the 1950s [13–17]. The mostdesirable feature of a density-matrix-based approach isthat it avoids the diagonalization of the Fock matrix,thereby facilitating the development of linearly scalingalgorithms [18–21]. An immediate drawback, however,is that the 1-RDM associated with the lowest possibleenergy does not correspond to any antisymmetrized N -electron wave function, let alone one comprised of a singleSlater determinant. To obtain a physically meaningful re-sult, one must explicitly consider the N -representabilityof the 1-RDM [22]. Typically, N -representability (andidempotency) in density-matrix-based Hartree-Fock isacheived through “purification” of an approximately N - representable 1-RDM [17]. Alternatively, recent work[23–25] has demonstrated the utility of semidefinite pro-gramming (SDP) techniques for this problem. Thepresent work explores the SDP-based strategy.Veeraraghavan and Mazziotti [23] recently recast theHartree-Fock problem as a constrained optimization overthe space of positive semidefinite two-body matrices. Theelectronic energy was expressed as a linear functional ofthe 1-RDM and a two-body matrix (denoted M in Ref.23) that is related to the 1-RDM through a contraction.By restricting the rank of M , one can obtain a rigor-ous upper bound to the globally optimal Hartree-Focksolution, and this “rank-constrained SDP”[23] optimiza-tion has the same formal O ( k ) cost as the Roothan-Hallform of the problem. Here, k represents the dimension ofthe one-electron basis set. Alternatively, a lower boundcan be obtained by imposing a relaxed set of idempo-tency constraints on the full-rank matrix. When separatesemidefinite optimizations under each set of constraintsyield the same 1-RDM, that solution is guaranteed to bethe globally optimal one. This guarantee is desirable; theprice paid, though, is that the solution of the lower-boundproblem requires at least O ( k ) floating-point operations.In this work, we consider an alternate representationof the rank-constrained SDP problem that eliminatesthe consideration of the two-body matrix, M . The re-sult is a similar SDP algorithm with a non-linear objec-tive function and non-linear constraints on the 1-RDM(in the case where the expectation value of ˆ S is con-strained). The algorithm retains its formal O ( k ) scaling,and, like other density-matrix-based Hartree-Fock imple-mentations, it avoids the repeated diagonalization of theFock matrix. We demonstrate that the SDP-based ap-proach can be applied to several flavors of Hartree-Focktheory, including restricted, unrestricted, and generalizedHartree-Fock (RHF, UHF, and GHF, respectively), de-pending on which spin symmetries are imposed on the1-RDM. We validate the implementation by exploringits spatial (point group) and spin ( ˆ S and ˆ S ) symmetrybreaking properties in several molecular systems. II. THEORYA. Density-matrix-based Hartree-Fock theory
The electronic energy for the ground state of a many-electron system is a function of the 1-RDM ( γ ) and thetwo-electron reduced-density matrix (2-RDM, Γ ) E = X pqrs Γ prqs ( pq | rs ) + X pq γ pq h pq . (1)Here, ( pq | rs ) represents a two-electron repulsion integralin Mulliken notation, h pq represents the sum of the one-electron kinetic energy and electron-nuclear potential en-ergy integrals, and the indices p , q , r , and s run over allspin orbitals. The 1-RDM and 2-RDM are defined as γ pq = h Ψ | ˆ a † p ˆ a q | Ψ i , (2)Γ prqs = 12 h Ψ | ˆ a † p ˆ a † r ˆ a s ˆ a q | Ψ i , (3)where ˆ a † and ˆ a represent the fermionic creation and an-nihilation operators of second quantization, respectively.Consider now the cumulant decomposition of the 2-RDM Γ = γ ∧ γ + ∆ , (4)where the cumulant matrix, ∆ , represents pure two-body correlations, and the symbol ∧ represents anantisymmetric tensor product (or Grassman product)[26]. By ignoring the cumulant matrix, we arrive at astatistically-independent description of electron motion;the resulting energy expression is equivalent to that fromHartree-Fock theory E = 12 X pqrs ( γ pq γ rs − γ ps γ rq )( pq | rs ) + X pq γ pq h pq . (5)To obtain the optimal Hartree-Fock 1-RDM, one caninvoke the variational principle and minimize the energygiven by Eq. (5) with respect to the elements of the 1-RDM. Note, however, that this optimization should becarried out under the constraint that the density matrixbe idempotent and have a trace equal to the number ofelectrons. The idempotency condition is a specific man-ifestation of a more general requirement that any physi-cally meaningful density matrix should correspond to anantisymmetrized N -electron wave function. Necessaryensemble N -representability conditions require that theeigenvalues of the 1-RDM [the natural orbital (NO) oc-cupation numbers] lie between zero and one and sum tothe total number of electrons [22]. The bounds on NOoccupation numbers can be enforced by requiring that the 1-RDM be positive semidefinite and related to theone-hole density matrix, ¯ γ , by γ pq + ¯ γ qp = δ pq . (6)This one-hole density matrix, which must also be positivesemidefinite, is defined in second-quantized notation as¯ γ pq = h Ψ | ˆ a p ˆ a † q | Ψ i . (7)The SDP procedure outlined in Sec. II B enforces theseensemble N -representability conditions, rather than theidempotency condition. It is important to note, however,that these conditions do not guarantee idempotency. For-tunately, as discussed in Sec. IV, extensive numericaltests indicate that the minimization of the electronic en-ergy given by Eq. (5) with respect to the elements of the1-RDM under ensemble N -representability constraints always yields an idempotent 1-RDM.For a non-relativistic Hamiltonian, the exact wavefunction should have a well-definied total spin ( S ) andprojection of spin ( M S ). Hence, one can impose ad-ditional conditions on the 1-RDM that fix the particlenumber and spin state for the system. Particle numberand M S can be fixed according toTr( γ αα ) = N α , (8)and Tr( γ ββ ) = N β , (9)where the subscripts α and β refer to electrons of α and β spin, and γ αα and γ ββ represent the spin-conservingblocks of the 1-RDM, the full structure of which is γ = (cid:18) γ αα γ αβ γ βα γ ββ (cid:19) . (10)Note that the γ αβ and γ βα blocks are zero if the 1-RDMcorresponds to a wave function that is an eigenfunctionof ˆ S . In this case, the total spin quantum number isrelated to an off-diagonal trace of the αβ – αβ spin-blockof the 2-RDM [27, 28]; for a statistically-independent pairdensity, we have S ( S + 1) = 12 ( N α + N β ) + 14 ( N α − N β ) − Tr( γ αα γ ββ ) . (11)An ensemble N -representable 1-RDM corresponding toa state that is an eigenfunction of ˆ S should satisfy Eqs.(6) and (8)-(9); for the 1-RDM to represent an ensemblespin state with total spin, S , Eq. (11) should also besatisfied.Various flavors of Hartree-Fock can be classified ac-cording to the symmetries that are preserved by the wavefunction or density matrix [29–31]. If we relax the spinconstraint given by Eq. (11), we arrive at density-matrix-based UHF. If we also lift the constraint that the Hartree-Fock wave function be an eigenfunction of ˆ S , we obtainGHF [31–34]. In GHF, the particle number constraintsof Eqs. (8) and (9) reduce to a single constraint fixingthe total particle numberTr( γ αα ) + Tr( γ ββ ) = N (12)and γ and ¯ γ are no longer guaranteed to have a block-diagonal spin structure. While S may no longer be a goodquantum number, we can still evaluate the expectationvalue of ˆ S according to h ˆ S i = 34 N + 14 [Tr( γ αα ) − Tr( γ ββ )] −
14 Tr( γ αα γ αα ) −
14 Tr( γ ββ γ ββ )+ 12 Tr( γ αβ )Tr( γ βα ) + Tr( γ αβ γ βα ) − Tr( γ αα γ ββ ) . (13) B. Semidefinite optimization
The minimization of the electronic energy given by Eq.(5) subject to the constraints outlined above constitutes anonlinear semidefinite optimization. We adopt a matrix-factorization-based approach to this problem based uponthe “RRSDP” algorithm described in Refs. 35 and 36.The 1-RDM and one-hole density matrix are expressedas contractions of auxiliary matrices, d and ¯ d , as γ pq = X Q d Qp d Qq , (14)and ¯ γ pq = X Q ¯ d Qp ¯ d Qq , (15)and are thus positive semidefinite by construction. Theauxiliary matrices serve as the actual variable quantitiesin the optimization. Note that the spin structures of γ and ¯ γ depend on what spin symmetry is imposed. Whenˆ S symmetry is enforced, the density matrices consist oftwo spin-conserving blocks which can be separately fac-torized as in Eqs. (14) and (15); in this case, d and ¯ d arealso block diagonal, each comprised of two k × k blocks,where k represents the number of spatial basis functions.When ˆ S symmetry is broken, this block structure is lost,and, like the density matrices, the auxiliary matrices arecomprised of a single 2 k × k block.To obtain the optimal d and ¯ d , we minimize the aug-mented Lagrangian function L ( d , ¯ d ) = E ( d ) + X i (cid:20) µ C i ( d , ¯ d ) − λ i C i ( d , ¯ d ) (cid:21) , (16)with respect to variations in their respective elements.Here, the sum runs over all constraints, i , the symbol, C i ( d , ¯ d ), represents the error in constraint i , the symbol, λ i , represents the corresponding Lagrange multiplier, and µ is a penalty parameter. The optimization proceedsaccording to a two-step scheme that is similar to thatemployed in Refs. 35 and 36: 1. For a set of Lagrange multipliers { λ ( n ) i } and penaltyparameter µ ( n ) , minimize Eq. (16) with respect tothe elements of d and ¯ d .2. Update the Lagrange multipliers λ ( n +1) i = λ ( n ) i − C ( n ) i ( d , ¯ d ) /µ ( n ) (17)and the penalty parameter µ ( n +1) = f µ ( n ) , (18)where f is defined as f = . , if max {| C ( n ) i ( d ) |} max {| C ( n − i ( d ) |} < . g, otherwise. (19)The parameter, g , is a random number that lies onthe interval [0.08:0.12].Steps 1 and 2 are repeated until the error in the con-straints ( || C ( d , ¯ d ) || ) falls below 10 − and the energychanges between iterations by less than 10 − E h . In atypical optimization, there are fewer than 20 of thesemacroiterations. The parameter, g , has been introducedinto the RRSDP algorithm because of the non-convexnature of the Hartree-Fock problem. In cases where thealgorithm identifies a solution that is not the global so-lution, g introduces some non-deterministic behavior tofacilitate the identifaction of additional solutions in sub-sequent computations.In the present implementation, the minimization inStep 1 is achieved using the Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method, as imple-mented in the library liblbfgs . The L-BFGS routinerequires the repeated evaluation of the electronic energy,the constraints, and the gradient of the energy and con-straints with respect to the elements of d and ¯ d . Evaluat-ing the constraints and the derivative of the constraintswith respect to the elements of d and ¯ d requires only O ( k ) floating-point operations. The gradient of the en-ergy with respect to the elements of d is given by thematrix product of d and the Fock matrix: ∂E∂d Qp = 2 X q d Qq F qp (20)Hence, the rate limiting step in the algorithm is theconstruction of the Fock matrix, which requires O ( k )floating-point operations. Note that, as with the densitymatrices, the spin structure of the Fock matrix dependson whether or not ˆ S symmetry is enforced.The optimization is performed in the orthonormal ba-sis defined by L¨owdin’s symmetric orthogonalization, andthe initial d and ¯ d matrices are seeded with random num-bers on the interval [-1:1]. At the beginning of an opti-mization, the L-BFGS step can require hundreds or thou-sands of Fock matrix builds, but this number significantlydecreases near convergence of the macroiterations. Obvi-ously, this semidefinite procedure requires that the Fockmatrix be built far more times than would be requiredby a conventional Hartree-Fock algorithm. However, itdoes have the nice property that the diagonalization ofthe Fock matrix is avoided completely, unless, of course,the orbital energies themselves are desired at the end ofthe optimization. Should one wish to devise a linear-scaling algorithm based upon the formalism presentedherein, it would be desirable to avoid the diagonalizationof the overlap matrix as well; in this case, the ensemble N -representability conditions for the 1-RDM should begeneralized for non-orthogonal orbitals [25]. III. RESULTSA. Broken spatial symmetry
Artifactual symmetry breaking problems often arisein Hartree-Fock-based descriptions of strongly-correlatedsystems in which more than one electronic configurationis important in the full configuration interaction (CI)wave function. Small et al. [37] recently demonstratedthat complex orbitals resolve this issue for some well-known RHF-based cases. An illustrative example con-sidered in Ref. 37 is the potential energy curve (PEC)for the C insertion of Be into H , which has also servedas a popular model for multireference correlation studies[37–47]. Figure 1(a) illustrates the PEC as a function ofthe Be-H distance, x , with the geometry (in units of ˚A)defined as Be: (0.0, 0.0, 0.0) , H: ( x , 1.344 - 0.46 x , 0.0) , H: ( x , -1.344 + 0.46 x , 0.0) . (21)The curve labeled SDP-RHF was generated using thepresent density-matrix-based algorithm while imposingˆ S [Eq. (11)] and ˆ S symmetry [Eqs. (8)-(9)]. Here, theone-electron basis set was cc-pVDZ, and we employedthe density fitting approximation to the two-electron re-pulsion integrals (using the def2-TZVP-JK auxiliary ba-sis set). Full CI computations employed conventional 4-index integrals.As discussed in Ref. 37, the dominant configuration inthe full CI ground state is (1a ) (2a ) (1b ) at x = 0 . ) (2a ) (3a ) at x = 2 . ). Small et al. demon-strated that complex orbitals yield a PEC that smoothlyinterpolates between the two RHF limits, avoiding thecusp where these two configurations become degenerate.We see here that the present algorithm, which employsreal orbitals and density matrices, also yields a smoothPEC, and inspection of the 1-RDM reveals that thissmoothness is achieved through a break in spatial symme-try. Consider the point in the multireference region of thePEC at x = 1 . γ αα + γ ββ ) from the present procedure in thebasis of symmetry-pure orbitals obtained from conven-tional RHF with the configuration, (1a ) (2a ) (3a ) ,the first 4 × φ φ φ φ φ .
00 0 .
00 0 .
00 0 . φ .
00 1 . − .
13 0 . φ . − .
13 0 .
68 0 . φ .
00 0 .
08 0 .
92 1 . and b symmetry, which suggests that the symmetry of the over-all wave function has been reduced to at most C s . Evenwhen fully neglecting spatial symmetry in a conventionalRHF computation (using, for example, the implemen-tation in Psi4 [48]), one in general locates either onespatially-pure state or another, rather than the spatiallycontaminated global solution.Figure 1(b) illustrates the same curves generated usingsecond-order perturbation theory (MP2). The MP2 com-putations built upon orbitals from density-matrix-basedHartree-Fock (which are the eigenfunctions of γ ) yield aPEC that interpolates between distinct RHF-based MP2curves with (1a ) (2a ) (1b ) and (1a ) (2a ) (3a ) reference functions. Here, however, the curve in the mul-tireference region is not smooth; there is a distinct kinknear x = 1 . . The RHF curves correspond to configura-tions that differ only in the occupation of the π (1b symmetry) and π ∗ (1b symmetry) orbitals. As above,the PEC labeled SDP-RHF was generated using thepresent approach while enforcing both ˆ S and ˆ S symme-try, and we observe the same spatial symmetry breakingbehavior reported elsewhere. At equilibrium and dissoci-ation, the density-matrix-based approach yields solutionswith energies that agree with those of one of the spatiallypure RHF configurations or the other. However, for N-Nbond lengths in the range of 1.5–2.5 ˚A, we obtain solu-tions that are lower in energy than both configurations;this stabilization is achieved via a break in the symmetryof the π / π ∗ orbitals. FIG. 1: Potential energy curves for the C insertion of Be into H in the cc-pVDZ basis set. The curves are shifted such thatthe energy of the colinear geometry is 0 kcal mol − . Results are provided at the (a) RHF [using conventional and density-matrix-based (SDP-RHF) approaches] and (b) MP2 levels of theory. MP2 results were obtained using orbitals from bothapproaches (indicated in square brackets).
0 20 40 60 80 100 120 140 160 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 R e l a t i v e E ne r g y [ kc a l m o l − ] Be−H distance [Å] (a) RHF, (1a ) (2a ) (1b ) RHF, (1a ) (2a ) (3a ) SDP−RHFfull CI 0 20 40 60 80 100 120 140 160 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 R e l a t i v e E ne r g y [ kc a l m o l − ] Be−H distance [Å] (b) MP2 [ (1a ) (2a ) (1b ) ]MP2 [ (1a ) (2a ) (3a ) ]MP2 [ SDP−RHF ]full CI FIG. 2: Potential energy curves for the dissociation of N inthe cc-pVDZ basis set. a −109.00−108.80−108.60−108.40−108.20−108.00−107.80−107.60−107.40 1.0 1.5 2.0 2.5 3.0 3.5 4.0 E ne r g y [ E h ] N−N distance [Å](1a g ) (1b ) (2a g ) (2b ) (3a g ) (1b ) (1b ) (1a g ) (1b ) (2a g ) (2b ) (1b ) (3a g ) (1b ) SDP−RHF a RHF and SDP-RHF computations employ the density fittingapproximation and the cc-pVDZ-JK auxiliary basis set.
B. Broken spin symmetry
While it is desirable for the Hartree-Fock wave functionto retain all of the symmetries of the exact wave function,it is often useful, as in the case of describing moleculardissociation, to lift constraints on its spin symmetry ( ˆ S ).For example, it is well known that RHF in general doesnot yield the correct dissociation limit for closed shellmolecules; UHF usually delivers size consistent resultsin such cases, at the expense of retaining S as a goodquantum number. Interestingly, even UHF does not yielda size-consistent dissociation limit for some closed-shellmolecules, such as CO [51]. If a truly size-consistentsingle-determinant method is desired, one must be will-ing to break additional symmetries in the wave function,such as ˆ S . In this Section, we explore the spin-symmetry breaking properties of density-matrix-driven UHF andGHF (SDP-UHF and SDP-GHF, respectively) for onewell-studied [51] case: the dissociation of molecular oxy-gen. The 1-RDMs generated from SDP-UHF satisfyensemble N -representability conditions and the ˆ S con-straints of Eqs. (8)-(9), while those from SDP-GHF sat-isfy only ensemble N -representability conditions. For theremainder of this Section, we use the terms UHF/SDP-UHF and GHF/SDP-GHF interchangeably.Following Ref. 51, the inability of UHF to yield a size-consistent dissociation curve for O is easily understoodfrom simple spin arguments. At equilibrium, the groundstate of molecular oxygen is the triplet state ( Σ − g ), andthe ground-state of the dissociation limit involves twooxygen atoms in their triplet states ( P). The UHF wavefunction for a given multiplicity is taken to be the high-spin determinant, so O at equilibrium and each dissoci-ated oxygen atom all have M S = 1. At dissociation, thetwo M S = 1 fragments can only couple to yield statesof M S = 2 (the quintet state) and M S = 0 (the sin-glet state). Hence, the UHF triplet state cannot connectthe ground-state at equilibrium to those (the singlet orquintet) at dissociation.Figure 3 provides PECs for the dissociation of molec-ular oxygen corresponding to the lowest-energy UHFsinglet, triplet, and quintet states, as well as that forGHF. All computations were performed within the cc-pVDZ basis set using the present density-matrix-basedapproach to Hartree-Fock. As expected, the lowest-energy UHF curve at equilibrium is that for the tripletstate; the singlet UHF curve lies about 30 mE h higherin energy. At dissociation, the singlet and quintet curvesbecome degenerate and are lowest in energy, while thedissociation limit for the triplet lies about 40 mE h higherin energy. Hence, as simple spin arguments suggest, theUHF triplet PEC does not connect the lowest energy so-lutions at equilibrium and dissociation. The present re- FIG. 3: SDP-UHF and SDP-GHF potential energy curves forthe dissociation of O in the cc-pVDZ basis set. a −149.65−149.60−149.55−149.50−149.45−149.400.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 E ne r g y [ E h ] O−O distance [Å] SDP−UHF S=0SDP−UHF S=1SDP−UHF S=2SDP−GHF a SDP-UHF and SDP-GHF computations employ the density fittingapproximation and the cc-pVDZ-JK auxiliary basis set. sults are in excellent agreement with those of Ref. 51 andprovide strong evidence for the numerical equivalence be-tween density-matrix-driven and conventional GHF.
FIG. 4: The region of the potential energy curve for the dis-sociation of O where SDP-GHF provides a lower-energy so-lution than SDP-UHF (top panel). The value of h ˆ S i forSDP-GHF follows those for the SDP-UHF triplet and singletstates at bond lengths less than 1.44 ˚A and greater than 1.47˚A, respectively (bottom panel). −149.568−149.567−149.566−149.565−149.564−149.563−149.562−149.561−149.5601.44 1.45 1.46 1.47 1.48 E ne r g y [ E h ] SDP−UHF S=0 SDP−UHF S=1SDP−GHF1.11.21.31.41.51.61.71.81.92.02.11.44 1.45 1.46 1.47 1.48 Æ S æ O−O distance [Å]
The UHF triplet and singlet states cross at an O–Odistance slightly less than 1.46 ˚A. The top panel of Fig.4 provides a more detailed illustration of the relevantPECs in the vicinity of this state crossing. We see that,in the region where the ground state predicted by UHFchanges from the triplet state to the singlet state, GHFoffers a solution that is lower in energy than either UHFsolution. As with the PECs in Fig. 3, this result is inexcellent agreement with that from Ref. 51. The lowerpanel of Fig. 4 illustrates how h ˆ S i changes as the GHFPEC interpolates between the UHF triplet and singletstate curves. We note that while the qualitative behavior of h ˆ S i in this region resembles that presented in Ref. 51,there are quantitative differences between our computedvalues of h ˆ S i and theirs. For example, our computationsindicate that the expecation value of ˆ S merges with thatof the UHF singlet curve at an O–O distance less than1.47 ˚A, while the GHF and UHF values of h ˆ S i in Ref.[51] merge closer to 1.48 ˚A. Nonetheless, the results arequite similar. IV. CONCLUSIONS
We have considered the properties of a nonlinear repre-sentation of the rank-constrained semidefinite program-ming algorithm [23] for the Hartree-Fock problem. Therate limiting step of the SDP-based approach is the for-mation of the Fock matrix, so it shares the same for-mal O ( k ) scaling of the usual Roothan-Hall formula-tion of the problem. We demonstrated that several fla-vors of Hartree-Fock theory (e.g. RHF, UHF, GHF) canbe implemented within this density-matrix-based formal-ism depending on which spin symmetry conditions wechoose to enforce, and we validated our implementationby studying several spatial and spin symmetry breakingproblems.Interestingly, each flavor of Hartree-Fock theory is re-covered without explicitly enforcing any constraints onthe idempotency of the 1-RDM. We only enforce ensem-ble N -representability conditions, which apply to bothcorrelated and uncorrelated 1-RDMs alike. This re-sult contrasts with conventional approaches to density-matrix-based Hartree-Fock theory which rely on purifi-cation strategies to obtain idempotent density matrices.We note that we have verified numerically that the op-timized 1-RDMs were indeed idempotent in every com-putation performed in this work. We can rationalize theidempotency of our 1-RDMs by considering the struc-ture of the corresponding statistically independent pairdensity ( Γ = γ ∧ γ ). It can be show that the trace ofsuch a 2-RDM is minimized in the case of an idempo-tent 1-RDM and, accordingly, an energy minimizationprocedure will favor idempotent 1-RDMs because theyminimize electron-electron repulsions.Alternatively, the idempotency of the 1-RDMs canbe understood from the perspective of pure-state N -representability in RDM functional theory [52]. Inthis context, Valone noted [53] that the distinctionbetween pure-state and ensemble N -representability isunnecessary in the case of the exact RDM functional.More broadly, pure-state N -representability can beachieved under ensemble N -representability conditions,provided that the RDM functional is an appropriate one(i.e. that the functional is pure-state N -representable).Indeed, our observations are consistent with this propo-sition. The Hartree-Fock energy functional is the RDMfunctional that arises for a 1-RDM derived from a singleSlater determinant; the Hartree-Fock energy functionalis thus pure-state N -representable. The pure-state N -representability of the functional leads to additionaldesirable properties. For example, the trace constraintfor GHF that defines the total particle number [Eq.(12)] technically only enforces the expectation value of N . To specify the particle number exactly, the variancein N , h ˆ N i − h ˆ N i , should vanish. In the case that(i) the 2-RDM is expressed as an antisymmetrizedproduct of the 1-RDM with itself and (ii) the 1-RDM isidempotent, this variance is exactly zero. As discussedin the Appendix, similar arguments can be made forthe exact specification of the M S and S ( S + 1). Hence,with the choice of the Hartree-Fock energy functionalas the RDM functional, the application of ensemble N -representability and ensemble spin constraints yieldspure-state N - and S -representable 1-RDMs. Appendix A: On the exact specification of N , M S ,and S ( S + 1) In the semidefinite-programming-based approach toHartree-Fock theory, the particle number and spin statefor the system are specified through constraints on theexpectation values of ˆ N , ˆ S , and ˆ S . The exact speci-fication of these quantities technically requires that thecorresponding variances (e.g. h ˆ N i − h ˆ N i ) be zero. Inthe case that the 1-RDM is idempotent and 2-RDM canbe constructed as an antisymmetrized product of the 1-RDM with itself, we can easily show that the variancesfor ˆ N and ˆ S exactly vanish. The expectation value ofthe number operator squared leads to h ˆ N i = h ˆ N i + h ˆ N i − Tr( γγ ) , (A1)where we have used the fact that, at the Hartree-Focklevel of theory, Γ = γ ∧ γ , and we have assumed thattrace relations such as Eqs. (8) and (9) or Eq. (12)are satisfied. If the 1-RDM is also idempotent, the firstand last terms in Eq. (A1) cancel, and the variance iszero. The same analysis holds for number operators cor-responding electrons of α and β spin individually ( ˆ N α and ˆ N β , respectively).The variance in M S is defined by h ˆ S i − h ˆ S i , whereˆ S = ( ˆ N α − ˆ N β ), and M S = h ˆ S i , when Eqs. (8) and(9) are satisfied. If the 2-RDM is expressible in terms ofthe 1-RDM, we have h ˆ S i = 14 (cid:18) h ˆ N α i − Tr( γ α γ α ) + h ˆ N β i − Tr( γ β γ β )+ h ˆ N α i + h ˆ N β i − h ˆ N α ih ˆ N β i (cid:19) . (A2)In the case that 1-RDM is idempotent, Eq. (A2) reducesto h ˆ S i = 14 (cid:18) h ˆ N α i − h ˆ N β i (cid:19) = h ˆ S i , (A3) FIG. 5: The (a) total energy (E h ), (b) expectation value ofˆ S , and (c) the Euclidean norm of the error in the pure spinconditions of Eq. (A6) for the dissociation of molecular ni-trogen in the cc-pVDZ basis set. a -109.0-108.8-108.6-108.4-108.2-108.0-107.8-107.6 1.0 1.5 2.0 2.5 3.0 E n e r g y [ E h ] N-N distance [˚A](a)0.00.51.01.52.02.53.0 1.0 1.5 2.0 2.5 3.0 h ˆ S i N-N distance [˚A](b)0.00.20.40.60.81.01.21.41.61.8 1.0 1.5 2.0 2.5 3.0 || h ˆ a † k ˆ a l ˆ S + i || N-N distance [˚A](c) SDP-RHFSDP-UHFSDP-RHFSDP-UHFSDP-RHFSDP-UHF a SDP-RHF and SDP-UHF computations employ the density fittingapproximation and the cc-pVDZ-JK auxiliary basis set. and the variance in M S vanishes.A similar analysis of the variance in S ( S + 1) is moreinvolved, as the expectation value of ˆ S depends upon thefour-particle RDM. Instead, we consider the two-particlepure-spin state conditions detailed in Ref. 54, the firstof which is a “contraction condition” requiring that anypure spin state, | Ψ SM S i , should satisfy h Ψ SM S | ˆ a † k ˆ a l ( N ˆ S − M S ˆ N ) | Ψ SM S i = 0 , (A4) FIG. 6: The (a) total energy (E h ), (b) expectation value ofˆ S , and (c) the Euclidean norm of the error in the pure spinconditions of Eq. (A6) for the dissociation of hydroxyl radicalin the cc-pVDZ basis set. a -75.40-75.35-75.30-75.25-75.20-75.15-75.10-75.05-75.000.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 E n e r g y [ E h ] O-H distance [˚A](a)0.801.001.201.401.601.800.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 h ˆ S i O-H distance [˚A](b)0.000.200.400.600.801.000.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 || h ˆ a † k ˆ a l ˆ S + i || O-H distance [˚A](c) SDP-ROHFSDP-UHFSDP-ROHFSDP-UHFSDP-ROHFSDP-UHF a SDP-ROHF and SDP-UHF computations employ the density fittingapproximation and the cc-pVDZ-JK auxiliary basis set. for all spin orbitals, k and l . As with the variances above,this condition is automatically satisfied in the case that the 2-RDM is represented as an antisymmetrized productof an idempotent 1-RDM with itself and that Eqs. (8)and (9) are satisfied. A less trival constraint requires thatthe maximal spin projection for a spin pure state satisfies h Ψ SM S | ˆ a † k ˆ a l ˆ S + | Ψ SM S i = 0 , (A5)for all spin orbitals, k and l . For a 1-RDM with ˆ S sym-metry and a 2-RDM expressible in terms of the 1-RDM,this set of constraints reduces to ∀ k, l : γ k β l β − X p γ p α l α γ k β p β = 0 , (A6)where the Greek subscripts denote the spin componentof the orbital label. These constraints may be satisfiedfor a singlet state, where γ αα = γ ββ , and the 1-RDMis idempotent. However, these constraints are not ob-viously satisfied for other spin states. Here, we demon-strate numerically that these constraints are satisfied inSDP-based Hartree-Fock through constraints on the ex-pectation values of ˆ S and ˆ S .Figure 5(a) illustrates the potential energy curve forthe dissociation of molecular nitrogen at the SDP-RHFand SDP-UHF levels of theory. The expectation valuesof ˆ S from RHF and UHF [Fig. 5(b)] diverge as therespective energies diverge, beyond the Coulson-Fischerpoint. Figure 5(c) shows the norm of the constraintsgiven by Eq. (A6); we can see that these constraintsare satisfied by SDP-RHF (the norm is zero), whereasthe norm computed from the SDP-UHF 1-RDM becomesquite large beyond the Coulson-Fischer point.Lastly, we consider a similar analysis for a non-singletcase: the dissociation of the hydroxyl radical. Figure6 provides results obtained using SDP-based restrictedopen-shell Hartree-Fock (SDP-ROHF) and SDP-UHF.The 1-RDM obtained from SDP-ROHF satisfies Eqs.(8) and (9) as well as the ˆ S constraint of Eq. (11)with S =0.5. Panels (b) and (c) of Fig. 6 demonstratethat SDP-ROHF yields the correct value of S ( S + 1) atall O-H distances, and the norm of the errors definedby Eq. (A6) are zero for the entire curve as well. Thecorresponding data for SDP-UHF show clear deviationsfrom the respective values for a pure spin state, asexpected. Acknowledgments:
This work was supported as part of the Center for Ac-tinide Science and Technology (CAST), an Energy Fron-tier Research Center funded by the U.S. Department ofEnergy, Office of Science, Basic Energy Sciences underAward No. de-sc0016568. [1] J. ˇC´ıˇzek, J. Chem. Phys. , 4256 (1966).[2] G. D. Purvis and R. J. Bartlett, J. Chem. Phys. , 1910(1982). [3] K. Raghavachari, G. W. Trucks, J. A. Pople, andM. Head-Gordon, Chem. Phys. Lett. , 479 (1989).[4] R. J. Bartlett and M. Musia l, Rev. Mod. Phys. , 291 (2007).[5] C. C. J. Roothaan, Rev. Mod. Phys. , 69 (1951).[6] G. G. Hall, Proc. Royal Soc. A , 541 (1951).[7] P. Pulay, Chem. Phys. Lett. , 393 (1980).[8] P. Pulay, J. Comput. Chem. , 556 (1982).[9] K. Yasuda, J. Comput. Chem. , 334 (2008).[10] I. S. Ufimtsev and T. J. Martnez, J. Chem. Theory Com-put. , 222 (2008).[11] N. Luehr, I. S. Ufimtsev, and T. J. Martnez, J. Chem.Theory Comput. , 949 (2011).[12] J. Kalinowski, F. Wennmohs, and F. Neese, J. Chem.Theory Comput. , 3160 (2017).[13] P.-O. L¨owdin, Phys. Rev. , 1490 (1955).[14] R. McWeeny, Proc. Royal Soc. A , 496 (1956).[15] R. McWeeny, Phys. Rev. , 1528 (1959).[16] R. McWeeny, Phys. Rev. , 1028 (1962).[17] R. McWeeny, Rev. Mod. Phys. , 335 (1960).[18] F. Mauri, G. Galli, and R. Car, Phys. Rev. B , 9973(1993).[19] X.-P. Li, R. W. Nunes, and D. Vanderbilt, Phys. Rev. B , 10891 (1993).[20] E. B. Stechel, A. R. Williams, and P. J. Feibelman, Phys.Rev. B , 10088 (1994).[21] J. Kussmann, M. Beer, and C. Ochsenfeld, WIREs Com-put. Mol. Sci. , 614 (2013), ISSN 1759-0884.[22] A. J. Coleman, Rev. Mod. Phys. , 668 (1963).[23] S. Veeraraghavan and D. A. Mazziotti, Phys. Rev. A ,010502 (2014).[24] S. Veeraraghavan and D. A. Mazziotti, J. Chem. Phys. , 124106 (2014).[25] S. Veeraraghavan and D. A. Mazziotti, Phys. Rev. A ,022512 (2015).[26] A. J. Coleman and I. Absar, Int. J. Quantum Chem. ,1279 (1980).[27] E. P´erez-Romero, L. M. Tel, and C. Valdemoro, Int. J.Quantum Chem. , 55 (1997).[28] G. Gidofalvi and D. A. Mazziotti, Phys. Rev. A ,052505 (2005).[29] H. Fukutome, Int. J. Quantum Chem. , 955 (1981).[30] J. Stuber and J. Paldus, Fundamental World of QuantumChemistry, A Tribute Volume to the Memory of Per-OlovL¨owdin , 67 (2003).[31] C. A. Jim´enez-Hoyos, T. M. Henderson, and G. E. Scuse-ria, J. Chem. Theory Comput. , 2667 (2011).[32] J. G. Valatin, Phys. Rev. , 1012 (1961). [33] A. K. Kerman and A. Klein, Phys. Rev. , 1326 (1963).[34] M. M. Mestechkin and G. E. Whyman, Int. J. QuantumChem. , 45 (1974).[35] D. A. Mazziotti, Phys. Rev. Lett. , 213001 (2004).[36] D. A. Mazziotti, ESAIM, Math. Model. Numer. Anal. , 249 (2007).[37] D. W. Small, E. J. Sundstrom, and M. Head-Gordon, J.Chem. Phys. , 024104 (2015).[38] G. D. Purvis, R. Shepard, F. B. Brown, and R. J.Bartlett, Int. J. Quantum Chem. , 835 (1983).[39] U. S. Mahapatra, B. Datta, B. Bandyopadhyay, andD. Mukherjee, Adv. Quantum Chem. , 163 (1998).[40] U. S. Mahapatra, B. Datta, and D. Mukherjee, J. Chem.Phys. , 6171 (1999).[41] F. A. Evangelista, M. Hanauer, A. K¨ohn, and J. Gauss,J. Chem. Phys. , 204108 (2012).[42] D. I. Lyakh, M. Musia l, V. F. Lotrich, and R. J. Bartlett,Chem. Rev , 182 (2012).[43] J. Brabec, H. J. van Dam, J. Pittner, and K. Kowalski,J. Chem. Phys. , 124102 (2012).[44] Z. Chen and M. R. Hoffmann, J. Chem. Phys. ,014108 (2012).[45] J. P. Coe, D. J. Taylor, and M. J. Paterson, J. Chem.Phys. , 194111 (2012).[46] O. Demel, S. Kedˇzuch, M. ˇSvaˇna, S. Ten-No, J. Pittner,and J. Noga, Phys. Chem. Chem. Phys. , 4753 (2012).[47] C. A. Jim´enez-Hoyos, R. Rodr´ıguez-Guzm´an, and G. E.Scuseria, J. Chem. Phys. , 224110 (2013).[48] R. M. Parrish, L. A. Burns, D. G. A. Smith, A. C. Sim-monett, A. E. DePrince, E. G. Hohenstein, U. Bozkaya,A. Y. Sokolov, R. Di Remigio, R. M. Richard, et al., J.Chem. Theory Comput. , 3185 (2017).[49] X. Li and J. Paldus, J. Chem. Phys. , 084110 (2009).[50] A. J. W. Thom and M. Head-Gordon, Phys. Rev. Lett. , 193001 (2008).[51] C. A. Jim´enez-Hoyos, T. M. Henderson, T. Tsuchimochi,and G. E. Scuseria, J. Chem. Phys. , 164109 (2012).[52] M. Piris, in Many-body approaches at different scales: Atribute to Norman H. March on the occasion of his 90thbirthday , edited by G. G. N. Angilella and C. Amovilli(Springer, 2018), chap. 22, pp. 283–300.[53] S. M. Valone, J. Chem. Phys. , 1344 (1980).[54] H. van Aggelen, B. Verstichel, P. Bultinck, D. V. Neck,and P. W. Ayers., J. Chem. Phys.136