Convex Hartree–Fock Reformulations

• Convex Hartree–Fock is an electronic structure method that replaces nonconvex idempotency with convex admissibility conditions on the one-particle density matrix.
• SDP formulations provide rigorous lower and upper energy bounds and global optimality certificates, with benchmark tests showing near-quantitative agreement in systems like CN and Cr2.
• Advanced techniques including CVX-HF and occupation-number polytope methods address symmetry breaking and extend HF to efficient multiconfigurational approximations.

Searching arXiv for the specified paper and closely related work on convex Hartree–Fock formulations. arXiv search query: id:([1712.07680](/papers/1712.07680)) OR id:([2508.21453](/papers/2508.21453)) OR id:([1404.5217](/papers/1404.5217)) OR id:([1711.09129](/papers/1711.09129)) OR id:([2209.10189](/papers/2209.10189)) Convex Hartree–Fock denotes a family of Hartree–Fock reformulations in which convex structure is imposed either on the admissible one-particle variables, on lifted reduced-density-matrix representations, or on the local orbital-rotation landscape near instabilities. In the literature represented by Lieb’s variational principle, semidefinite-programming (SDP) Hartree–Fock, generalized-Pauli-constraint extensions, and the recent CVX-HF treatment of conical intersections, the common objective is to replace or control the nonconvex aspects of the conventional self-consistent-field problem without abandoning Hartree–Fock as the underlying reference theory (Bach, 2022).

1. Convex admissible sets and the Hartree–Fock variational principle

A central starting point is the one-particle reduced density matrix (1-RDM). For a Slater determinant

$\Phi(f_1,\dots,f_N)=f_1\wedge\cdots\wedge f_N,$

the associated 1-RDM is the rank-$N$ orthogonal projection

$\gamma_f=\sum_{i=1}^N |f_i\rangle\langle f_i|, \qquad \gamma_f^2=\gamma_f,\;0\le \gamma_f\le 1,\;\Tr \gamma_f=N.$

Lieb’s variational principle identifies the Hartree–Fock ground-state energy with the infimum of the Hartree–Fock functional over the convex set

$\mathcal K_N=\Bigl\{\gamma=\gamma^*\in\mathcal L^1(\mathfrak h)\;\Big|\;0\le \gamma\le 1,\;\Tr\gamma=N\Bigr\},$

rather than only over rank-$N$ projections (Bach, 2022).

For Coulomb systems, the extended Hartree–Fock functional is

$\mathcal E_{\rm HF}[\gamma] = \Tr(h\gamma) +\frac12\iint V(x,y)\bigl[\rho_\gamma(x)\rho_\gamma(y)-|\gamma(x,y)|^2\bigr]\,dx\,dy,$

with $h=-\tfrac12\Delta-\sum_k Z_k|x-R_k|^{-1}$, $V(x,y)=|x-y|^{-1}$, and $\rho_\gamma(x)=\gamma(x,x)$. The exact statement is that

$$E_{\rm HF}(N)=\inf_{\gamma\in\mathcal K_N}\mathcal E_{\rm HF}[\gamma].$$

This does not mean that the Hartree–Fock objective becomes globally convex in $N$0. The same review emphasizes that the direct Hartree term is convex in $N$1, the exchange term is concave when $N$2, and therefore $N$3 is a difference-of-convex functional. A common misconception is thus to identify “convex Hartree–Fock” with a convex energy functional; the exact reformulation instead places the search over a convex admissible set while preserving a nontrivial exchange contribution (Bach, 2022).

This convex-set viewpoint is foundational for later SDP formulations. It isolates the operator inequalities $N$4 and $N$5 as the essential admissibility conditions and makes explicit that idempotency is not needed to characterize the Hartree–Fock infimum at the variational level, even though it reappears in constructive formulations of single-determinant solutions.

2. Semidefinite programming formulations and global Hartree–Fock optimization

In finite orbital bases, SDP-based Hartree–Fock casts the problem directly in terms of reduced density matrices. One 1-RDM formulation writes the Hartree–Fock energy as

$N$6

where $N$7, $N$8 are the usual two-electron Coulomb integrals, and $N$9 in an orthonormal spin-orbital basis (Nascimento et al., 2017).

