Convex Hartree–Fock theory for modeling ground state conical intersections

TL;DR

This paper introduces a modified Hartree–Fock method to better model ground state conical intersections in nonadiabatic molecular dynamics.

Contribution

The novel Convex Hartree–Fock framework optimizes reference states in a tailored subspace to capture ground state conical intersections efficiently.

Findings

Convex Hartree–Fock successfully captures ground state conical intersections where traditional methods fail.

The method is validated across multiple test cases and compared to time-dependent Hartree–Fock benchmarks.

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.

Accurate modeling of conical intersections is essential for describing the nonadiabatic molecular dynamics behind photoinduced processes, but conventional single-reference electronic structure methods fail to capture the associated ground state intersections. In this work, the authors 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, offering an alternative to more expensive multiconfigurational approaches.

Figures

Potential energy surfaces of S_0_ and S_1_ in ammonia.Comparison of TDA-TDHF (a) and CVX-HF (b) potential energy surfaces of S_0_ and S_1_ in ammonia. The basis set is aug-cc-pVDZ. For each point the energies are plotted in eV relative to the average energy $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{2}$$\end{document}$ (*E*_0_ + *E*_1_) of the states. In (d), the CVX-HF total energy potential energy

Analysis of the orbitals and potential energy curves of states S_0_ and S_1_ in ammonia.The basis set is 6-31G*. One N–H bond is stretched while maintaining a fixed out-of-plane angle of *α* = 89. 5^∘^ (see Fig. [1](#Fig1) for details). All energies are reported in Hartrees. CVX-HF is shown as black and gray crosses, TDA-TDHF as blue and red lines, and beyond *r* = 2.65 Å, a second Hartree–Fock solution TDA-TDHF 2 is shown with green and orange lines. The HOMO-1, HOMO, and LUMO orbitals are shown for all methods at three representative bond lengths: before (*r* = 1.385 Å), near (*r* = 2.416 Å)

Potential energy surfaces of S_0_ and S_1_ in 2,4-cyclohexadien-1-ylamine.Comparison of TDA-TDHF (a) and CVX-HF (b) potential energy surfaces of S_0_ and S_1_ in 2,4-cyclohexadien-1-ylamine. The basis set is cc-pVDZ. For each point, the energies are plotted in eV relative to the average energy $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{2}$$\end{document}$ (*E*_0_ + *E*_1_) of the states. In (d), the

Potential energy surfaces of S_0_ and S_1_ in HBDI^−^.Comparison of TDA-TDHF (a) and CVX-HF (b) potential energy surfaces of S_0_ and S_1_ in HBDI^−^.The basis set is 6-31G*. For each point the energies are plotted in eV relative to the average energy $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{2}$$\end{document}$ (*E*_0_ + *E*_1_) of the states. In (d), the CVX-HF total energy potential energy surfac

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^. In this discussion, we limited ourselves to the case of linear-response TDDFT in the adiabatic approximation, which is the standard formulation in practical applications. Beyond this approximation, non-linear terms can explicitly introduce a coupling between ground and excited states^4,5^.

The stability of Hartree–Fock solutions has been studied in great detail previously^6–8^. 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^9,10^. The non-analytical behavior has been discussed by Čížek and Paldus^8^ and was recently shown to be particularly pronounced in areas around ground state conical intersections^11^. 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^12^, which redefines the reference state to access crossings between ground and excited states by allowing spin-changing excitations. Other approaches include dual-functional TDA (DF-TDA)^13^, and TDDFT-1D^14^, where specific excited configurations are incorporated into the eigenvalue problem. Configuration interaction-corrected DFT (CIC-TDA^15^) 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^16^), multiconfigurational short-range DFT^17^, and ensemble-based DFT methods adapted to describe excited states with near-degeneracies (SI-SA-REKS^18^). In addition, mode-following techniques^19^ have been employed to directly locate excited state saddle points by following specific eigenvectors of the Hessian. Most recently, phase-space electronic structure theory^20^ 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^21^.

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^22,23^, 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert {\Phi }_{0}\right\rangle$$\end{document} constructed from orbitals obtained from the superposition of atomic densities (SAD)^24,25^, as this ensures a well-behaving starting point. The Hartree–Fock wave function is parameterized in terms of a single global orbital rotation, such that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\vert {\rm{HF}}\rangle =\exp (\sum _{ai}{\kappa }_{ai}{E}_{ai}^{-})\vert {\Phi }_{0}\rangle ,$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E}_{ai}^{-}={E}_{ai}-{E}_{ia}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E}_{ai}=\sum _{\sigma }{a}_{a\sigma }^{\dagger }{a}_{i\sigma }$$\end{document} , and κa**i are the orbital rotation parameters. We use a, b to label virtual orbitals and i, j for occupied orbitals. The Hartree–Fock energy \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}_{{\rm{HF}}}$$\end{document} is obtained as the expectation value of the electronic Hamiltonian H,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{E}}}_{{\rm{HF}}}=\langle {\rm{HF}}\vert H\vert {\rm{HF}}\rangle ,$$\end{document}

