Convex Hartree-Fock theory: A simple framework for ground state conical intersections

TL;DR

Convex Hartree-Fock introduces a modified approach that efficiently models ground state conical intersections, overcoming limitations of traditional methods by optimizing within a tailored subspace and enabling accurate excited state calculations.

Contribution

This work presents Convex Hartree-Fock, a novel framework that improves modeling of conical intersections by optimizing the reference state within a specific subspace, reducing computational costs.

Findings

Successfully models ground state conical intersections.

Outperforms traditional Hartree-Fock and TDHF in benchmarks.

Validates approach across multiple test cases.

Abstract

Accurate modeling of conical intersections is crucial in nonadiabatic molecular dynamics, as these features govern processes such as radiationless transitions and photochemical reactions. Conventional electronic structure methods, including Hartree-Fock, density functional theory, and their time-dependent extensions, struggle in this regime. Due to their single reference nature and separate treatment of ground and excited states, they fail to capture ground state intersections. Multiconfigurational approaches overcome these limitations, but at a prohibitive computational cost. In this work, we propose a modified Hartree-Fock framework, referred to as Convex Hartree-Fock, that optimizes the reference within a tailored subspace by removing projections along selected Hessian eigenvectors. The ground and excited states are then obtained through subsequent Hamiltonian diagonalization. We validate the approach across several test cases and benchmark its performance against time-dependent Hartree-Fock within the Tamm-Dancoff approximation.

Equations

$$|\mathrm{HF}⟩=\exp (\sum _{ai}{\kappa }_{ai}{E}_{ai}^{-})|{\Phi }_0⟩,$$

$${\mathcal{E}}_{\text{HF}}=⟨\mathrm{HF}|H|\mathrm{HF}⟩,$$

$$\mathcal{E}(\mathbf{\gamma })=⟨\mathrm{HF}|\exp (-\gamma )H\exp (\gamma )|\mathrm{HF}⟩,$$

$$G_{ai}^{(0)}(0)$$

$$G_{ai,bj}^{(1)}(0)$$

$$\left(\begin{array}{cc}\mathbf{A}& \mathbf{B}\\ \mathbf{B}& \mathbf{A}\end{array}\right)\left(\begin{array}{c}\mathbf{X}\\ \mathbf{Y}\end{array}\right)=\omega \left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)\left(\begin{array}{c}\mathbf{X}\\ \mathbf{Y}\end{array}\right),$$

$$A_{ia,jb}={\delta }_{ij}{\delta }_{ab}({\epsilon }_a-{\epsilon }_i)+2g_{aijb}-g_{abji}$$

$$B_{ia,jb}=2g_{aibj}-g_{ajbi},$$

$$G_{ai}^{(0)}(\mathbf{\gamma })=G_{ai}^{(0)}(0)+\sum _{bj}{G}_{ai,bj}^{(1)}(0){\gamma }_{bj}+\dots =0.$$

$$G_n^{(0)}(\mathbf{\gamma })=G_n^{(0)}(0)+{\lambda }_n{\gamma }_n+\dots =0.$$

$${\tilde{\mathbf{G}}}^{(0)}={\mathbf{G}}^{(0)}-{\mathbf{r}}_1({\mathbf{r}}_1^T{\mathbf{G}}^{(0)}).$$

$$\kappa =\sum _{ai}{\kappa }_{ai}{E}_{ai}^{-}=\sum _n{\kappa }_n{R}_n$$

$${\tilde{\kappa }}_{ai}={\kappa }_{ai}-r_{ai,1}\sum _{bj}{r}_{bj,1}{\kappa }_{bj},$$

$${\tilde{\mathbf{G}}}^{(0)}=0$$

$${\mathbf{G}}^{(1)}{\mathbf{r}}_1={\lambda }_1{\mathbf{r}}_1,$$

$${\mathbf{H}}^{\mathrm{FS}}$$

$${\mathbf{H}}^{\mathrm{FS}}{\mathbf{x}}_n={\mathcal{E}}_n{\mathbf{S}}^{\mathrm{FS}}{\mathbf{x}}_n$$

$${\mathbf{S}}^{\mathrm{FS}}$$

$$|{\psi }_A⟩=|{\text{HF}}_A⟩x_0^A+|R_A⟩x_1^A$$

$$|{\psi }_B⟩=|{\text{HF}}_B⟩x_0^B+|R_B⟩x_1^B.$$