Ensemble $\gamma_f=\sum_{i=1}^N |f_i\rangle\langle f_i|, \qquad \gamma_f^2=\gamma_f,\;0\le \gamma_f\le 1,\;\Tr \gamma_f=N.$0-representability is enforced by introducing a one-hole matrix $\gamma_f=\sum_{i=1}^N |f_i\rangle\langle f_i|, \qquad \gamma_f^2=\gamma_f,\;0\le \gamma_f\le 1,\;\Tr \gamma_f=N.$1 and imposing

$\gamma_f=\sum_{i=1}^N |f_i\rangle\langle f_i|, \qquad \gamma_f^2=\gamma_f,\;0\le \gamma_f\le 1,\;\Tr \gamma_f=N.$2

From these relations it follows that $\gamma_f=\sum_{i=1}^N |f_i\rangle\langle f_i|, \qquad \gamma_f^2=\gamma_f,\;0\le \gamma_f\le 1,\;\Tr \gamma_f=N.$3 as matrices and hence all natural-orbital eigenvalues lie in $\gamma_f=\sum_{i=1}^N |f_i\rangle\langle f_i|, \qquad \gamma_f^2=\gamma_f,\;0\le \gamma_f\le 1,\;\Tr \gamma_f=N.$4. To recover exactly the single-determinant Hartree–Fock solution one must additionally enforce the nonconvex condition

$\gamma_f=\sum_{i=1}^N |f_i\rangle\langle f_i|, \qquad \gamma_f^2=\gamma_f,\;0\le \gamma_f\le 1,\;\Tr \gamma_f=N.$5

which forces $\gamma_f=\sum_{i=1}^N |f_i\rangle\langle f_i|, \qquad \gamma_f^2=\gamma_f,\;0\le \gamma_f\le 1,\;\Tr \gamma_f=N.$6 (Nascimento et al., 2017).

A complementary SDP construction introduces an auxiliary matrix $\gamma_f=\sum_{i=1}^N |f_i\rangle\langle f_i|, \qquad \gamma_f^2=\gamma_f,\;0\le \gamma_f\le 1,\;\Tr \gamma_f=N.$7 with entries $\gamma_f=\sum_{i=1}^N |f_i\rangle\langle f_i|, \qquad \gamma_f^2=\gamma_f,\;0\le \gamma_f\le 1,\;\Tr \gamma_f=N.$8, so that the Hartree–Fock energy becomes linear in $\gamma_f=\sum_{i=1}^N |f_i\rangle\langle f_i|, \qquad \gamma_f^2=\gamma_f,\;0\le \gamma_f\le 1,\;\Tr \gamma_f=N.$9. The resulting convex “lb-SDP” yields a rigorous lower bound $\mathcal K_N=\Bigl\{\gamma=\gamma^*\in\mathcal L^1(\mathfrak h)\;\Big|\;0\le \gamma\le 1,\;\Tr\gamma=N\Bigr\},$0, while a rank-constrained SDP with $\mathcal K_N=\Bigl\{\gamma=\gamma^*\in\mathcal L^1(\mathfrak h)\;\Big|\;0\le \gamma\le 1,\;\Tr\gamma=N\Bigr\},$1 yields an upper bound $\mathcal K_N=\Bigl\{\gamma=\gamma^*\in\mathcal L^1(\mathfrak h)\;\Big|\;0\le \gamma\le 1,\;\Tr\gamma=N\Bigr\},$2. Equality of the upper- and lower-bound energies guarantees that the computed solution is the globally optimal solution of Hartree–Fock theory (Veeraraghavan et al., 2014).