where we assume \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert {\rm{HF}}\right\rangle$$\end{document} is normalized. When minimizing the energy with respect to the κa**i parameters, we need to calculate the electronic gradient and Hessian. These are conveniently obtained by considering the Hartree–Fock energy function

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{E}}({\boldsymbol{\gamma }})=\langle {\rm{HF}}\vert \exp (-\gamma )H\exp (\gamma )\vert {\rm{HF}}\rangle ,$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\exp (\gamma )$$\end{document} is an orbital rotation. The electronic gradient and Hessian with respect to the γa**i parameters is given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${G}_{ai}^{(0)}({\bf{0}})={\frac{\partial {\mathcal{E}}}{\partial {\gamma }_{ai}}\bigg| }_{{\boldsymbol{\gamma }} = {\bf{0}}}= \langle {\rm{HF}}\vert [H,{E}_{ai}^{-}]\vert {\rm{HF}}\rangle,$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${G}_{ai,bj}^{(1)}({\bf{0}})={\frac{{\partial }^{2}{\mathcal{E}}}{\partial {\gamma }_{ai}\partial {\gamma }_{bj}}\bigg| }_{{\boldsymbol{\gamma }} = {\bf{0}}} = \frac{1}{2}{P}_{ai,bj}\langle {\rm{HF}}\vert [[H,{E}_{ai}^{-}],{E}_{bj}^{-}]\vert {\rm{HF}}\rangle ,$$\end{document}

where the permutation operator is defined as Pa**i,*b**j *Aa**i,b**j = Aa**i,b**j + Ab**j,a**i. The κa**i 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 ω are obtained from the non-Hermitian TDHF eigenvalue problem

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left(\begin{array}{ll}{\bf{A}}&{\bf{B}}\\ {\bf{B}}&{\bf{A}}\end{array}\right)\left(\begin{array}{l}{\bf{X}}\\ {\bf{Y}}\end{array}\right)=\omega \left(\begin{array}{ll}{\bf{1}}&{\bf{0}}\\ {\bf{0}}&-{\bf{1}}\end{array}\right)\left(\begin{array}{l}{\bf{X}}\\ {\bf{Y}}\end{array}\right),$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${A}_{ia,jb}={\delta }_{ij}{\delta }_{ab}({\varepsilon }_{a}-{\varepsilon }_{i})+2{g}_{aijb}-{g}_{abji},$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${B}_{ia,jb}=2{g}_{aibj}-{g}_{ajbi},$$\end{document}

gpqr**s are the two-electron integrals, and εp are the orbital energies. A comprehensive analysis of the stability of the TDHF eigenvalue problem has been given by Jørgensen and Simons^26^. In the Tamm-Dancoff approximation, the B matrix in Eq. (6) is discarded and the eigenvalue problem simply reads A****x = ωx^27^. This equation is equivalent to the formulation of CIS^28^, but we will refer to them as TDA-TDHF to maintain a clearer parallel with the LR-TDDFT framework.

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

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${G}_{ai}^{(0)}({\boldsymbol{\gamma }})={G}_{ai}^{(0)}({\bf{0}})+\sum _{bj}{G}_{ai,bj}^{(1)}({\bf{0}}){\gamma }_{bj}+\ldots =0.$$\end{document}

Transforming to the eigenbasis of the Hessian, G^(1)^rn = λnrn, we obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${G}_{n}^{(0)}({\boldsymbol{\gamma }})={G}_{n}^{(0)}({\bf{0}})+{\lambda }_{n}{\gamma }_{n}+\ldots =0.$$\end{document}

When approaching a conical intersection, the lowest positive eigenvalue ω1 in Eq. (6) becomes close to zero. And as zero-energy eigenvectors of Eq. (6) are also zero-energy eigenvectors of the Hessian^29^, this implies that λ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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{{\bf{G}}}}^{(0)}$$\end{document} , where the projection along the first eigenvector of the Hessian r1 is removed, i.e.,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{{\bf{G}}}}^{(0)}={{\bf{G}}}^{(0)}-{{\bf{r}}}_{1}({{\bf{r}}}_{1}^{T}{{\bf{G}}}^{(0)}).$$\end{document}