$$|\psi ⟩=|\text{HF}⟩x_0+|R_1⟩x_1+|R_2⟩x_2,$$

$$|\text{HF}⟩=|{\text{HF}}_A⟩|{\text{HF}}_B⟩$$

$$|R_1⟩=|R_A⟩|{\text{HF}}_B⟩$$

$$|R_2⟩=|{\text{HF}}_A⟩|R_B⟩.$$

Full text

Introduction

The electronic structure of molecules in regions where two or more states approach degeneracy presents one of the most complex challenges in quantum chemistry. In particular, the behavior of systems near ground state conical intersections demands a theoretical description capable of capturing pronounced multireference character. Traditional single-determinant frameworks like Hartree–Fock (HF) and density functional theory (DFT), despite their widespread use and computational efficiency, are fundamentally inadequate for this task. Both Hartree-Fock and conventional Kohn-Sham DFT rely on a mean-field approximation that implicitly assumes a dominant electronic configuration. This assumption breaks down in the vicinity of such conical intersections, where the ground state wave function becomes strongly mixed and cannot be described by a single Slater determinant. The results are discontinuities in the energies, incorrect topology of the potential energy surfaces, and convergence failures [1]. For DFT, the problem is further intensified by the fact that most exchange–correlation functionals are developed for a non-degenerate ground state, leading to unreliable predictions when this is no longer the case [2, 3].

The stability of Hartree-Fock solutions has been studied in great detail previously [4, 5, 6]. Bifurcations of the solutions are well known to appear frequently on potential energy surfaces, for instance, the Coulson-Fischer point in the hydrogen molecule for an unrestricted Hartree-Fock wave function [7, 8]. The non-analytical behavior has been discussed by Čížek and Paldus [6] and was recently shown to be particularly pronounced in areas around ground state conical intersections [9]. These limitations directly impact the accuracy of simulations in photochemistry and ultrafast spectroscopy. Due to the critical importance of ground state conical intersections, a range of computational solutions has been developed over the years to address these degeneracies within time-dependent mean-field methods. One example is spin-flip TDDFT [10], which redefines the reference state to access crossings between ground and excited states by allowing spin-changing excitations. Other approaches include modified exchange-correlation kernels to account for double excitation effects [11], and TDDFT-1D [12], where specific excited configurations are incorporated into the eigenvalue problem. Configuration interaction-corrected DFT (CIC-TDA [13]) restores the correct dimensionality of conical intersections by augmenting the Tamm–Dancoff approximation (TDA) with a perturbative configuration interaction scheme. Further advances include orthogonality constrained DFT (OCDFT [14]), multiconfigurational short-range DFT [15], and ensemble-based DFT methods adapted to describe excited states with near-degeneracies (SI-SA-REKS [16]). In addition, mode-following techniques [17] have been employed to directly locate excited state saddle points by following specific eigenvectors of the Hessian. Most recently, phase-space electronic structure theory [18] reformulates electronic structure in terms of phase-space variables, offering new insights into conical intersections by explicitly incorporating electronic momentum. These and other methods to address conical intersection problems have recently been reviewed by Matsika [19].

To our knowledge, no general procedure to resolve the aforementioned complications has been developed by correcting the wave function parameterization in Hartree-Fock theory. In this work, we introduce the Convex Hartree–Fock (CVX-HF) method, in which the stationarity conditions and the Hessian eigenvalue equations are solved self-consistently within a subspace where the Hessian is positive definite. This framework offers excellent convergence properties and introduces the necessary coupling elements in the Hamiltonian matrix to yield eigenstates that correctly capture conical intersections and the geometric phase effect. Since TDA is typically used to avoid complex excitation energies arising from the non-Hermitian TDHF eigenvalue problem [20], we compare the performance of CVX-HF to TDHF with TDA in a variety of molecular systems.

Results

Convex Hartree-Fock theory

We start from a reference determinant $|\Phi_{0}\rangle$ constructed from orbitals obtained from the superposition of atomic densities (SAD) [21, 22]. The Hartree-Fock wave function is parameterized in terms of a single global orbital rotation, such that

[TABLE]

where $E_{ai}^{-}=E_{ai}-E_{ia}$, $E_{ai}=\sum_{\sigma}a^{\dagger}_{a\sigma}a_{i\sigma}$, and $\kappa_{ai}$ are the orbital rotation parameters. We use $a,b$ to label virtual orbitals and $i,j$ for occupied orbitals. The Hartree-Fock energy $\mathcal{E}_{\text{HF}}$ is obtained as the expectation value of the electronic Hamiltonian $H$,