The ROHF applications reported for this framework are explicitly targeted at strongly correlated systems. For $\mathcal K_N=\Bigl\{\gamma=\gamma^*\in\mathcal L^1(\mathfrak h)\;\Big|\;0\le \gamma\le 1,\;\Tr\gamma=N\Bigr\},$3 in 6-31G$\mathcal K_N=\Bigl\{\gamma=\gamma^*\in\mathcal L^1(\mathfrak h)\;\Big|\;0\le \gamma\le 1,\;\Tr\gamma=N\Bigr\},$4, lb-SDP is within $\mathcal K_N=\Bigl\{\gamma=\gamma^*\in\mathcal L^1(\mathfrak h)\;\Big|\;0\le \gamma\le 1,\;\Tr\gamma=N\Bigr\},$5 a.u. of rc-SDP over most of the curve; for CN in cc-pVDZ, lb-SDP lies within $\mathcal K_N=\Bigl\{\gamma=\gamma^*\in\mathcal L^1(\mathfrak h)\;\Big|\;0\le \gamma\le 1,\;\Tr\gamma=N\Bigr\},$6 a.u.; for $\mathcal K_N=\Bigl\{\gamma=\gamma^*\in\mathcal L^1(\mathfrak h)\;\Big|\;0\le \gamma\le 1,\;\Tr\gamma=N\Bigr\},$7 in a TZV basis, the rc-SDP branch beyond $\mathcal K_N=\Bigl\{\gamma=\gamma^*\in\mathcal L^1(\mathfrak h)\;\Big|\;0\le \gamma\le 1,\;\Tr\gamma=N\Bigr\},$8 Å is up to $\mathcal K_N=\Bigl\{\gamma=\gamma^*\in\mathcal L^1(\mathfrak h)\;\Big|\;0\le \gamma\le 1,\;\Tr\gamma=N\Bigr\},$9 a.u. lower than the DIIS solution and lb-SDP confirms this branch as globally optimal to within $N$0 a.u.; for $N$1 bend in cc-pVDZ, lb-SDP certifies the global minimum within $N$2 a.u. The same study states that lb-SDP typically scales like $N$3, rc-SDP like $N$4, and that no symmetry guesses are needed (Veeraraghavan et al., 2014).

Within this SDP literature, “convex Hartree–Fock” therefore often refers not to a single algorithm but to a convex relaxation architecture: a convex feasible region supplies lower bounds and certificates, while nonconvex rank conditions recover the exact single-determinant problem. This suggests a precise distinction between convex admissibility and nonconvex determinantal representability.

3. Spin and spatial symmetry breaking in SDP-based Hartree–Fock

The paper on spatial and spin symmetry breaking develops an equivalent optimization over positive semidefinite 1-RDMs that retains the nonconvexity of the Hartree–Fock energy expression while allowing ensemble spin-state conditions to be imposed directly on $N$5 (Nascimento et al., 2017). When a well-defined third component of spin $N$6 is enforced, $N$7 is written in spin blocks,

$N$8

with additional constraints

$N$9

and

$\mathcal E_{\rm HF}[\gamma] = \Tr(h\gamma) +\frac12\iint V(x,y)\bigl[\rho_\gamma(x)\rho_\gamma(y)-|\gamma(x,y)|^2\bigr]\,dx\,dy,$0

If one further requires $\mathcal E_{\rm HF}[\gamma] = \Tr(h\gamma) +\frac12\iint V(x,y)\bigl[\rho_\gamma(x)\rho_\gamma(y)-|\gamma(x,y)|^2\bigr]\,dx\,dy,$1, the nonlinear pure-spin condition is

$\mathcal E_{\rm HF}[\gamma] = \Tr(h\gamma) +\frac12\iint V(x,y)\bigl[\rho_\gamma(x)\rho_\gamma(y)-|\gamma(x,y)|^2\bigr]\,dx\,dy,$2

All of these constraints are convex in $\mathcal E_{\rm HF}[\gamma] = \Tr(h\gamma) +\frac12\iint V(x,y)\bigl[\rho_\gamma(x)\rho_\gamma(y)-|\gamma(x,y)|^2\bigr]\,dx\,dy,$3 except the last one, which is bilinear in the spin blocks (Nascimento et al., 2017).

The treatment of spatial symmetry is especially significant. In a conventional point-group-adapted calculation one would require $\mathcal E_{\rm HF}[\gamma] = \Tr(h\gamma) +\frac12\iint V(x,y)\bigl[\rho_\gamma(x)\rho_\gamma(y)-|\gamma(x,y)|^2\bigr]\,dx\,dy,$4 to be block diagonal according to the irreducible representations of the molecular point group. In the SDP framework, one may simply omit that linear block-diagonality constraint. The global minimizer of the convex region defined by the ensemble conditions may then spill between irreps and break spatial symmetry spontaneously (Nascimento et al., 2017).