This new gradient allows us to solve along all the other eigenvectors and to implicitly define an effective Hessian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{{\bf{G}}}}^{(1)}$$\end{document} that is positive definite, forcing the energy function to be convex. We can express the κ operator in the basis of the Hessian eigenvectors

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa =\sum _{ai}{\kappa }_{ai}{E}_{ai}^{-}=\sum _{n}{\kappa }_{n}{R}_{n}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${R}_{n}=\sum _{ai}{r}_{ai,n}{E}_{ai}^{-}.$$\end{document} Since we employ a modified gradient, we remove the corresponding component from κ

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\kappa }}_{ai}={\kappa }_{ai}-{r}_{ai,1}\sum _{bj}{r}_{bj,1}{\kappa }_{bj},$$\end{document}

which guarantees \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{\kappa }}_{1}=0$$\end{document} . The final set of Convex Hartree–Fock equations is given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{{\bf{G}}}}^{(0)}=0 ,$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{G}}}^{(1)}{{\bf{r}}}_{1}={\lambda }_{1}{{\bf{r}}}_{1},$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{{\bf{G}}}}^{(0)}$$\end{document} and G^(1)^ are evaluated with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{{\boldsymbol{\kappa }}}$$\end{document} .

A detailed description of the employed algorithm can be found in the Methods 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{\left\vert {\rm{HF}}\right\rangle ,\left\vert {R}_{1}\right\rangle ,\left\vert \tilde{\nu }\right\rangle =\left\vert \nu \right\rangle -\left\vert {R}_{1}\right\rangle \left\langle {R}_{1}| \nu \right\rangle \}$$\end{document} . Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert {R}_{1}\right\rangle$$\end{document} is defined as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{\sqrt{2}}\sum _{ai}{r}_{ai,1}{E}_{ai}\left\vert {\rm{HF}}\right\rangle$$\end{document} . The Hamiltonian expressed in this full space basis is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{H}}}^{{\rm{FS}}}=\left(\begin{array}{lll}\left\langle {\rm{HF}}\right\vert H\left\vert {\rm{HF}}\right\rangle &\left\langle {\rm{HF}}\right\vert H\left\vert {R}_{1}\right\rangle &0\\ \left\langle {R}_{1}\right\vert H\left\vert {\rm{HF}}\right\rangle \hfill&\left\langle {{R}}_{1}\right\vert H\left\vert {R}_{1}\right\rangle &\left\langle {R}_{1}\right\vert H\left\vert \tilde{\nu }\right\rangle \\ 0 &\left\langle \tilde{\mu }\right\vert H\left\vert {R}_{1}\right\rangle &\left\langle \tilde{\mu }\right\vert H\left\vert \tilde{\nu }\right\rangle \hfill\end{array}\right),$$\end{document}

where we used that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\langle \tilde{\mu }\right\vert H\left\vert {\rm{HF}}\right\rangle =\left\langle {\rm{HF}}\right\vert H\left\vert \tilde{\nu }\right\rangle =0$$\end{document} . The associated full space generalized eigenvalue problem is defined as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{H}}}^{{\rm{FS}}}{{\bf{x}}}_{n}={{\mathcal{E}}}_{n}{{\bf{S}}}^{{\rm{FS}}}{{\bf{x}}}_{n},$$\end{document}

which provides ground and excited state energies. The full space overlap matrix S^FS^ is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{S}}}^{{\rm{FS}}}=\left(\begin{array}{lll}1 &\left\langle {\rm{HF}}| {R}_{1}\right\rangle &\left\langle {\rm{HF}}| \tilde{\nu }\right\rangle \\ \left\langle {R}_{1}| {\rm{HF}}\right\rangle \hfill&1 &\left\langle {R}_{1}| \tilde{\nu }\right\rangle \\ \left\langle \tilde{\mu }| {\rm{HF}}\right\rangle \hfill&\left\langle \tilde{\mu }| {R}_{1}\right\rangle \hfill&\left\langle \tilde{\mu }| \tilde{\nu }\right\rangle \end{array}\right)=\left(\begin{array}{lll}1&0&0\\ 0&1&0\\ 0&0&\left\langle \tilde{\mu }| \tilde{\nu }\right\rangle \end{array}\right) ,$$\end{document}

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

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^15,30,31^.