[TABLE]

where we assume $|\mathrm{HF}\rangle$ is normalized. When minimizing the energy with respect to the $\kappa_{ai}$ parameters, we need to calculate the electronic gradient and Hessian. These are conveniently obtained by considering the Hartree-Fock energy function

[TABLE]

where $\exp(\gamma)$ is an orbital rotation. The electronic gradient and Hessian with respect to the $\gamma_{ai}$ parameters is given by

[TABLE]

where the permutation operator is defined as $P_{ai,bj}A_{ai,bj}=A_{ai,bj}+A_{bj,ai}$. The $\kappa_{ai}$ parameters are determined by requiring the gradient to be zero and the Hessian to be positive definite, if possible. Once the Hartree-Fock state is determined, the excited states can be calculated using time-dependent Hartree-Fock (TDHF). The excitation energies $\omega$ are obtained from the non-Hermitian TDHF eigenvalue problem

[TABLE]

where

[TABLE]

$g_{pqrs}$ are the two-electron integrals, and $\varepsilon_{p}$ are the orbital energies. A comprehensive analysis of the stability of the TDHF eigenvalue problem has been given by Jørgensen and Simons [23]. In the Tamm-Dancoff approximation, the $\mathbf{B}$ matrix in Eq. 6 is discarded and the eigenvalue problem simply reads $\mathbf{A}\mathbf{x}=\omega\mathbf{x}$ [24].

When optimizing the Hartree-Fock energy, we may encounter convergence problems or multiple solutions to the ground state equations ($\mathbf{G}^{(0)}=\mathbf{0}$). To analyze this situation, we Taylor expand the gradient condition as

[TABLE]

Transforming to the eigenbasis of the Hessian, $\mathbf{G}^{(1)}\mathbf{r}_{n}=\lambda_{n}\mathbf{r}_{n}$, we obtain

[TABLE]

When approaching a conical intersection, the lowest positive eigenvalue $\omega_{1}$ in Eq. 6 becomes close to zero. And as zero-energy eigenvectors of Eq. 6 are also zero-energy eigenvectors of the Hessian [25], this implies that $\lambda_{1}$ will also approach zero. The first-order contribution is thus small, and the higher-order terms become important, giving rise to bifurcations in the energy. We define a modified gradient $\tilde{\mathbf{G}}^{(0)}$, where the projection along the first eigenvector of the Hessian $\mathbf{r}_{1}$ is removed, i.e.

[TABLE]

This new gradient allows us to solve along all the other eigenvectors and to implicitly define an effective Hessian $\tilde{\mathbf{G}}^{(1)}$ that is positive definite, forcing the energy function to be convex. We can express the $\kappa$ operator in the basis of the Hessian eigenvectors

[TABLE]

where $R_{n}=\sum_{ai}r_{ai,n}E^{-}_{ai}\,.$ Since we employ a modified gradient, we remove the corresponding component from $\bm{\kappa}$

[TABLE]

which guarantees $\tilde{\kappa}_{1}=0$. The final set of convex Hartree-Fock equations is given by

[TABLE]

where $\tilde{\mathbf{G}}^{(0)}$ and $\mathbf{G}^{(1)}$ are evaluated with $\tilde{\bm{\kappa}}$.

A detailed description of the employed algorithm can be found in the Method section. Since the equations are solved in a specific subspace, we avoid convergence issues related to degeneracies. The projected component can be introduced in a final diagonalization of the Hamiltonian matrix, written in the basis of $\{|\mathrm{HF}\rangle,|R_{1}\rangle,|\tilde{\nu}\rangle=|\nu\rangle-|R_{1}\rangle\langle R_{1}|\nu\rangle\}$. Here, $|R_{1}\rangle$ is defined as $\tfrac{1}{\sqrt{2}}\sum_{ai}r_{ai,1}E_{ai}|\mathrm{HF}\rangle$. The Hamiltonian expressed in this full space basis is

[TABLE]

where we used that $\langle\tilde{\mu}|H|\mathrm{HF}\rangle=\langle\mathrm{HF}|H|\tilde{\nu}\rangle=0$. The associated full space generalized eigenvalue problem is defined as

[TABLE]

The full space overlap matrix $\mathbf{S}^{\mathrm{FS}}$ is

[TABLE]