The Be–H$\mathcal E_{\rm HF}[\gamma] = \Tr(h\gamma) +\frac12\iint V(x,y)\bigl[\rho_\gamma(x)\rho_\gamma(y)-|\gamma(x,y)|^2\bigr]\,dx\,dy,$5 insertion path is the exemplary case. Along this path, the two distinct $\mathcal E_{\rm HF}[\gamma] = \Tr(h\gamma) +\frac12\iint V(x,y)\bigl[\rho_\gamma(x)\rho_\gamma(y)-|\gamma(x,y)|^2\bigr]\,dx\,dy,$6 RHF solutions with dominant configurations $\mathcal E_{\rm HF}[\gamma] = \Tr(h\gamma) +\frac12\iint V(x,y)\bigl[\rho_\gamma(x)\rho_\gamma(y)-|\gamma(x,y)|^2\bigr]\,dx\,dy,$7 or $\mathcal E_{\rm HF}[\gamma] = \Tr(h\gamma) +\frac12\iint V(x,y)\bigl[\rho_\gamma(x)\rho_\gamma(y)-|\gamma(x,y)|^2\bigr]\,dx\,dy,$8 become degenerate. Enforcing only $\mathcal E_{\rm HF}[\gamma] = \Tr(h\gamma) +\frac12\iint V(x,y)\bigl[\rho_\gamma(x)\rho_\gamma(y)-|\gamma(x,y)|^2\bigr]\,dx\,dy,$9 and $h=-\tfrac12\Delta-\sum_k Z_k|x-R_k|^{-1}$0 but not $h=-\tfrac12\Delta-\sum_k Z_k|x-R_k|^{-1}$1 symmetry, the SDP solver finds a single smooth $h=-\tfrac12\Delta-\sum_k Z_k|x-R_k|^{-1}$2 that mixes $h=-\tfrac12\Delta-\sum_k Z_k|x-R_k|^{-1}$3 and $h=-\tfrac12\Delta-\sum_k Z_k|x-R_k|^{-1}$4 character continuously, thereby removing the cusp in the RHF potential-energy curve (Nascimento et al., 2017). The same work also demonstrates numerically that, upon relaxation of $h=-\tfrac12\Delta-\sum_k Z_k|x-R_k|^{-1}$5 and $h=-\tfrac12\Delta-\sum_k Z_k|x-R_k|^{-1}$6 symmetry constraints, the RDM-based approach is equivalent to real-valued generalized Hartree–Fock theory.

From the algorithmic standpoint, this framework updates $h=-\tfrac12\Delta-\sum_k Z_k|x-R_k|^{-1}$7 directly via semidefinite optimization, or its factored variant $h=-\tfrac12\Delta-\sum_k Z_k|x-R_k|^{-1}$8, and therefore does not require diagonalization of the Fock matrix in the inner loop. The stated practical implication is that linearly-scaling techniques can be more readily incorporated (Nascimento et al., 2017).

4. CVX-HF and ground-state conical intersections

A distinct recent use of the label “Convex Hartree–Fock” appears in the framework proposed for ground-state conical intersections. Here the problem is not global optimization over 1-RDMs but the loss of positive definiteness of the Hartree–Fock orbital-rotation Hessian near a would-be conical intersection (Rossi et al., 29 Aug 2025). Standard Hartree–Fock writes the energy as

$h=-\tfrac12\Delta-\sum_k Z_k|x-R_k|^{-1}$9

and requires stationarity with respect to the orbital-rotation generator $V(x,y)=|x-y|^{-1}$0. Near a conical intersection, the local quadratic form develops a zero-curvature direction in orbital-rotation space, so Hessian eigenvalues approach or cross zero, multiple HF solutions appear, energies become discontinuous, and the topology of the $V(x,y)=|x-y|^{-1}$1–$V(x,y)=|x-y|^{-1}$2 seam is wrong (Rossi et al., 29 Aug 2025).