In Fig. 1a, 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 (α = 90^∘^), 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. 1b) projecting the first eigenvector of the Hessian, we obtain continuous energy surfaces, and the intersection is limited to a single point. In Figs. S1–S2 of the Supplementary Information we provide a comparison of CVX-HF, TDA-TDHF, GCCSD^32^, CCSD and FCI, when stretching one bond at a fixed α = 89.5^∘^.Fig. 1. Potential energy surfaces of S_0_ and S_1_ in ammonia.Comparison of TDA-TDHF (a) and CVX-HF (b) potential energy surfaces of S_0_ and S_1_ in ammonia. The basis set is aug-cc-pVDZ. For each point the energies are plotted in eV relative to the average energy \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{2}$$\end{document} (E0 + E1) of the states. In (d), the CVX-HF total energy potential energy surfaces for the same system are shown, expressed in Hartree. At the bottom of each plot, a colormap shows the value of E1 − E0. The geometry is shown in (c), with the bond length r1 and the out-of-plane angle α indicated in red. One N–H bond is stretched while the remaining two are fixed at 1.04 Å. The out-of-plane angle α is defined as the angle between the vector trisecting the three N–H bonds and any of the N–H bonds. For α = 90^∘^ the geometry is planar, with all H–N–H angles equal to 120^∘^. These angles are kept equal as α changes.

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^∘^ and study the changes in energy and character of the σ HOMO-1, pN HOMO and σ^^ 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.Fig. 2. Analysis of the orbitals and potential energy curves of states S_0_ and S_1_ in ammonia.The basis set is 6-31G. One N–H bond is stretched while maintaining a fixed out-of-plane angle of α = 89. 5^∘^ (see Fig. 1 for details). All energies are reported in Hartrees. CVX-HF is shown as black and gray crosses, TDA-TDHF as blue and red lines, and beyond r = 2.65 Å, a second Hartree–Fock solution TDA-TDHF 2 is shown with green and orange lines. The HOMO-1, HOMO, and LUMO orbitals are shown for all methods at three representative bond lengths: before (r = 1.385 Å), near (r = 2.416 Å), and after (r = 3.133 Å) the avoided crossing. At r = 1.385 Å, only CVX-HF orbitals are shown, as those from TDA-TDHF are visually indistinguishable. Similarly, at r = 3.133 Å, only the CVX-HF orbitals are shown, since they closely match those from TDA-TDHF 2.

We now consider a larger system and study an S0/S1 intersection in 2,4-cyclohexadien-1-ylamine. The g and h vectors are obtained with a CCSD algorithm^33^ at the geometry used in previous works^32,34^ and reported in Table S9 of the Supplementary Information. The potential energy surfaces for TDA-TDHF are shown in Fig. 3a. 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. 3b show no convergence issues and the surfaces are continuous, retaining the expected conical shape.Fig. 3. Potential energy surfaces of S_0_ and S_1_ in 2,4-cyclohexadien-1-ylamine.Comparison of TDA-TDHF (a) and CVX-HF (b) potential energy surfaces of S_0_ and S_1_ in 2,4-cyclohexadien-1-ylamine. The basis set is cc-pVDZ. For each point, the energies are plotted in eV relative to the average energy \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{2}$$\end{document} (E0 + E1) of the states. In (d), the CVX-HF total energy potential energy surfaces for the same system are shown, expressed in Hartree. At the bottom of each plot, a colormap shows the value of E1 − E0 in eV. The geometry at the CVX-HF CI is shown in (c).

Lastly, we examine the S0/S1 conical intersection in the anionic form of p-hydroxybenzylidene-2,3-dimethylimidazolinone (HBDI^−^)^35^. 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^−^^36,37^. The torsional angle ϕP is defined as the dihedral angle between the carbon atoms of the methine bridge and C_1_, and ϕI is defined as the dihedral angle between N_1_ and the carbon atoms of the methine bridge (see Fig. 4). The initial structure corresponds to P90 minimum energy conical intersection (MECI)^35^, but with modified dihedral angles as described in Table S14 of the Supplementary Information. In Fig. 4a, a region is observed where HF fails to converge. In addition, the behavior of the potential energy surfaces for ϕP > 72. 5^∘^ is significantly different from that at smaller ϕP values. In the latter case, the energy shows little variation as ϕI changes. Again, the CVX-HF in Fig. 4b retains the correct topology of the intersection.Fig. 4. Potential energy surfaces of S_0_ and S_1_ in HBDI^−^.Comparison of TDA-TDHF (a) and CVX-HF (b) potential energy surfaces of S_0_ and S_1_ in HBDI^−^.The basis set is 6-31G*. For each point the energies are plotted in eV relative to the average energy \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{2}$$\end{document} (E0 + E1) of the states. In (d), the CVX-HF total energy potential energy surfaces for the same system are shown, expressed in Hartree. At the bottom of each plot, a colormap shows the value of E1 − E0 in eV. The geometry at the CVX-HF CI is shown in (c).

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 Supplementary Notes 4–5, 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 projected (see Tables S2–S5). 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^11^. 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^32^. 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 κ 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^38,39^. 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

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert {\psi }_{A}\right\rangle =\left\vert {{\rm{HF}}}_{A}\right\rangle {x}_{0}^{A}+\left\vert {R}_{A}\right\rangle {x}_{1}^{A} ,$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert {\psi }_{B}\right\rangle =\left\vert {{\rm{HF}}}_{B}\right\rangle {x}_{0}^{B}+\left\vert {R}_{B}\right\rangle {x}_{1}^{B}.$$\end{document}

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

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert \psi \right\rangle =\left\vert {\rm{HF}}\right\rangle {x}_{0}+\left\vert {R}_{1}\right\rangle {x}_{1}+\left\vert {R}_{2}\right\rangle {x}_{2}\,,$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert {\rm{HF}}\right\rangle =\left\vert {{\rm{HF}}}_{A}\right\rangle \left\vert {{\rm{HF}}}_{B}\right\rangle ,$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert {R}_{1}\right\rangle =\left\vert {R}_{A}\right\rangle \left\vert {{\rm{HF}}}_{B}\right\rangle ,$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert {R}_{2}\right\rangle =\left\vert {{\rm{HF}}}_{A}\right\rangle \left\vert {R}_{B}\right\rangle .$$\end{document}