but it can be replaced with the identity matrix, as shown in the Supplementary Information. The above framework can easily be extended to an arbitrary number of states included in the projector operator. Details can be found in the Supplementary Information.

Applications

To illustrate the properties of the CVX-HF method, we consider the ammonia molecule. A conical intersection is encountered when one N-H bond is stretched and the molecule is in a nearly planar configuration [26, 13, 27].

In Fig. 1 (a) we show a region close to the intersection for TDA-TDHF, where the HF energy is higher than the excited state energy, resulting in a negative excitation energy shown by a blue region in the colormap at the bottom of the plot. When the geometry is planar ($\alpha=90^{\circ}$), the negative excitation energy extends to larger bond lengths, as can be seen from the blue segment in the colormap. Looking at the potential energy surfaces, it is clear that the surfaces in this region are not continuous with the surroundings. Using CVX-HF (Fig. 1 (b)) projecting the first eigenvector of the Hessian, we obtain continuous energy surfaces, and the intersection is limited to a single point. In the Supplementary Information, we provide a comparison of CVX-HF, TDA-TDHF, GCCSD [28], CCSD and FCI, when stretching one bond at a fixed $\alpha=89.5^{\circ}$.

It is instructive to analyze the molecular orbitals in the CVX-HF model. In Fig. 2 we consider the stretching of a single N-H bond at a fixed out-of-plane angle of $89.5^{\circ}$ and study the changes in energy and character of the $\sigma$ HOMO-1, $p_{N}$ HOMO and $\sigma^{*}$ LUMO orbitals. For TDA-TDHF we observe two state crossings: at $r=2.37$ Å and $r=2.65$ Å. Beyond this second point, a new HF solution (TDA-TDHF 2) appears with a negative excitation energy, corresponding to a state that is energetically lower than the original HF solution. In contrast, CVX-HF shows a single avoided crossing at $r=2.37$Å. The energy profile shows a smooth and continuous connection from the initial TDA-TDHF solution to the TDA-TDHF 2 after the avoided crossing. The character of the CVX-HF orbitals is consistent across the entire range of bond lengths. This is not the case for TDA-TDHF, as the $p_{N}$ orbital gets mixed after the first crossing. The original picture is retrieved with TDA-TDHF 2, only after the second crossing point.

We now consider a larger system and study an $S_{0}$/$S_{1}$ intersection in 2,4-cyclohexadien-1-ylamine. The $\mathbf{g}$ and $\mathbf{h}$ vectors are obtained with a CCSD algorithm [29] at the geometry in Ref. 30. Additional information can be found in the Supplementary Information and in Ref. 28. The potential energy surfaces for TDA-TDHF are shown in Fig. 3 (a). A region in which HF did not converge is present for $g$ greater than 2.25 and $h=2.75$. The two surfaces are discontinuous across this region, resulting in an overall non-conical shape. In contrast, the results for CVX-HF in Fig. 3 (b) show no convergence issues and the surfaces are continuous, retaining the expected conical shape.

Lastly, we examine the $S_{0}$/$S_{1}$ conical intersection in the anionic form of p-hydroxybenzylidene-2,3-dimethylimidazolinone (HBDI-) [31]. This molecule is the chromophore of the green fluorescent protein (GFP), a widely used fluorescent marker in biological imaging. Its bright green emission arises from a photoinduced process closely linked to the electronic structure of HBDI- [32, 33]. The torsional angle $\phi_{P}$ is defined as the dihedral angle between the carbon atoms of the methine bridge and C1, and $\phi_{P}$ is defined as the dihedral angle between N1 and the carbon atoms of the methine bridge (see Fig. 4). The initial structure corresponds to the P90 minimum energy conical intersection (MECI) in Ref. 31, but with modified dihedral angles as described in the Supplementary Information. In Fig. 4 (a), a region is observed where HF fails to converge. In addition, the behavior of the potential energy surfaces for $\phi_{P}>72.5^{\circ}$ is significantly different from that at smaller $\phi_{P}$ values. In the latter case, the energy shows little variation as $\phi_{I}$ changes. Again, the CVX-HF in Fig. 4 (b) retains the correct topology of the intersection.

Discussion