CVX-HF restores a strictly convex optimization landscape for the ground-state determinant by projecting out the orbital-rotation directions along which the Hessian first loses positive definiteness. With orbital gradient and Hessian

$V(x,y)=|x-y|^{-1}$3

one solves

$V(x,y)=|x-y|^{-1}$4

orders the eigenvalues $V(x,y)=|x-y|^{-1}$5, and constructs

$V(x,y)=|x-y|^{-1}$6

The projected gradient is $V(x,y)=|x-y|^{-1}$7, and orbital updates are restricted to the image of $V(x,y)=|x-y|^{-1}$8. In the projected subspace, $V(x,y)=|x-y|^{-1}$9 is strictly positive definite when the problematic low-curvature mode is removed (Rossi et al., 29 Aug 2025).

The self-consistent algorithm begins from a simple SAD guess $\rho_\gamma(x)=\gamma(x,x)$0 and $\rho_\gamma(x)=\gamma(x,x)$1. At each iteration it builds the Hessian, extracts the lowest mode, forms $\rho_\gamma(x)=\gamma(x,x)$2, solves

$\rho_\gamma(x)=\gamma(x,x)$3

in the Krylov/subspace sense, and updates

$\rho_\gamma(x)=\gamma(x,x)$4

At convergence one obtains a unique convexified HF determinant with no instabilities in the remaining subspace (Rossi et al., 29 Aug 2025).

The removed mode is then reintroduced through a final small Hamiltonian diagonalization in a basis $\rho_\gamma(x)=\gamma(x,x)$5, where

$\rho_\gamma(x)=\gamma(x,x)$6

Solving

$\rho_\gamma(x)=\gamma(x,x)$7

restores multistate behavior, and because the HF–$\rho_\gamma(x)=\gamma(x,x)$8 coupling is reintroduced, $\rho_\gamma(x)=\gamma(x,x)$9 and $$E_{\rm HF}(N)=\inf_{\gamma\in\mathcal K_N}\mathcal E_{\rm HF}[\gamma].$$0 repel, creating the characteristic double cone (Rossi et al., 29 Aug 2025).

The reported benchmarks include ammonia in aug-cc-pVDZ, a one-dimensional N–H stretch at $$E_{\rm HF}(N)=\inf_{\gamma\in\mathcal K_N}\mathcal E_{\rm HF}[\gamma].$$1, 2,4-cyclohexadien-1-ylamine in cc-pVDZ, and the GFP chromophore HBDI$$E_{\rm HF}(N)=\inf_{\gamma\in\mathcal K_N}\mathcal E_{\rm HF}[\gamma].$$2 in 6-31G$$E_{\rm HF}(N)=\inf_{\gamma\in\mathcal K_N}\mathcal E_{\rm HF}[\gamma].$$3. Across these tests, TDHF-TDA shows negative excitation energies, multiple HF solution branches, nonconvergence, or distorted $$E_{\rm HF}(N)=\inf_{\gamma\in\mathcal K_N}\mathcal E_{\rm HF}[\gamma].$$4 surfaces, whereas CVX-HF yields smooth surfaces meeting at a single point for NH$$E_{\rm HF}(N)=\inf_{\gamma\in\mathcal K_N}\mathcal E_{\rm HF}[\gamma].$$5, a single avoided crossing at approximately $$E_{\rm HF}(N)=\inf_{\gamma\in\mathcal K_N}\mathcal E_{\rm HF}[\gamma].$$6 Å in the one-dimensional scan, a proper two-dimensional conical seam for 2,4-cyclohexadien-1-ylamine, and the expected double cone centered on the known P90 MECI geometry for HBDI$$E_{\rm HF}(N)=\inf_{\gamma\in\mathcal K_N}\mathcal E_{\rm HF}[\gamma].$$7 (Rossi et al., 29 Aug 2025).

This formulation remains a single-reference method, so absolute excitation energies still carry the usual HF error. The same paper states, however, that the orbitals are free from HF instabilities and may be used in correlated post-HF methods, that the framework is orbital-invariant and size-intensive when projecting a fixed number of states, and that size-extensivity can be fully restored by including products of projected modes in $$E_{\rm HF}(N)=\inf_{\gamma\in\mathcal K_N}\mathcal E_{\rm HF}[\gamma].$$8 (Rossi et al., 29 Aug 2025).