This wave function is not size-extensive because the component given by the combined excitation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\vert {R}_{1}{R}_{2}\right\rangle =\left\vert {R}_{A}\right\rangle \left\vert {R}_{B}\right\rangle$$\end{document} 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 Supplementary Note 3. 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^40,41^. 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 Supplementary Note 1.

Algorithm 1

CVX-HF algorithm

1: C0 ← SAD guess

2: Initialize n = 0, κ^[0]^ = 0 and C^[0]^ = C0

3: while ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n < {n}_{\max }$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\parallel {\tilde{{\bf{G}}}}^{(0)}{\parallel }_{{L}_{2}} > $$\end{document} threshold) do

4:  Solve the eigenvalue problem \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{G}}}^{(1)[n]}{{\bf{r}}}_{i}^{[n]}={\lambda }_{i}^{[n]}{{\bf{r}}}_{i}^{[n]}$$\end{document} , i = 1, …, nexcited

5:  Construct G^(0)[n]^ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{P}}^{[n]}={\mathcal{I}}-\mathop{\sum }\limits_{i=1}^{{n}_{{\rm{proj}}}}{{\bf{r}}}_{i}^{[n]}{{\bf{r}}}_{i}^{[n]T}$$\end{document}

6:   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tilde{{\bf{G}}}}^{(0)[n]}\leftarrow {\hat{P}}^{[n]}{{\bf{G}}}^{(0)[n]}$$\end{document}

7:  Solve \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hat{P}}^{[n]}{{\bf{G}}}^{(1)[n]}{\hat{P}}^{[n]}{\rm{\Delta }}{{\boldsymbol{\kappa }}}^{[n]}=-{\tilde{{\bf{G}}}}^{(0)[n]}$$\end{document} for Δκ^[n]^

8:   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{\kappa }}}^{[n+1]}\leftarrow {\hat{P}}^{[n]}({{\boldsymbol{\kappa }}}^{[n]}+{\rm{\Delta }}{{\boldsymbol{\kappa }}}^{[n]})$$\end{document}

9:   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C}^{[n+1]}\leftarrow {C}_{0}\,\,\exp ({{\boldsymbol{\kappa }}}^{[n+1]})$$\end{document}

10:  n ← n + 1

11: end while

12: Construct and diagonalize the reduced matrix H^RS^

13: Solve the final eigenvalue problem \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{H}}}^{{\rm{FS}}}{{\bf{x}}}_{j}={{\mathcal{E}}}_{j}{{\bf{x}}}_{j},$$\end{document} j = 1, …, nexcited + 1

The first set of orbital coefficients C0 is obtained diagonalizing the Fock matrix built using the SAD density. We point out that the quality of the initial orbital coefficients is important as the optimization is restricted to a subspace, and along one or more directions the initial guess is not explicitly optimized. 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\min (1{0}^{-2},\parallel {\tilde{{\bf{G}}}}^{(0)[n]}{\parallel }_{{L}_{2}})$$\end{document} 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 κ 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 C0 coefficients with respect to those of the previous calculation. In this way, any reordering of the orbitals is taken into account and all the elements of κ corresponds to the correct orbital pair. A comparison of timings and number of iterations between CVX-HF and conventional Hartree–Fock can be found in Tables S6-S7 of the Supplementary Information.

Supplementary information

Supplementary Information