In this work, we introduced the Convex Hartree-Fock framework as an efficient method for modeling ground state conical intersections. By reformulating the Hartree-Fock optimization in a subspace where it is convex, the approach yields smooth potential energy surfaces and resolves failures inherent to conventional mean-field methods. The CVX-HF framework has several important fundamental properties. The method is orbital invariant with respect to rotations within either the occupied or the virtual orbital spaces. Furthermore, as shown in the Supplementary Information, the excitation energies are size-intensive when additional molecules are introduced, such as in solvated systems. The calculated excitation energies are rather insensitive to the number of states that are being projected (see Tables S2 and S3). Therefore, the choice of the number of projected states is not critical and can be adjusted to the number of states involved, for instance, in the dynamics studied.

Despite its advantages, CVX-HF remains fundamentally at the Hartree–Fock level of accuracy, which can limit the reliability for quantitative predictions of excitation energies and electron correlation effects. To overcome this limitation, two strategies can be employed. First, the orbitals determined in the CVX-HF framework can be used as a basis for more accurate methods, such as coupled cluster models and perturbation theory, ensuring that discontinuities at the HF level do not propagate to the correlated model. This is particularly important close to ground state conical intersections, where coupled cluster displays divergent behavior due to the inability to describe the geometric phase effect [9]. The CVX-HF orbitals can also be combined with the generalized coupled cluster framework that was developed to give a correct description of ground state conical intersections [28]. Second, the framework presented in this work can be applied to TDDFT, using the same exponential parameterization to define the orbitals in the Kohn-Sham equations. The component along the first Hessian eigenvector, determined during the ground state optimization, is removed from the gradient and the $\kappa$ matrix, and later introduced in the final diagonalization step.

A promising feature of the method is that additional excitations can be included in the final diagonalization procedure. A notable example involves systems with multiple chromophores, as in the case of pigments in photosynthetic light-harvesting complexes [34, 35]. When the molecules are infinitely separated, the eigenvectors of the Hessian of the combined system are intensive and will be localized on the individual molecules. These eigenvectors may be used to construct a size-extensive parameterization. As an example, consider the case of two non-interacting systems, $A$ and $B$. The wave function of the combined system should match the product of the wave functions of the isolated systems. For CVX-HF, the wave functions of the two isolated systems in the reduced space can be expressed as

[TABLE]

The wave function for the combined system in the reduced space, when projecting 2 states, is given by

[TABLE]

where

[TABLE]

This wave function is not size-extensive because the component given by the combined excitation $|R_{1}R_{2}\rangle=|R_{A}\rangle|R_{B}\rangle$ is missing. In order to retrieve the correct product form, this term can be included in the final reduced space eigenvalue problem, only changing the prefactor of the overall CVX-HF scaling. The size-extensivity of the combined system in the extended space is shown in the Supplementary Information. The same idea can be generalized to a larger number of chromophores. Finally, we note that, in some cases, double excitations are needed in order to describe the conical intersection. These can be added in a similar manner.

Methods

All calculations reported in this work were performed using an implementation of the CVX-HF framework in a local development version of the eT program [36]. The algorithm presented below outlines the iterative procedure employed in CVX-HF to self-consistently solve the Hessian eigenvalue problem and optimize the ground state within the restricted subspace. The full space eigenvalue problem is solved once at the end of the procedure. Additional details on the definitions of vectors and matrices are provided in the Supplementary Information.

The eigenvalue problem in step 4 is solved using Davidson’s algorithm, avoiding the explicit calculation of second-derivative terms. In each iteration, the threshold is updated as $\text{min}(10^{-2},\lVert\tilde{\mathbf{G}}^{(0)[n]}\rVert_{L_{2}})$ to match the convergence level of the ground state equations. This approach ensures that the computational cost of this step is distributed across the entire cycle, instead of unnecessarily solving the problem to full accuracy at every iteration. Step 7 is solved using a trust-region algorithm, again avoiding the explicit calculation of second-derivative terms. The final eigenvalue problem in step 13 is also solved using Davidson’s algorithm, initialized with the eigenvectors obtained from the direct diagonalization in step 12. When multiple calculations are performed for the same system at similar geometries, the final value of $\bm{\kappa}$ from a previous calculation can be used as the initial guess for the next calculation. This requires consistency among the SAD guesses of both calculations, which is ensured through a diabatization step of the new $C_{0}$ coefficients with respect to those of the previous calculation.

Acknowledgments

We thank Eirik F. Kjønstad for enlightening discussions. This work was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Program (grant agreement No. 101020016).

Author contributions

F.R. and H.K. conceived the CVX-HF framework, analyzed the data and wrote the paper. F.R. developed the implementation in eT and performed all the calculations. H.K. supervised the project.

Competing interests

The authors declare no competing interests.