5. Occupation-number polytopes and the natural extension of Hartree–Fock

Another line of work uses convexity at the level of natural occupation numbers rather than directly through SDP constraints. Beyond Pauli’s exclusion principle $$E_{\rm HF}(N)=\inf_{\gamma\in\mathcal K_N}\mathcal E_{\rm HF}[\gamma].$$9, the natural occupation numbers $N$00 of a pure $N$01-fermion state satisfy generalized Pauli constraints

$N$02

which, together with $N$03 and the ordering $N$04, define a convex polytope

$N$05

This polytope is a proper subset of the Pauli simplex whenever nontrivial generalized Pauli constraints are present (Benavides-Riveros et al., 2017).

When the exact occupation vector lies on a facet

$N$06

one has

$N$07

In a Slater-determinant expansion

$N$08

this yields the selection rule

$N$09

For the Borland–Dennis case $N$10, pinning of the nontrivial generalized Pauli constraint reduces the state to

$N$11

namely three Slater determinants only (Benavides-Riveros et al., 2017).

The resulting variational ansatz is formulated as a convex optimization over $N$12: $N$13 where, in the natural-orbital basis,

$N$14

and

$N$15

The energy is

$N$16

This is presented as a natural multiconfigurational generalization of Hartree–Fock in which the geometry of the occupation-number polytope controls the size of the active configurational subspace (Benavides-Riveros et al., 2017).

The same framework establishes geometric bounds on the recovered correlation energy. With $N$17 the exact ground-state energy, $N$18 the Hartree–Fock energy, $N$19 the facet-pinned ansatz energy, $N$20, and $N$21, one has

$N$22

where $N$23 is the $N$24-distance to the Hartree–Fock vertex. Whenever the exact occupation numbers are quasi-pinned so that $N$25, the ansatz recovers nearly all of the true correlation energy (Benavides-Riveros et al., 2017).

6. Related convex mean-field models, guarantees, and limitations

Reduced Hartree–Fock (rHF) provides an important neighboring model in which convexity is exact at the density-functional level. In periodic notation, the rHF energy is

$N$26

with $N$27, $N$28, and $N$29, so the entire nonlinear piece $N$30 is strictly convex in $N$31 (Bordignon et al., 2024). This model is not identical to full Hartree–Fock because the exchange term is removed, but it shows how convex mean-field structure supports fully guaranteed and computable a posteriori energy bounds.

For periodic Kohn–Sham equations with convex density functionals, the bound derived in this setting is

$N$32

and this error can be decomposed into discretization and SCF contributions,

$N$33

The reported numerical illustrations include a Silicon crystal and a Hydrogen Fluoride molecule simulated with the rHF model (Bordignon et al., 2024). A plausible implication is that convex Hartree–Fock methodologies sit naturally within a broader program of certified nonlinear mean-field computation.

Across these formulations, the principal limitations are explicit in the source literature. The exact Hartree–Fock functional on the convex set $N$34 remains a difference-of-convex functional rather than a globally convex one (Bach, 2022). SDP relaxations recover exact single-determinant Hartree–Fock only when the requisite rank conditions are imposed or when upper and lower bounds coincide (Veeraraghavan et al., 2014). In the 1-RDM symmetry-breaking formulation, the pure-spin constraint is bilinear and therefore nonconvex (Nascimento et al., 2017). In CVX-HF for conical intersections, the method remains single-reference and therefore does not remove the usual Hartree–Fock error in absolute excitation energies (Rossi et al., 29 Aug 2025).

Taken together, these results show that convex Hartree–Fock is best understood as a family of strategies for isolating, relocating, or regularizing Hartree–Fock nonconvexity. In one branch, convex admissible sets and SDP liftings supply rigorous bounds, global-optimality certificates, and a flexible treatment of broken symmetry. In another, the Hessian-based CVX-HF construction restores a strictly convex local optimization landscape near conical intersections while recovering multistate topology through a final diagonalization. In a third, occupation-number geometry converts facet pinning in a convex polytope into a controlled multiconfigurational extension of the Hartree–Fock ansatz.