# Calculations of Excited Electronic States by Converging on Saddle Points
Using Generalized Mode Following

Yorick L. A. Schmerwitz  Science Institute and Faculty of Physical Sciences,
University of Iceland, Reykjavík, Iceland  Gianluca Levi<EMAIL_ADDRESS>Science
Institute and Faculty of Physical Sciences, University of Iceland, Reykjavík,
Iceland  Hannes Jónsson<EMAIL_ADDRESS>Science Institute and Faculty of Physical
Sciences, University of Iceland, Reykjavík, Iceland

###### Abstract

Variational calculations of excited electronic states are carried out by
finding saddle points on the surface that describes how the energy of the
system varies as a function of the electronic degrees of freedom. This
approach has several advantages over commonly used methods especially in the
context of density functional calculations, as collapse to the ground state is
avoided and yet, the orbitals are variationally optimized for the excited
state. This optimization makes it possible to describe excitations with large
charge transfer where calculations based on ground state orbitals are
problematic, as in linear response time-dependent density functional theory. A
generalized mode following method is presented where an $n$th-order saddle
point is found by inverting the components of the gradient in the direction of
the eigenvectors of the $n$ lowest eigenvalues of the electronic Hessian
matrix. This approach has the distinct advantage of following a chosen excited
state through atomic configurations where the symmetry of the single
determinant wave function is broken, as demonstrated in calculations of
potential energy curves for nuclear motion in the ethylene and dihydrogen
molecules. The method is implemented using a generalized Davidson algorithm
and an exponential transformation for updating the orbitals within a
generalized gradient approximation of the energy functional. Convergence is
found to be more robust than for a direct optimization approach previously
shown to outperform standard self-consistent field approaches, as illustrated
here for charge transfer excitations in nitrobenzene and N-phenylpyrrole,
involving calculations of 4th\- and 6th-order saddle points, respectively.
Finally, calculations of a diplatinum and silver complex are presented,
illustrating the applicability of the method to excited state energy curves of
large molecules.

##  1 Introduction

Calculations of photochemical processes require electronic structure methods
that can describe a wide range of excitations and give the variation of the
energy of the excited states as a function of the atomic coordinates. Linear-
response time-dependent density functional theory1, 2 in the adiabatic
approximation (henceforth referred to as TDDFT) is a commonly used method for
the calculation of excited electronic states, due to its relatively low
computational cost. However, TDDFT has some important limitations3. Firstly,
due to the linear-response formalism, TDDFT with local and semi-local
functionals cannot adequately describe excitations involving large changes in
the electron density, such as charge transfer, core, and Rydberg
excitations.4, 5, 6. Secondly, as the atomic configuration changes, the
reference state of a TDDFT calculation can abruptly change character, leading
to a discontinuity in the energy of the excited state, thereby making the
evaluation of atomic forces ill-defined. This behavior is, e.g., observed in
bond stretching7, 8 and in the vicinity of electronic degeneracy,9 such as a
conical intersection. The latter can play a major role in photochemistry, as
the Born-Oppenheimer approximation breaks down, and population transfer can
readily occur between the electronic states. The topology of conical
intersections is described incorrectly in TDDFT10 because doubly excited
configurations are missing in the response state, and the degeneracy is
thereby not lifted along one of the branching space atomic displacements.
Several modifications of TDDFT have been proposed to improve the description
of excited states, such as a configuration interaction using one doubly
excited configuration,11, 12 the dual-functional approach,13, 14 the spin-flip
TDDFT15, and the particle-particle16 and hole-hole17 Tamm-Dancoff
approximations. The description of charge transfer states can be improved by
using higher-level functional approximations, such as double hybrid18, 19 and
optimally tuned20, 21, 22 functionals, but these enhancements involve a
significant increase in computational effort. Moreover, optimal tuning is
system- and geometry-specific and is, therefore, of limited applicability in
simulations of atomic dynamics in excited states.

Variational density functional calculations of excited states23, 6, 24, 25, 26
are emerging as an attractive alternative to TDDFT. They typically involve
similar computational cost as ground state calculations and can better
describe long-range charge transfer6, 9, 27, Rydberg28, 29, core-level30, 31,
and other excitations6, 32 where a significant change in the electron density
occurs. As the calculations are variational, they provide atomic forces that
can be used in excited state geometry optimization and classical dynamics
simulations33, 34, 35, 36. In a variational calculation within a mean-field
approximation, an excited state is found as a solution of higher energy than
the ground state and corresponds to an optimal single Slater determinant with
non-aufbau orbital occupation. An important feature of an excited state
stationary solution is that it typically corresponds to a saddle point on the
surface that describes the variation of the energy as a function of the
electronic degrees of freedom, while the ground state corresponds to a
minimum.

On the exact electronic energy landscape, provided by full configuration
interaction, the $n$th excited state above the ground state corresponds to a
saddle point of order $n$37, 38. Accordingly, state-specific
multiconfigurational self-consistent field (MCSCF) calculations of excited
states have employed optimization strategies that specifically search for
saddle points39, 40, 41 on the MCSCF energy surface. For mean-field
approximations, excited states do not always correspond to saddle points, but
it is usually observed that the number of negative eigenvalues of the
electronic Hessian at a stationary solution increases with the energy
(although not necessarily monotonically 37), so the mean-field approximate
excited states are typically saddle points6, 24, 32, 42. Commonly used methods
for variational mean-field calculations of excited states have, however, not
been specifically designed to find saddle points on the electronic energy
surface. The most common approach is the SCF method based on Hamiltonian
matrix diagonalization. Since the excitation energy is computed as the
difference relative to the ground state energy, the term “$\Delta$SCF” is
often used in the literature for such calculations 3, 23, 33, 34, 43, 35, 36.
SCF calculations of excited states typically use iterative eigensolvers
developed for ground state calculations, such as the Davidson44 method or the
direct inversion in the iterative subspace (DIIS)45, 46. Additionally, a
maximum overlap method47, 48, 9 (MOM) is often used to reduce the likelihood
of converging on the ground state, i.e. variational collapse (we will refer to
this approach of using an SCF procedure with MOM as SCF-MOM). Although methods
based on SCF can converge on saddle points, they are better suited for the
minimization needed for ground state calculations. Moreover, they can have
convergence problems close to electronic near-degeneracies, e.g. in the
vicinity of conical intersections or points at which symmetry-broken solutions
arise, as are often encountered in bond breaking processes.

Recently, a direct optimization (DO) method has been developed49, 50, 24 based
on an exponential transformation and a quasi-Newton algorithm for finding the
pairwise orbital rotations that make the energy stationary. This approach can
converge on saddle points on the electronic energy surface when a
preconditioned quasi-Newton method is used together with MOM (referred to as
DO-MOM method) to maximize the energy along some electronic degree(s) of
freedom corresponding to directions of negative curvature. Combined with the
robustness typical of approximate second-order optimization in the vicinity of
an electronic degeneracy24, 51, 52, these attributes make DO-MOM a useful tool
for variational excited state calculations. For example, by using DO-MOM, it
has been shown53 that symmetry-broken ground and excited state solutions
provide a qualitatively correct conical intersection and avoided crossing in
the ethylene molecule, even when a semi-local generalized gradient
approximation functional is used.54, 55 The application of explicit self-
interaction correction improves the calculated results and gives nearly
quantitative agreement with high-level quantum chemistry calculations.53 A
previous study using an SCF-MOM approach turned out to be less successful.35

Despite these advancements in the algorithms for variational density
functional calculations of excited states, some limitations remain which
hamper widespread application of the methodology. One issue is related to the
fact that multiple solutions providing qualitatively different variation of
the energy as a function of the atomic coordinates (the molecular potential
energy surface) can exist for a given excited state. For example, when
photoexcitation involves bond breaking, different electronic configurations
approach degeneracy, leading to increased static correlation and the emergence
of multiple mean-field solutions, some conserving and others breaking the
symmetries of the Hamiltonian.56, 57, 58, 59, 60 We will refer to atomic
configurations where multiple mean-field solutions start to appear as
symmetry-breaking onsets (SBOs). Typically, symmetry-broken solutions provide
an adequate description of the energy surface beyond SBOs, due to their
implicit multiconfigurational character53, 61, 62, 63, 64, 65, 66, 67, 68.
Other solutions, however, give qualitatively incorrect potential energy
surfaces for atom dynamics, for instance corresponding to diabatic surfaces53,
which can lead to unphysical state switching in geometry optimizations and
classical dynamics simulations. As the ultimate goal of variational excited
state calculations is often the simulation of atomic dynamics in the excited
state 69, 33, 34, 35, it is important to ensure that the calculations converge
on the solution providing the appropriate value of the energy and atomic
forces. While the ground state solution can, by the minimum energy principle,
be identified as the one with the lowest energy,70 an analogous method cannot
be used for tracking excited state solutions.53, 38 In previous work53, we
demonstrated that the mean-field solution for the doubly excited state of
ethylene giving the best description of the energy curve along the double bond
torsion can be identified based on the conservation of the saddle point order
as the wave function undergoes symmetry breaking. This approach provides a
guiding principle for targeting a specific solution during the exploration of
an excited state potential energy surface for atomic dynamics. A second
remaining issue is problematic convergence on a high-order saddle point,
especially when wave function relaxation leads to large rearrangements of the
orbitals, as can be the case for charge transfer excitations. Although the DO-
MOM approach offers an improvement over SCF-based methods49, 24, 51,
variational collapse can still be an issue, as will be demonstrated below.

These challenges call for the development of a practical and robust method
specifically designed to find a saddle point of a given order on the
electronic energy surface. In the context of finding the mechanism and rate of
atomic rearrangements, such as chemical reactions and diffusion events,
several methods have been developed for finding 1st-order saddle points on
energy surfaces that describe the variation of the ground electronic state
energy as a function of atomic coordinates.71, 72, 73 Convergence on a 1st-
order saddle point can be achieved by taking a step uphill in energy in the
direction of the eigenvector corresponding to the lowest eigenvalue of the
Hessian matrix and downhill in all perpendicular directions. This approach can
be generalized to find higher-order saddle points by taking an uphill step
along more than one eigenvector. If the Newton-Raphson algorithm is used,
second-order local convergence can be achieved 74, 75, 76. Such algorithms
have been developed for state-specific MCSCF calculations of excited states,
typically involving the trust-region method for the Newton-Raphson step 77,
78, 39, 41. This method is based on the evaluation of second derivatives of
the energy, thereby involving large computational effort making it impractical
for large systems. The step control with trust radius for excited states can,
furthermore, be problematic39. Quasi-Newton methods use approximate Hessian
update techniques, such as the symmetric rank-one (SR1),79 Powell,80 and
Bofill81 updates, and require less computational effort. This approach has
been used in DO-MOM calculations of excited states, 49, 24, 51 but it relies
on having a good initial guess for the Hessian to construct a preconditioner.

An alternative approach for finding 1st-order saddle points that does not
involve constructing or estimating the Hessian and still provides fast
convergence is the minimum mode following method 82, 83, 84, 85. There, only
the eigenvector of the Hessian corresponding to the lowest eigenvalue
(referred to as the minimum mode) needs to be evaluated, which can be
accomplished, for example, by using the dimer 82, Lánczos 86, 83 or Davidson
44, 85 methods. The projection of the gradient on the direction of the minimum
mode is then inverted and the optimization performed using any gradient-based
minimization technique.

Here, the minimum mode following method is generalized to target an $n$th-
order saddle point by following the eigenvectors of the $n$ lowest eigenvalues
instead of just the minimum mode. The needed eigenvectors are evaluated using
a generalized Davidson algorithm87, thereby requiring only first-order
derivatives of the energy. The resulting generalized mode following (GMF)
approach is combined with the exponential transformation direct optimization
of the orbitals to yield the DO-GMF method. It can be used to calculate
excited states corresponding to saddle points without the need of MOM or other
special procedures to reduce the risk of variational collapse. To illustrate
the methodology and demonstrate its performance, DO-GMF calculations of
excited states with the PBE generalized gradient approximation functional54
are presented here. The convergence properties are compared to those of the
DO-MOM method, which has previously been shown to outperform SCF-MOM-based
strategies49, 24, 51. First, the energy curves of the doubly excited states of
the dihydrogen and ethylene molecules along the single bond stretching and
double bond torsion, respectively, are calculated. These calculations are
performed with DO-GMF by targeting a 2nd-order saddle point, while DO-MOM
converges on a lower-energy solution after the SBO, giving a qualitatively
incorrect description of the potential energy surface for atomic dynamics. The
performance of DO-GMF is further assessed by calculating charge transfer
excited states of the nitrobenzene and twisted N-phenylpyrrole molecules as
well as the bond dissociation energy curve of the first excited state of a
large diplatinum and silver complex, [AgPt2(P2O5H2)4]3- (AgPtPOP). The DO-GMF
approach is found to be more robust than the DO-MOM method in calculations of
excited states using density functionals and furthermore, ensures convergence
to a solution providing potential energy curves consistent with the target
excited state for atomic dynamics at SBOs.

The present article is structured as follows. In section 2, the methodology is
presented, including a summary of direct orbital optimization and the saddle
point search algorithm. The correspondence between excited states and saddle
points on the electronic energy surface is illustrated for the minimal-basis
H2 molecule. Section 3 gives information on the implementation and parameters
used in the density functional calculations. Results are presented in section
4, where the ability of DO-GMF to track excited state solutions through SBOs
is demonstrated by calculations of potential energy curves for the dihydrogen
and ethylene molecules. Additionally, the performance of DO-GMF and DO-MOM is
compared in calculations of charge transfer excitations of nitrobenzene and
twisted N-phenylpyrrole as well as in calculations of the energy curve of the
first excited state of AgPtPOP. The article concludes with a discussion of the
results and a future outlook in section 5.

##  2 Methodology

###  2.1 Direct Orbital Optimization

DO24, 88 is an alternative to the SCF approach. In DO, one seeks a unitary
matrix $\mathbf{U}$ transforming an initial guess (or reference) electronic
wave function $\bm{\psi}_{0}$ into the wave function
$\bm{\psi}_{\mathrm{stat}}$ that represents a stationary point on the
electronic energy surface. Within the linear combination of atomic orbitals
(LCAO) formalism, the molecular orbitals,
$\bm{\psi}=\left(\ket{\psi_{1}},\ldots,\ket{\psi_{M}}\right)^{T}$ , are
expanded in a set of $M$ basis functions,
$\bm{\phi}=\left(\ket{\phi_{1}},\ldots,\ket{\phi_{M}}\right)^{T}$ , as

$\bm{\psi}=\mathbf{C}\bm{\phi}\,,$ (1)

where $\mathbf{C}$ is the $M\times M$ matrix of coefficients. Thus, the
objective is to find the optimal orbital coefficients
$\mathbf{C}_{\mathrm{stat}}$ that make the energy stationary by applying a
unitary transformation to a set of reference orbitals $\mathbf{C}_{0}$ ,

$\textbf{C}_{\mathrm{stat}}=\mathbf{C}_{0}\mathbf{U}\,.$ (2)

In the Hartree-Fock and Kohn-Sham (KS) approaches, the ground state is the
global minimum of the energy as a function of $\mathbf{U}$, while excited
electronic states are typically saddle points.37, 6, 24, 42

####  2.1.1 Exponential Transformation

The unitary matrix can conveniently be expressed as the exponential of a
matrix $\mathbf{\Theta}$ containing pairwise orbital rotation angles. The
orbital orthonormality constraints are imposed by requiring $\mathbf{\Theta}$
to be anti-Hermitian,

$\mathbf{U}=e^{\mathbf{\Theta}},\mathrm{\ subject\ to\
}\mathbf{\Theta}=-\mathbf{\Theta}^{{\dagger}}\,.$ (3)

A stationary point in the linear space of anti-Hermitian matrices directly
corresponds to a stationary point in the non-linear space of the corresponding
exponentials. Therefore, the problem is recast into making the energy
stationary as a function of the elements of $\mathbf{\Theta}$, which can be
decomposed into three unique blocks containing rotations among the occupied
orbitals (oo), rotations among the virtual orbitals (vv), and rotations
between the occupied and virtual orbitals (ov),

$\mathbf{\Theta}=\begin{pmatrix}\mathbf{\Theta}_{\mathrm{oo}}&\mathbf{\Theta}_{\mathrm{ov}}\\\
-\mathbf{\Theta}_{\mathrm{ov}}^{{\dagger}}&\mathbf{\Theta}_{\mathrm{vv}}\end{pmatrix}\,.$
(4)

The energy gradient with respect to rotations in the vv block is always 0,
meaning that this block is constant during an optimization and can without
loss of generality be set to 0. The oo block can be set to 0 if the density
functional is unitary invariant. However, if orbital density dependence is
incorporated into the functional, e.g. by including the Perdew-Zunger self-
interaction correction89, the functional is no longer unitary invariant,
meaning that the optimal $\mathbf{\Theta}_{\text{oo}}$ is a unique anti-
Hermitian matrix.90 The ov block always has to be optimized. The size of the
oo, ov and vv blocks is $N\times N$, $N\times(M-N)$ and $(M-N)\times(M-N)$,
respectively, $N$ being the number of occupied orbitals, though only one
triangular part of the oo block is unique for orbital density dependent
functionals.

The optimization is performed in the linear space of anti-Hermitian matrices
by applying standard numerical optimization techniques, such as quasi-Newton
methods. The optimization method used for excited state calculations must be
able to converge on saddle points. The components of the gradient are

$\frac{\partial
E}{\partial\Theta_{ij}}=\frac{2-\delta_{ij}}{2}\left[\int_{0}^{1}dte^{t\mathbf{\Theta}}\mathbf{L}e^{-t\mathbf{\Theta}}\right]_{ji}=\frac{2-\delta_{ij}}{2}\left[\mathbf{L}+\frac{1}{2!}\left[\mathbf{\Theta},\mathbf{L}\right]+\frac{1}{3!}\left[\mathbf{\Theta},\left[\mathbf{\Theta},\mathbf{L}\right]\right]+...\right]_{ji}$
(5)

with

$L_{ij}=\left(f_{i}-f_{j}\right)H_{ij}\,,$ (6)

where $H_{ij}$ are the elements of the Hamiltonian matrix in the basis of
molecular orbitals and $f_{i}$ are the orbital occupation numbers. The right-
hand side of eq 5 is a special case of the Baker-Campbell-Hausdorff formula91,
92, 93, 94 and makes use of the fact that the anti-Hermitian matrices form the
Lie-algebra corresponding to the Lie-group of unitary matrices. If the norm of
$\mathbf{\Theta}$ is small (${\|\mathbf{\Theta}\|}\ll 1$), the series can be
truncated at the first term ($\frac{\partial
E}{\partial\Theta_{ij}}\approx\frac{2-\delta_{ij}}{2}L_{ji}$), providing an
efficient way of evaluating the gradient. The elements of $\mathbf{\Theta}$
can be kept small by updating the reference orbitals and setting all orbital
rotations to 0 at regular step intervals. For a unitary invariant functional,
the following diagonal approximation to the Hessian can be used as
preconditioner,

$\frac{\partial^{2}E}{\partial\Theta_{ij}^{2}}\approx
2\left(f_{j}-f_{i}\right)\left(\epsilon_{i}-\epsilon_{j}\right)\,,$ (7)

where $f_{i}$ and $\epsilon_{i}$ are the orbital occupation numbers and
eigenvalues of the Hamiltonian matrix, respectively.

An excited state calculation is usually initialized using the ground state
orbitals after swapping the occupation numbers of the occupied and unoccupied
orbitals involved in the excitation. For example, in a HOMO-LUMO excitation,
the HOMO gets assigned an occupation number of 0 while the LUMO gets assigned
an occupation number of 1. Eq 7 is then used with the resulting non-aufbau
occupations and the energies of the ground state orbitals. The preconditioner
in an excited state calculation has one negative eigenvalue for each pair of
occupied-unoccupied orbitals where the energy of the unoccupied orbital is
smaller than the energy of the occupied orbital.

###  2.2 Excited states and saddle points on the energy surface

The correspondence between saddle points on the electronic energy surface and
excited electronic states is illustrated in Figure 1 using the H2 molecule
with a minimal basis set as an example. The calculations are carried out with
the PBE functional54 in the spin-unrestricted formalism.

Figure 1:  Energy of the three states with one electron in each spin channel,
$\alpha$ and $\beta$, that can be obtained for the minimal-basis H2 molecule
as a function of the atom distance. The orbitals are related by rotation
angles $\sigma$ and $\tau$ with respect to the orbitals of the symmetry-pure
ground state solution, S${}_{0}^{\mathrm{S}}$, with configuration
$\sigma_{g}^{2}\sigma_{u}^{*0}$:
$\psi^{\alpha}_{0,\mathrm{ref}}=\psi^{\beta}_{0,\mathrm{ref}}=\sigma_{g}$ and
$\psi^{\alpha}_{1,\mathrm{ref}}=\psi^{\beta}_{1,\mathrm{ref}}=\sigma_{u}^{*}$,
where $\sigma_{g}$ and $\sigma_{u}^{*}$ are symmetry-adapted bonding and
antibonding molecular orbitals, respectively. The subscripts $0$ and $1$
indicate occupied and unoccupied orbitals, respectively. The ground state
corresponds to a minimum on the electronic energy surface shown for the
optimal bond length, $r_{e}$, in the left contour graph, while the first and
doubly excited state solutions, S${}_{1}^{\mathrm{BS}}$ and
S${}_{2}^{\mathrm{S}}$, correspond to 1st\- and 2nd-order saddle points,
respectively. Beyond a distance of 1.1 Å, the spatial symmetry can break in
the second excited state to give a higher-energy solution,
S${}_{2}^{\mathrm{BS}}$, corresponding to charge transfer between the two H
atoms. Beyond a distance of 1.7 Å, the spin symmetry can break in the ground
state to give a lower-energy solution, S${}_{0}^{\text{BS}}$, corresponding to
spin polarization between the two H atoms. The right contour graph shows the
energy surface at $r_{e}+1.5$ Å, which is beyond the two symmetry breaking
points. There, the extrema have each split into two, with the minima and
maxima corresponding to symmetry-broken ground and doubly excited state
solutions, S${}_{0}^{\mathrm{BS}}$ and S${}_{2}^{\mathrm{BS}}$, respectively,
while the symmetry-pure solutions, S${}_{0}^{\mathrm{S}}$ and
S${}_{2}^{\mathrm{S}}$, correspond to 1st-order saddle points.

There are only two molecular orbitals for each spin channel in this case for
the two electrons. The figure shows the potential energy curves for all
solutions that can be obtained with one electron in each spin channel. Contour
graphs of the electronic energy surface are also shown to illustrate how the
energy depends on the rotation angles $\tau$ and $\sigma$ that mix the
occupied and unoccupied orbitals in each spin channel. One contour graph is
shown for the equilibrium bond length and another for a stretched bond where
symmetry has been broken in both the ground state and the doubly excited
state. Each point on this electronic energy surface is obtained according to
the transformation

$\begin{pmatrix}\vspace{1pt}\psi_{0}^{\alpha}\\\
\vspace{1pt}\psi_{1}^{\alpha}\\\ \vspace{1pt}\psi_{0}^{\beta}\\\
\vspace{1pt}\psi_{1}^{\beta}\end{pmatrix}=\begin{pmatrix}\cos{\tau}&\sin{\tau}&0&0\\\
-\sin{\tau}&\cos{\tau}&0&0\\\ 0&0&\cos{\sigma}&\sin{\sigma}\\\
0&0&-\sin{\sigma}&\cos{\sigma}\end{pmatrix}\begin{pmatrix}\vspace{1pt}\psi_{0,\,\mathrm{ref}}^{\alpha}\\\
\vspace{1pt}\psi_{1,\,\mathrm{ref}}^{\alpha}\\\
\vspace{1pt}\psi_{0,\,\mathrm{ref}}^{\beta}\\\
\vspace{1pt}\psi_{1,\,\mathrm{ref}}^{\beta}\end{pmatrix}\,,$

where the indices 0 and 1 correspond to occupied and unoccupied orbitals,
respectively, and the indices $\alpha$ and $\beta$ correspond to the two spin
channels. The reference orbitals (subscript “ref”) are the orbitals of the
symmetry-pure ground state solution, S${}_{0}^{\text{S}}$, with configuration
$\sigma_{g}^{2}\sigma_{u}^{*0}$, where $\sigma_{g}$ and $\sigma_{u}^{*}$ are
the symmetry-adapted bonding and antibonding molecular orbitals, respectively.
At the ground state equilibrium geometry, S${}_{0}^{\text{S}}$ corresponds to
the global minimum on the electronic energy surface. There are two other types
of solutions at this geometry: a solution with broken spin symmetry,
S${}_{1}^{\text{BS}}$, representing a singly excited open-shell state with
configuration $\sigma_{g}^{1}\sigma_{u}^{*1}$ corresponding to a 1st-order
saddle point obtained by $\pm 90^{\circ}$ rotation in one of the spin
channels, and a symmetry-pure solution representing a doubly excited state
with configuration $\sigma_{g}^{0}\sigma_{u}^{*2}$ corresponding to a 2nd-
order saddle point obtained by $\pm 90^{\circ}$ rotation in both spin
channels. For each symmetry-pure state, there are two degenerate solutions
related by a sign change, while there are four degenerate solutions for each
symmetry-broken state additionally related by a symmetry operation.

As the bond is stretched, SBOs appear from which two additional types of
solutions emerge. One is a solution with broken spin symmetry,
S${}_{0}^{\text{BS}}$, representing a ground state with diradical character
(H↑H↓/H↓H↑). The SBO for this symmetry-broken ground state solution is often
referred to as a Coulson-Fischer point60. S${}_{0}^{\text{BS}}$ corresponds to
the global minimum on the electronic energy surface, and as the H-H distance
tends to infinity, it can be obtained from the symmetry-pure
S${}_{0}^{\text{S}}$ solution by applying $45^{\circ}$ orbital rotations of
opposite sign in the two spin channels. The other symmetry-broken solution,
S${}_{2}^{\text{BS}}$, represents a doubly excited state with ionic character
(H+H-/H-H+) where the spatial symmetry is lost. It is a 2nd-order saddle point
accessible through $45^{\circ}$ orbital rotations of the same sign in both
spin channels as the H-H distance tends to infinity. The symmetry-broken
solutions, S${}_{0}^{\text{BS}}$ and S${}_{2}^{\text{BS}}$, arise as
degenerate pairs related by a swap of the spin channels and inversion at the
center of the molecule, respectively. They thereby cease to be eigenfunctions
of those symmetry operations.

The symmetry-broken states are the preferred solutions in a mean-field
calculation as they effectively take static correlation into account61, 66. In
a calculation of an energy curve representing the potential energy profile for
atomic dynamics, it is important to identify which solution to follow at SBOs.
As can be seen from Figure 1, the saddle point order for the symmetry-broken
solutions is the same as that of the symmetry-pure counterpart before the
SBOs. Following a saddle point of a given order, therefore, smoothly tracks
the symmetry-broken solution through an SBO. The symmetry-pure
S${}_{0}^{\text{S}}$ and S${}_{2}^{\text{S}}$ solutions, however, become 1st-
order saddle points beyond the SBOs. Therefore, the SBOs represent inflection
points on the energy curves of the symmetry-pure solutions (a sign change of
the first and second Hessian eigenvalues for S${}_{0}^{\text{S}}$ and
S${}_{2}^{\text{S}}$, respectively), while no inflection points are present
for the symmetry-broken solutions (persistent sign of the first/second Hessian
eigenvalue). The symmetry-broken solutions provide potential energy curves for
the bond dissociation with qualitatively correct shape and asymptotic
behavior95, while the symmetry-pure solutions overestimate the energy of the
ground state and underestimate the energy of the doubly excited state, the two
becoming degenerate at infinite H-H distance.

####  2.2.1 Computing excited state energy surfaces

There are two ways of mapping out excited state potential energy surfaces for
the dynamics of atoms in molecules, such as those illustrated in Figure 1. The
strategy that is most commonly employed in geometry optimizations and
classical dynamics simulations is a sequential point acquisition scheme. Here,
one first performs a ground state calculation at the initial geometry and uses
the obtained orbitals with occupation numbers changed to reflect the character
of the target excited state as an initial guess for the excited state wave
function optimization. Subsequent points are generated by displacing the atoms
and using the orbitals of the previous converged solution as the initial
guess. The second strategy is a separate point acquisition scheme, where for
each configuration of the atoms the ground state orbitals are used as the
initial guess for the excited state wave function calculation. The sequential
scheme is more efficient, as it does not require a ground state calculation
for each geometry and, moreover, can converge faster if the atomic
displacements are sufficiently small. However, as mentioned above, a problem
can arise when multiple solutions emerge from one at an SBO, if the
calculation converges to a solution that gives a qualitatively incorrect
potential energy surface. As illustrated by calculations of the excited state
of H2 discussed above and also by our previously published calculations of the
torsional potential energy curves for the ethylene molecule53, the multiple
excited state solutions correspond to saddle points of different orders. For
excited electronic states, one cannot simply select the state that has lower
energy, as for the ground state, but a possibility is to select the state that
corresponds to the same saddle point order before and after the SBO. The
mapping of excited state potential energy surfaces, therefore, requires a
method for reconverging on a saddle point of a given order after the atomic
configuration has been changed. Below we present such a method, the
generalized mode following approach.

###  2.3 Generalized Mode Following

####  2.3.1 Minimum Mode Following

Minimum mode following is a saddle point optimization technique used to
determine 1st-order saddle points on the energy surface in atomic
configuration space 82, 83, 85. The method recasts the challenging saddle
point search as a minimization by inverting the projection of the gradient
vector $\bm{g}$ on the unstable Hessian mode $\bm{v}_{1}$ at each iteration of
the optimization,

$\displaystyle\bm{g}^{\,\parallel}$
$\displaystyle=\bm{v}_{1}\bm{v}_{1}^{\mathrm{T}}\bm{g}\,,$ (8)
$\displaystyle\bm{g}^{\mathrm{\,mod}}$
$\displaystyle=\bm{g}-2\bm{g}^{\,\parallel}$ $\displaystyle\mathrm{\ if\
}\lambda_{1}<0\,,$ (9) $\displaystyle\bm{g}^{\mathrm{\,mod}}$
$\displaystyle=-\bm{g}^{\,\parallel}$ $\displaystyle\mathrm{\ if\
}\lambda_{1}\geq 0\,.$ (10)

The perpendicular component of the force is not included (eq 10) if the
curvature $\lambda_{1}$ along the lowest eigenvector is positive, i.e. if the
function is convex in this direction. This procedure yields a gradient that
leads to a minimum in a revised objective function where the original function
has a 1st-order saddle point, so any gradient-based minimization technique can
then be used to converge on the saddle point. However, since the modified
objective function is not known, line search techniques commonly used to
accelerate convergence cannot straightforwardly be applied.

The lowest curvature Hessian mode can be obtained with a variety of partial
diagonalization methods. The most common ones are the dimer82, 84, 96,
Davidson44, 87 and Lánczos86 methods, where the relevant parts of the Hessian
are evaluated by a finite difference approximation. The Davidson method can be
regarded as a preconditioned Lánczos method. The preconditioner is typically
diagonal but a full approximate Hessian can also be used.

####  2.3.2 Generalization to saddle points of arbitrary order

While the estimation of rates of atomic rearrangements involves finding 1st-
order saddle points, the calculation of excited electronic states can require
finding saddle points of higher order. The mode following method can be
generalized to find a saddle point of order $n$ by identifying the modes
corresponding to the $n$ lowest eigenvalues and following them simultaneously.
The gradient projections along the $n$ modes are summed up to yield the
parallel gradient which is then inverted in the usual way,

$\displaystyle\bm{g}^{\,\parallel}$
$\displaystyle=\sum_{i=1}^{n}\bm{v}_{i}\bm{v}_{i}^{\mathrm{T}}\bm{g}\ \ \
{\mathrm{and}}\ \ \ \bm{g}^{\mathrm{\,mod}}=\bm{g}-2\bm{g}^{\,\parallel}$
$\displaystyle\mathrm{\ if\ }\lambda_{n}<0\,,$ (11)
$\displaystyle\bm{g}^{\,\parallel}$
$\displaystyle=\sum_{\begin{subarray}{c}i=1\\\\[4.0pt] \lambda_{i}\geq
0\end{subarray}}^{n}\bm{v}_{i}\bm{v}_{i}^{\mathrm{T}}\bm{g}\ \ \
{\mathrm{and}}\ \ \ \bm{g}^{\mathrm{\,mod}}=-\bm{g}^{\,\parallel}$
$\displaystyle\mathrm{\ if\ }\lambda_{n}\geq 0\,.$ (12)

Note that the perpendicular component of the gradient is only included in eq.
11 if eigenvalue $n$ is negative, i.e. if the electronic energy surface is
concave along all the target eigenvectors. If this condition is not satisfied,
only the target eigenvectors corresponding to non-negative eigenvalues are
followed (eq. 12).

####  2.3.3 Implementation of GMF

The implementation makes use of the generalized Davidson method presented in
ref 87 to obtain the lowest $n$ eigenvectors of the electronic Hessian when
targeting a saddle point of order $n$. At the first step of the wave function
optimization, a diagonal approximation $\mathbf{D}$ to the Hessian is computed
according to eq 7. Then, the generalized Davidson algorithm starts by defining
an initial matrix $\mathbf{K}$ having $n$ unit vectors $\bm{k}_{i}$ as columns
along the orbital rotations corresponding to the lowest $n$ eigenvalues of the
approximate Hessian (a Krylov subspace). A small random perturbation is
applied to this initial Krylov subspace and the subspace is orthonormalized
using the modified Gram-Schmidt approach. At each iteration of the Davidson
algorithm, the effect of the electronic Hessian matrix $\mathscrbf{H}$ on the
Krylov subspace is evaluated by a forward finite difference approximation

$\mathscrbf{H}\bm{k}_{j}\approx\frac{\nabla
E\left(\mathbf{C}e^{h\mathbf{\Theta}\left[\bm{k}_{j}\right]}\right)-\nabla
E\left(\mathbf{C}\right)}{h}\,,$ (13)

where $\mathbf{\Theta}\left[\bm{k}_{j}\right]$ is the anti-Hermitian matrix
having the elements of the $j$th vector of the Krylov subspace, $\bm{k}_{j}$,
in the upper triangular part, $h$ is the finite difference step size and
$\nabla E\left(\mathbf{C}\right)$ is the energy gradient vector $\bm{g}$
evaluated using the current LCAO coefficient matrix $\mathbf{C}$. Next, the
Rayleigh matrix $\mathbf{K}^{\mathrm{T}}\mathscrbf{H}\mathbf{K}$ is
diagonalized to obtain the lowest $n$ eigenpairs ($\lambda_{i}$,
$\bm{y}_{i}$). The Ritz vectors $\bm{x}_{i}=\mathbf{K}\bm{y}_{i}$ represent
approximations to the target eigenvectors. The residual vectors
$\bm{r}_{i}=(\lambda_{i}\mathbf{I}-\mathscrbf{H})\bm{x}_{i}$ are multiplied by
a preconditioner $\mathbf{P}_{i}$, and the Krylov subspace is extended by
incorporating the resulting $n$ vectors $\mathbf{P}_{i}\bm{r}_{i}$. The
preconditioner presented in ref 97 is used,

$\mathbf{P}_{i}=(\lambda_{i}\mathbf{I}-\mathbf{D})^{-1}\,,$ (14)

where $\mathbf{D}$ is the diagonal approximation to the Hessian given by eq 7.
A threshold of $-0.1$ Ha is applied to the components of the preconditioner.
If it is exceeded for a given element, the element is set to this threshold to
ensure that the preconditioner is negative definite and that the Davidson
method converges to the lowest eigenpairs. At a given iteration, only the
effect of the Hessian on the at most $n$ vectors added to the Krylov subspace
in the previous iteration needs to be evaluated. If the dimensionality of the
Krylov subspace becomes too large, a new Krylov subspace is formed by
including only the $n$ approximate eigenvectors $\bm{x}_{i}$ and
preconditioned residual vectors so that the cost of the subspace
eigendecomposition is always small. The residuals $\bm{r}_{i}$ are also used
to check for convergence. In the present implementation, the convergence
threshold is a maximum component of 0.01 Ha. If a target eigenpair is
converged in a given iteration, the Krylov subspace is not extended in the
preconditioned direction of the residual. Once all eigenpairs have converged
and a step in the wave function optimization has been taken using the modified
gradient (eq 11 or 12), the next Davidson cycle is accelerated by using the
eigenvectors found at the previous wave function optimization step to form the
initial Krylov subspace.

Figure 2: Flowchart of the DO-GMF algorithm. The approach consists of a direct
optimization outer loop (green) using the exponential transformation and a
quasi-Newton step with a modified gradient determined in a partial Hessian
diagonalization inner loop (blue) using the generalized Davidson method.

As GMF recasts the saddle point search as a minimization, there is no risk of
convergence on a lower-order saddle point or the ground state minimum (i.e. a
variational collapse) and therefore, methods designed to reduce such problems,
in particular MOM, are not needed. Additionally, there is no need for a quasi-
Newton Hessian update with the flexibility to develop an indefinite matrix.
Instead, one can use robust update formulas for minimization, such as the
efficient limited-memory BFGS (L-BFGS) algorithm. The implementations of the
exponential transformation and the quasi-Newton algorithms, including L-BFGS,
for updating the orbitals in the LCAO representation are presented in refs 24
and 50. An overview of the DO-GMF method is provided in the flowchart in
Figure 2.

The DO-GMF method has been implemented in a development branch of the grid-
based projector augmented wave (GPAW)98, 99 software using the exponential
transformation24 to update the orbitals represented in an LCAO basis100.

###  2.4 Specifics of the electronic structure calculations

The calculations make use of the generalized gradient approximation (GGA)
functional PBE, apart from the calculations on the AgPtPOP complex, which
instead use the BLYP functional. Valence electrons are represented by an LCAO
basis set consisting of primitive Gaussian functions taken from the aug-cc-
pVDZ101, 102, 103 (dihydrogen, ethylene), def2-TZVPD104 (nitrobenzene), or cc-
pVDZ (N-phenylpyrrole) sets augmented with a single set of numeric atomic
orbitals105, 100. For the AgPtPOP complex, a double zeta polarized (dzp) basis
set of numeric atomic orbitals100 is used. Real orbitals are used. The frozen
core approximation is used within the PAW approach106.

The ground state calculations as well as the DO-MOM and DO-GMF excited state
calculations are carried out with the exponential transformation direct
optimization method 24, 49, 50 using a limited-memory BFGS (L-BFGS) algorithm
with inexact line search for the ground state50, a limited-memory SR1 (L-SR1)
algorithm with a maximum step length of 0.2 for the excited state DO-MOM
calculations24 (default maximum step length in GPAW), and L-BFGS with a
maximum step length of 0.2 for the excited state DO-GMF calculations, unless
otherwise stated. The DO-MOM calculations use the MOM method as presented in
ref 9, where at each wave function optimization step, the occupied orbitals
are chosen as those with the largest projections 24

$\omega_{j}=\sqrt{\sum_{i=1}^{N}|O_{ij}|^{2}}\,,$ (15)

where $O_{ij}$ are the elements of the overlap matrix between the orbitals at
the current step and at the initial step

$\mathbf{O}=\mathbf{C}_{0}^{{\dagger}}\mathbf{S}\mathbf{C}$ (16)

with $\mathbf{S}$ being the overlap matrix of basis functions. No
orthogonality constraints to lower-energy states are enforced, so the
calculations are fully variational. The DO-MOM and DO-GMF calculations are
converged to a precision of $10^{-7.4}\,\mathrm{eV}^{2}$ per valence electron
in the squared residual of the KS equations and $10^{-7}$ eV per valence
electron in the three-step energy change. The atomic configuration of the
nitrobenzene molecule is taken from ref 32, while the N-phenylpyrrole molecule
is from ref 107. The atomic configuration of the ethylene molecule and the
AgPtPOP complex are first optimized in the ground state at the same level of
theory as used in the excited state calculations. The potential energy curves
are then evaluated by incrementing a structural parameter while keeping the
values of all other parameters the same as in the ground state optimal
configuration. All open-shell singlet excited states are spin-mixed states,
and spin purification is not applied to the energy. The calculations are
performed with the GPAW software98, 99 and LIBXC108 version 4.0.4. The grid
spacings are 0.14 Å (dihydrogen with minimal basis set, AgPtPOP), 0.2 Å
(dihydrogen, ethylene) and 0.15 Å (nitrobenzene, N-phenylpyrrole), while the
dimensions of the simulation cell are according to the default cutoff of the
numeric representation of the basis functions105.

###  2.5 Advantage of DO-GMF for finding symmetry-broken solutions

The advantage of the DO-GMF method over, for example, DO-MOM is illustrated in
Figure 3. There, two calculations of the doubly excited state of minimal-basis
H2 using the PBE functional are shown. The first corresponds to the separate
point acquisition scheme, and the second corresponds to the sequential one
when generating a potential energy surface for atomic dynamics.

In the first calculation, illustrated in Figure 3(a), a single point
calculation for an H-H distance of 0.4 Å larger than the minimum energy
distance $r_{e}$ is carried out. This distance is beyond the SBO, as can be
seen from the contour graph of the electronic energy surface, where the
maximum has split up into two maxima. The calculation is initialized using the
symmetry-pure ground state solution, S${}_{0}^{\text{S}}$, by promoting both
electrons from the occupied to the unoccupied orbitals, i.e. 90∘ orbital
rotations, and then a 2nd-order saddle point is targeted. The initial guess
turns out to be located at the 1st-order saddle point to within numerical
accuracy, thereby corresponding to the symmetric doubly excited state,
S${}_{2}^{\text{S}}$. Despite the initial close proximity to a 1st-order
saddle point, the DO-GMF method climbs up the energy surface, converging on
one of the two equivalent 2nd-order saddle points that correspond to symmetry-
broken S${}_{2}^{\text{BS}}$ solutions and have higher energy than the
symmetric S${}_{2}^{\text{S}}$ state. However, a DO-MOM calculation using an
L-SR1 quasi-Newton optimizer24 converges right away on the 1st-order saddle
point when starting from this initial guess and is not able to climb up to the
2nd-order saddle point.

Figure 3: Illustration of the way DO-GMF calculations can converge on
symmetry-broken solutions and provide a potential energy curve for atomic
dynamics simulations, while DO-MOM converges on the symmetry-pure solutions.
Results of calculations on the minimal-basis H2 molecule in the doubly excited
singlet state are shown using the PBE functional. (a) Single point
calculations with DO-GMF and DO-MOM at $r_{e}+0.4$ Å where symmetry breaking
can occur. The initial guess (cross) is generated by double excitation from
the symmetry-pure ground state solution, S${}_{0}^{\text{S}}$. The DO-MOM
calculation (red line) converges to the 1st-order saddle point (red diamond)
corresponding to the symmetry-pure solution, S${}_{2}^{\text{S}}$, as it is
nearly at the same location on the electronic energy surface as the initial
guess, while the DO-GMF calculation climbs up (purple line) to one of the two
equivalent second order saddle points (purple diamond), corresponding to a
symmetry-broken solution, S${}_{2}^{\text{BS}}$. The origin of the contour
graph of the electronic energy surface is at the symmetry-pure solution,
S${}_{2}^{\text{S}}$. The lower graph shows the energy curves for
S${}_{2}^{\text{BS}}$ (purple curve) and S${}_{2}^{\text{S}}$ (red curve). For
short H-H distance, only the symmetry-pure solution exists (orange curve). (b)
Sequential calculations of six points spaced by 0.1 Å along the energy curves
shown in (a) using DO-GMF (purple) and DO-MOM (red). The initial guess at the
first point at $r_{e}$ is obtained as in (a), while subsequent calculations
use the orbitals of the previous H-H distance (indicated by dashed lines).
Even after the onset of symmetry breaking at $r_{e}+0.1$ Å, DO-MOM keeps
converging on the symmetry-pure solution, S${}_{2}^{\text{S}}$, a 1st-order
saddle point, while DO-GMF converges consistently on 2nd-order saddle points
corresponding to symmetry-broken solutions, S${}_{2}^{\text{BS}}$, and thereby
provides a more accurate energy curve for atomic dynamics in the excited
state.

The second example, illustrated in Figure 3(b), is a sequential calculation of
six points along the potential energy curve for the second excited state of
the minimal-basis H2 molecule. The H-H distance is increased by 0.1 Å at each
step starting from the optimal bond length in the ground electronic state,
$r_{e}$. Again, the calculations target a 2nd-order saddle point to find the
second excited state. The calculation is initialized by performing 90∘ orbital
rotations from the symmetry-pure ground state solution, S${}_{0}^{\text{S}}$
at $r_{e}$. This initial guess is close to the 2nd-order saddle point
corresponding to the symmetry-pure solution, S${}_{2}^{\text{S}}$, the only
doubly excited state solution that exists before the SBO. Both DO-GMF and DO-
MOM converge on this stationary point. Then, the H-H distance is incremented
to $r_{e}+0.1$ Å which is beyond the SBO, and two types of solutions emerge: a
1st-order saddle point corresponding to the symmetry-pure doubly excited state
solution, S${}_{2}^{\text{S}}$, and a pair of equivalent 2nd-order saddle
points corresponding to the spatially symmetry-broken solution,
S${}_{2}^{\text{BS}}$. The 1st-order saddle point on the energy surface for
$r_{e}+0.1$ Å is located at the same point on the electronic energy surface as
the 2nd-order saddle point on the energy surface for $r_{e}$ before the SBO.
Therefore, the initial guess in the calculation for $r_{e}+0.1$ Å is located
right at the 1st-order saddle on the electronic energy surface, so the DO-MOM
calculation converges there. However, DO-GMF is able to move away from the
1st-order saddle point and converges on one of the two 2nd-order saddle points
corresponding to a symmetry-broken solution. As the H-H distance is increased
further, the two equivalent 2nd-order saddle points move farther away from
each other, and DO-GMF moves on the electronic energy surface, while the DO-
MOM calculations repeatedly converge on the 1st-order saddle point
corresponding to the symmetry-pure solution.

These examples illustrate how convergence to a specific excited state solution
can be achieved with DO-GMF, thereby making it possible to generate a
potential energy surface for atomic dynamics systematically making use of the
advantage symmetry-broken solutions provide. Methods that do not guarantee
convergence on a saddle point of a given order have a tendency to converge on
the stationary solution closest to the initial guess, thereby producing the
less accurate symmetry-pure states.

##  3 Results

###  3.1 Energy curve for excited ethylene

A more challenging example of the way the DO-GMF method can produce a
potential energy curve for atomic dynamics in the presence of symmetry
breaking is given below in calculations of a doubly excited state of ethylene.
At a certain value of the torsional angle, symmetry-broken solutions emerge,
and they can be obtained by converging consistently on a 2nd-order saddle
point on the electronic energy surface.

Figure 4 shows the energy curve of ethylene calculated with both the DO-GMF
and DO-MOM methods as a function of the C=C torsional angle.

Figure 4: Calculated energy of a doubly excited state of ethylene as a
function of the torsional angle $\theta$ calculated with DO-GMF and DO-MOM
using sequential point acquisition. A double excitation from the ground state
is performed to initialize the excited state calculation at the first
geometry, $\theta=0^{\circ}$. The DO-GMF calculation targets a 2nd-order
saddle point. The points on the curves are colored according to the value of
the second eigenvalue of the electronic Hessian, $\lambda_{2}$, while the gray
vertical line marks where symmetry-broken solutions appear. Before that, both
DO-MOM and DO-GMF converge on the 2nd-order saddle point corresponding to the
symmetry-pure solution, $\pi^{0}\pi^{*2}$. After that, DO-MOM converges on a
1st-order saddle point corresponding to the symmetry-pure solution giving an
incorrect potential energy curve. Instead, the DO-GMF calculations keep
converging on a 2nd-order saddle point corresponding to a symmetry-broken
solution with ionic character (H2C+C-H2/H2C-C+H2), thereby providing a more
accurate potential energy curve.

The potential energy curves are calculated using sequential point acquisition
starting from the planar geometry of the ethylene molecule. The excited state
calculation at the first geometry is initialized by constructing a double
excitation from the ground state solution (90∘ HOMO-LUMO rotation in both spin
channels). The DO-GMF calculations always converge on a 2nd-order saddle point
on the electronic energy surface. To characterize the saddle point order of
the obtained solutions, Figure 4 shows the value of the second electronic
Hessian eigenvalue, $\lambda_{2}$, for each point along the potential energy
curve. As the torsional angle is increased towards 90∘, the double bond in
ethylene is broken. At the SBO, $\lambda_{2}$ becomes zero. Thereafter, DO-MOM
converges to a symmetry-pure solution with positive $\lambda_{2}$, thereby
corresponding to a 1st-order saddle point. The solution obtained with DO-MOM
loses another instability at a torsional angle of 120∘, where it coalesces
with a ground state solution with broken spin symmetry53. The symmetry-pure
solution obtained using DO-MOM beyond the SBO has covalent character,
$\pi^{0}\pi^{*2}$. This solution gives a diabatic energy surface unlike the
energy curve obtained from high-level multireference calculations53, 109, 95,
while DO-GMF converges on a 2nd-order saddle point corresponding to a solution
with ionic character (H2C+C-H2/H2C-C+H2), analogous to the
S${}_{2}^{\mathrm{BS}}$ solution of minimal-basis H2 shown in Figure 1. The
spatial symmetry is broken, and the potential energy curve displays a minimum
at a C=C torsional angle of $\theta=90^{\circ}$ in agreement with the
multireference results95. Therefore, by converging on a saddle point of the
appropriate order using DO-GMF, the calculations can give a qualitatively
correct potential energy curve for atomic dynamics.

###  3.2 Charge transfer excitations

Excitations involving large changes in electron density are especially prone
to variational collapse and can lead to convergence problems because they
often correspond to high-order saddle points on the electronic energy surface,
even when the calculation is initialized by an excitation involving ground
state orbitals close to the HOMO and LUMO.49, 24 As a proof of principle, the
performance of DO-GMF is first demonstrated in calculations of the open-shell
singlet $\pi^{*}\leftarrow\pi^{\prime}$ charge transfer excited state of
nitrobenzene (the notation is taken from ref. 110), which has successfully
been calculated previously with DO-MOM,24 but is known to be problematic for
SCF-MOM.49, 24, 51, 32, 110 Then, a more challenging calculation is presented,
a charge transfer excitation in orthogonally twisted N-phenylpyrrole107, where
DO-MOM is prone to convergence failure but DO-GMF performs well.

####  3.2.1 Nitrobenzene

Figure 5(a) shows the orbitals obtained in a ground state calculation of
nitrobenzene as well as the change in orbital occupation corresponding to a
singlet $\pi^{*}\leftarrow\pi^{\prime}$ charge transfer excited state. Charge
is transferred from the aromatic $\pi$-system to the nitro group.

Figure 5: (a) Ground state orbitals of the nitrobenzene molecule with
occupation numbers chosen to provide an initial guess for calculations of the
$\pi^{*}\leftarrow\pi^{\prime}$ charge transfer excited state, obtained using
the PBE functional and def2-TZVPD+sz basis set. (b) Orbitals of the converged
excited state obtained by converging on a 4th-order saddle point on the
electronic energy surface. The orbitals are rendered for an isovalue of $\pm
0.1\,\mathrm{\AA}^{-1.5}$.

Based on the ground state orbitals, the excitation is from the HOMO-2 to the
LUMO, indicating that an electron hole is formed below three occupied
orbitals. The excited state at the initial guess, before orbital relaxation,
thereby, appears to correspond to a 3rd-order saddle point, according to a
diagonal approximation of the Hessian. The preconditioner for the quasi-Newton
algorithm (eq 7) indeed has three negative diagonal elements. However, the
electronic Hessian at this initial guess based on the ground state orbitals
turns out to have 18 negative eigenvalues, many of them being close to zero
(five negative eigenvalues are larger than -1 eV). A better estimate of the
saddle point order corresponding to the target excited state is obtained by
freezing the HOMO-2 and the LUMO of the ground state and relaxing all other
orbitals. The resulting Hessian has five negative eigenvalues, one of them
being close to zero (larger than -1 eV). The charge transfer induces a
reordering of the orbitals, as the orbitals that are localized at the phenyl
group are stabilized, while orbitals localized at the nitro group are
destabilized. The final, converged excited state turns out to correspond to a
4th-order saddle point with the Hessian having four negative eigenvalues. The
orbitals for the converged excited state are shown in Figure 5(b). There, the
hole is below four occupied orbitals, rather than three as judged from the
ground state calculation.

Figure 6 shows the convergence of the energy and residual of the KS equations
in the DO-MOM and DO-GMF calculations starting from the same initial guess
(see also ref 24 for the DO-MOM calculations).

Figure 6: Convergence of DO-MOM and DO-GMF calculations of the open-shell
singlet $\pi^{*}\leftarrow\pi^{\prime}$ charge transfer excited state of
nitrobenzene measured in terms of the excitation energy in (a) and squared
residual of the KS equations in (b). The calculations use a maximum step
length of 0.2 in the quasi-Newton Hessian update (a default setting in GPAW).

DO-MOM with L-SR1 is able to converge on the excited state, meaning that the
L-SR1 Hessian update is in this case able to develop an additional negative
Hessian eigenvalue over the course of the optimization. DO-GMF using the
L-BFGS update and set to target a 4th-order saddle point converges on the same
excited state, but converges more smoothly than the DO-MOM calculation and in
slightly fewer iterations.

####  3.2.2 N-Phenylpyrrole

A second example of a charge transfer excitation is shown in Figure 7. The
promotion of an electron from the HOMO to the LUMO of the ground state wave
function of twisted N-phenylpyrrole corresponds to transfer of charge from the
pyrrole group ($\pi_{\mathrm{py}}$) to the phenyl group
($\pi^{*}_{\mathrm{ph}}$) as can be seen from the rendering of the orbitals.

Figure 7: (a) Ground state orbitals of N-phenylpyrrole with changed
occupations used as initial guess for DO-MOM and DO-GMF calculations of the
$\pi^{*}_{\mathrm{ph}}\leftarrow\pi_{\mathrm{py}}$ charge transfer excited
state using the PBE functional and a cc-pVDZ+sz basis set. (b) Orbitals of the
target charge transfer solution, which is reached by converging on a 6th-order
saddle point. (c) Orbitals of the 1st-order saddle point solution obtained in
a calculation with DO-MOM using the default maximum step length of 0.2 for the
quasi-Newton Hessian update. All orbitals are shown at an isovalue of $\pm
0.05\,\mathrm{\AA}^{-1.5}$.

Since this excitation appears to be a HOMO to LUMO transition based on the
ground state orbitals, the preconditioner at the initial guess has only one
negative element. Yet, the electronic Hessian has 28 negative eigenvalues,
with five being close to zero (larger than -1 eV). A relaxation of the
orbitals with the HOMO and LUMO frozen leads to a Hessian with seven negative
eigenvalues, one being close to zero (larger than -1 eV), indicating that the
saddle point on the energy surface corresponding to the converged excited
state is of 6th order. This prediction turns out to be correct, as the hole
gets stabilized during the variational optimization of the orbitals dropping
below three other occupied orbitals in energy, and furthermore, two other
unoccupied orbitals get stabilized with respect to the orbital to which the
electron is excited. The converged excited state has 6 negative Hessian
eigenvalues, thereby corresponding to a 6th-order saddle point.

Figure 8 shows the rate of convergence of the DO-GMF calculation and a
comparison with a DO-MOM calculation starting from the same initial guess for
the $\pi^{*}_{\mathrm{ph}}\leftarrow\pi_{\mathrm{py}}$ charge transfer
excitation.

Figure 8: Convergence of the excitation energy (a) and squared residual of the
KS equations (b) during DO-MOM and DO-GMF calculations of the open-shell
singlet $\pi^{*}_{\mathrm{ph}}\leftarrow\pi_{\mathrm{py}}$ charge transfer
excited state of N-phenylpyrrole. Both DO-MOM and DO-GMF calculations use a
maximum step length of 0.2 for the quasi-Newton Hessian update (default in
GPAW). While DO-MOM displays erratic convergence behavior and ultimately
collapses to a lower-energy 1st-order saddle point solution with small dipole
moment, DO-GMF rapidly converges to the target 6th-order saddle point
solution.

DO-GMF converges in a robust way in only 19 optimization steps and with a
monotonically decreasing energy and squared residual of the KS equations. The
DO-MOM calculation, however, using the default maximum step length of 0.2
collapses to a 1st-order saddle point and fails to converge on the targeted
charge transfer solution. This collapse is a consequence of the mismatch
between the number of concave directions estimated by the preconditioner and
that found at the converged charge transfer solution and an inability of the
L-SR1 update to develop the missing negative Hessian eigenvalues. The orbitals
of the 1st-order saddle point solution obtained with DO-MOM are visualized in
Figure 7c. The frontier orbitals contain a pair of occupied-unoccupied
orbitals delocalized over the entire molecule, which leads to a small dipole
moment. These delocalized orbitals arise from $\sim$45∘ mixing between the
$\pi_{\mathrm{py}}$ hole localized on the pyrrole group and a
$\pi_{\mathrm{ph}}$ occupied orbital localized on the phenyl group, giving
linear combination orbitals with characters
$\frac{1}{\sqrt{2}}(\pi_{\mathrm{py}}+\pi_{\mathrm{ph}})$ and
$\frac{1}{\sqrt{2}}(\pi_{\mathrm{py}}-\pi_{\mathrm{ph}})$ . MOM by
construction cannot prevent variational collapse for rotational angles smaller
than 45∘.

###  3.3 Bond dissociation in AgPtPOP

As an example of an application to a larger molecule, Figure 9 shows the
results of a calculation for a diplatinum and silver complex, AgPtPOP.

Figure 9: Energy of the AgPtPOP complex in the ground state (GS, blue) and in
the singlet excited states (ES, green) corresponding to a LUMO$\leftarrow$HOMO
excitation as a function of the Ag-Pt distance. The excited state points are
obtained using DO-GMF by targeting a 1st-order saddle point, using either a
sequential (green crosses) or separate acquisition (green circles). Results of
DO-MOM calculations (purple stars for the sequential acquisition and red
squares for the separate acquisition) are also shown. When the points are
calculated separately, the initial guess for an excited state calculation at a
given point is given by the orbitals of the ground state solution obtained for
the same Ag-Pt distance. In the sequential calculations, the initial guess at
a given point is given by the excited state orbitals obtained at the previous
point, starting from an Ag-Pt distance of 2.5 Å and then increasing and
decreasing the distance. The complex is visualized in the inset, showing the
relevant Ag-Pt distance with a dashed line, and the color code for the atoms
being: Ag = silver, Pt = blue, P = orange, O = red, H = white. The
calculations use the BLYP functional and a dzp basis set of numeric atomic
orbitals.

This compound derives from binding of an Ag+ ion with the d8-d8 diplatinum
complex [Pt2(P2O5H2)4]4– (PtPOP)111. The latter has attracted considerable
interest in recent years due to its photocatalytic properties as well as rich
photophysics, which has led to fundamental studies of photoinduced bond
vibration dynamics and energy relaxation112, 113, 36, 114. Upon
photoexcitation, an electron is promoted from an antibonding d$\sigma^{*}$
HOMO to a bonding p$\sigma$ LUMO metal-metal orbital, leading to Pt-Pt
contraction and subsequent coherent vibrations. Formation of an exciplex with
an additional metal ion, as in AgPtPOP115, 116, gives an opportunity to obtain
insights into metal-metal interactions and dynamics by excitation to the
lowest excited state. Computational studies of these dynamical processes
require the potential energy surface corresponding to the lowest excited
electronic state.

Figure 9 shows the calculated energy for the ground state and open-shell
singlet excited state corresponding to HOMO-LUMO excitation as a function of
the Ag-Pt distance. The Ag atom is placed along the Pt-Pt axis, and the
potential energy curves are generated using either a separate point
acquisition scheme where an excited state calculation for a given Ag-Pt
distance is initialized using the ground state orbitals at that distance, or
by a sequential scheme where the excited state calculation for an Ag-Pt
distance of 2.5 Å is initialized using the ground state orbitals, and then,
the converged excited state orbitals of the previous distance are used as the
initial guess for the next calculation, both increasing the Ag-Pt distance and
decreasing it. The DO-GMF calculations of the excited state are set to target
a 1st-order saddle point and they produce a smooth curve, with either
acquisition scheme, converging to the same solution for a given point on the
curve. However, the DO-MOM calculations with the separate acquisition scheme
converge on several different solutions as evidenced by the saddle point
orders ranging from 1 to larger than 25, resulting in a discontinuous energy
curve. Also, a few points on the curve could not be converged within 500
optimization steps. With the sequential point acquisition scheme, DO-MOM
converges to a consistent solution where the saddle point order is 3 with the
exception of the first point on the curve which corresponds to a 2nd-order
saddle point.

##  4 Discussion

The DO-GMF method presented here is found to be reliable and robust for
variational density functional calculations of excited electronic states by
converging on a saddle point of a specified order on the electronic energy
surface. The method has the distinct advantage that no precaution is needed to
prevent collapse to the ground state, such as MOM. It outperforms the DO-MOM
method which has, in turn, previously been shown to have superior convergence
properties than commonly used SCF-MOM approaches53, 49, 24.

Variational excited state calculations have a tendency to converge on the
stationary solution closest to the initial guess, irrespective of what saddle
point order it corresponds to. Therefore, if multiple electronic states emerge
as an atomic configuration is changed slightly, possibly breaking symmetries,
the variational calculation tends to converge on the state that is closest to
the initial guess, thereby conserving symmetry. This behavior can lead to
unphysical potential energy curves for atomic dynamics and give incorrect
excited states, as exemplified by the calculations of bond breaking in the H2
and ethylene molecules. DO-GMF is able to follow solutions that move away from
the original position on the electronic energy surface as the atomic
configuration changes, as the saddle point order of the converged solution is
conserved. This property is illustrated by the calculations of the symmetry-
broken doubly excited states of the H2 and ethylene molecules. The potential
energy curves calculated with DO-GMF are in qualitative agreement with
multireference results.

The DO-GMF calculation, however, requires as input the saddle point order of
the targeted excited states. The diagonal approximation to the Hessian using
the ground state orbitals with non-aufbau occupation numbers tends to
underestimate the saddle point order for excitations with large rearrangement
of the orbitals, as exemplified by the charge transfer excited states of
nitrobenzene and N-phenylpyrrole. Also, eigendecomposition of the full
electronic Hessian at the initial guess based on the ground state orbitals
does not provide the correct number of negative eigenvalues, as the electronic
energy surface at the initial guess does not have the same curvature as at the
converged excited state. In the cases of nitrobenzene and N-phenylpyrrole, the
number of negative eigenvalues obtained this way is much larger than the
saddle point order of the converged excited state. A constrained optimization
where the orbitals corresponding to the electron-hole pairs are fixed is found
to provide a better estimate of the saddle point order due to the partial
inclusion of relaxation effects. In the cases of nitrobenzene and
N-phenylpyrrole, the saddle point order estimated this way overestimates the
true saddle point order by only one, due to a negative eigenvalue that is
close to zero (larger than -1 eV). An improved strategy could be to re-
evaluate the diagonal approximation of the Hessian after the first few steps
of the wave function optimization, when the orbitals have started rearranging
in response to the new electron distribution. A localization transformation of
the initial orbitals may also be useful. Future efforts should aim to assess
and compare these possible strategies to estimate the saddle point order of
the target excited state.

The DO-MOM method would also benefit from a preconditioner that accurately
estimates the degrees of freedom along which the energy needs to be maximized.
If the relaxed solution is a saddle point of higher order than estimated by
the diagonal Hessian approximation at the initial guess, the DO-MOM
optimization can start collapsing to a lower-order saddle point before the
quasi-Newton algorithm is able to develop additional negative eigenvalues, as
seen in the calculations of the charge transfer excitation in N-phenylpyrrole.
There, MOM cannot avoid the variational collapse as the orbitals in the lower-
energy solution mix by less than 45∘. This failure shows that the commonly
used MOM method does not guarantee that variational collapse is prevented.
However, DO-GMF is less dependent on the preconditioner for the quasi-Newton
step, and variational collapse is excluded by construction, thus convergence
on a saddle point of the target order is guaranteed.

The generalized Davidson method used in DO-GMF requires an initial guess for
each target eigenvector. At the first wave function optimization step, the
most straightforward guess consists of the frontier orbital rotations up to
the target saddle point order in order of increasing corresponding element of
the diagonal Hessian approximation. However, the electronic Hessian
eigenvectors often contain significant contributions from multiple orbital
rotations. Therefore, such an initial guess may be far from the actual Hessian
eigenvectors. Subsequent partial Hessian diagonalizations typically converge
faster since the converged eigenvectors of the previous optimization step can
be used as the initial guess. In order to accelerate the first partial Hessian
diagonalization, it could be worthwhile to use a generalized version of the
dimer method to converge more than one of the lowest eigenpairs at the same
time. Such a generalized dimer method would optimize multiple orthogonal
dimers at the same time, meaning that multiple optimal dimer rotations are
found simultaneously. While the lack of a preconditioner for a generalized
dimer method could make it difficult to ensure convergence to the lowest
eigenvectors, the possibility of using the L-BFGS optimizer to find the
optimal dimer rotations could accelerate convergence significantly for the
first partial diagonalization step.

The computational effort of DO-GMF is greater than that of DO-MOM as DO-MOM
requires one energy/gradient evaluation per wave function optimization step,
while DO-GMF requires one additional energy/gradient calculation per target
eigenpair per step if a forward finite difference approximation is used in the
partial Hessian diagonalization. As the accuracy of the forward finite
difference approximation scales linearly in the finite difference step size,
the latter should be chosen as small as possible, while ensuring numerical
stability. A finite difference step size of $10^{-3}$ has been found to be
numerically stable and accurate enough for the calculations presented here.
Quadratic scaling can be obtained by using a central finite difference
approximation performing twice the amount of energy/gradient evaluations, but
we have not found it necessary in the calculations presented here.

##  5 Conclusion

The DO-GMF method is presented for variational calculations of excited
electronic states where a stationary solution is reached by converging on a
saddle point of a given order. The method involves a generalization of the
minimum mode following approach where an $n$th-order saddle point is found by
inverting the components of the gradient in the direction of the eigenvectors
of the $n$ lowest eigenvalues of the Hessian. The method recasts the saddle
point search as a minimization using the modified gradient, so variational
collapse to the ground state is prevented by construction. An implementation
is presented where the exponential transformation direct optimization is used
in combination with the L-BFGS algorithm to perform the orbital optimization
and the generalized Davidson method is used for the partial diagonalization of
the electronic Hessian. The performance of DO-GMF is demonstrated in
calculations using generalized gradient approximation density functionals of
excited states that are challenging for a previously proposed method based on
direct optimization and the maximum overlap method. DO-GMF is demonstrated to
track excited states through atomic configurations where the symmetry of the
wave function breaks in calculations of the potential energy curves of H2 and
ethylene by systematically targeting a saddle point of a given order. The
method has, furthermore, been shown to converge in a robust way in
calculations of challenging charge transfer excited states of nitrobenzene and
N-phenylpyrrole as well as the lowest excited state of a large diplatinum and
silver complex, AgPtPOP, as the Ag-Pt distance is varied.

The results presented here indicate that DO-GMF is a promising tool for Born-
Oppenheimer and non-adiabatic simulations of atomic dynamics involving
sequential excited state calculations. Additional tests on more complicated
systems are, of course, needed to verify that convergence on a saddle point of
a given order can be used as a guiding principle when symmetry-broken
solutions provide improved estimate of the energy of the system. Another issue
is that symmetry-broken solutions introduce undesirable features, such as
spurious dipole and magnetic moments, which can, in turn, lead to inaccurate
atomic forces along degrees of freedom other than bond breaking coordinates.
Quite likely, their symmetry-related counterparts will need to be taken into
account to accurately evaluate atomic forces over the entire energy surface.

This work was supported by the Icelandic Research Fund (grant agreements nos.
217751, 217734, 196070). The calculations were carried out at the Icelandic
High Performance Computing Center (IHPC). The authors thank Asmus O. Dohn,
Christoffer H. Egeberg and Kristoffer Haldrup for useful discussions and
support with the calculations on the AgPtPOP complex, and Elli Selenius for
support with the calculations on twisted N-phenylpyrrole.

The authors confirm that the data supporting the findings of this study are
available within the article and/or its supplementary materials.

Second Hessian eigenvalue along doubly excited state energy curve of H2;
Convergence of the $\pi_{\mathrm{ph}}^{*}\leftarrow\pi_{\mathrm{py}}$ charge
transfer excited state of twisted N-phenylpyrrole with DO-MOM and DO-GMF when
a smaller optimization step size of 0.1 is used; Optimized geometries of all
systems.

## References

  * Runge and Gross 1984 Runge, E.; Gross, E. K. U. Density-Functional Theory for Time-Dependent Systems. _Phys. Rev. Lett._ 1984, _52_ , 997
  * Casida 1995 Casida, M. E. _Recent Advances in Computational Chemistry_ ; 1995; pp 155–192
  * Herbert 2022 Herbert, J. M. Density Functional Theory for Electronic Excited States. 2022; https://arxiv.org/abs/2204.10135
  * Dreuw and Head-Gordon 2004 Dreuw, A.; Head-Gordon, M. Failure of Time-Dependent Density Functional Theory for Long-Range Charge-Transfer Excited States: The Zincbacteriochlorin-Bacteriochlorin and Bacteriochlorophyll-Spheroidene Complexes. _J. Am. Chem. Soc._ 2004, _126_ , 4007–4016
  * Dreuw and Head-Gordon 2005 Dreuw, A.; Head-Gordon, M. Single-reference ab initio methods for the calculation of excited states of large molecules. _Chem. Rev._ 2005, _105_ , 4009–4037
  * Hait and Head-Gordon 2021 Hait, D.; Head-Gordon, M. Orbital Optimized Density Functional Theory for Electronic Excited States. _J. Phys. Chem. Lett._ 2021, _12_ , 4517–4529
  * Hait et al. 2019 Hait, D.; Rettig, A.; Head-Gordon, M. Well-behaved versus ill-behaved density functionals for single bond dissociation: Separating success from disaster functional by functional for stretched H2. _J. Chem. Phys._ 2019, _150_
  * Hait et al. 2019 Hait, D.; Rettig, A.; Head-Gordon, M. Beyond the Coulson-Fischer point: Characterizing single excitation CI and TDDFT for excited states in single bond dissociations. _Phys. Chem. Chem. Phys._ 2019, _21_ , 21761–21775
  * Barca et al. 2018 Barca, G. M.; Gilbert, A. T.; Gill, P. M. Simple Models for Difficult Electronic Excitations. _J. Chem. Theory Comput._ 2018, _14_ , 1501–1509
  * Levine et al. 2006 Levine, B. G.; Ko, C.; Quenneville, J.; Martínez, T. J. Conical intersections and double excitations in time-dependent density functional theory. _Mol. Phys._ 2006, _104_ , 1039–1051
  * Athavale et al. 2021 Athavale, V.; Teh, H.-H.; Subotnik, J. On The Inclusion of One Double Within CIS and TD-DFT. _J. Chem. Phys._ 2021, _155_ , 154105
  * Teh and Subotnik 2019 Teh, H. H.; Subotnik, J. E. The Simplest Possible Approach for Simulating S0-S1 Conical Intersections with DFT/TDDFT: Adding One Doubly Excited Configuration. _J. Phys. Chem. Lett._ 2019, _10_ , 3426–3432
  * Shu et al. 2017 Shu, Y.; Parker, K. A.; Truhlar, D. G. Dual-Functional Tamm-Dancoff Approximation with Self-Interaction-Free Orbitals: Vertical Excitation Energies and Potential Energy Surfaces near an Intersection Seam. _J. Phys. Chem._ 2017, _121_ , 9728–9735
  * Shu et al. 2017 Shu, Y.; Parker, K. A.; Truhlar, D. G. Dual-Functional Tamm–Dancoff Approximation: A Convenient Density Functional Method that Correctly Describes S1/S0 Conical Intersections. _J. Phys. Chem. Lett._ 2017, _8_ , 2107–2112
  * Shao et al. 2003 Shao, Y.; Head-Gordon, M.; Krylov, A. I. The spin-flip approach within time-dependent density functional theory: Theory and applications to diradicals. _J. Chem. Phys._ 2003, _118_ , 4807–4818
  * Yang et al. 2016 Yang, Y.; Shen, L.; Zhang, D.; Yang, W. Conical Intersections from Particle-Particle Random Phase and Tamm-Dancoff Approximations. _J. Phys. Chem. Lett._ 2016, _7_ , 2407–2411
  * Bannwarth et al. 2020 Bannwarth, C.; Yu, J. K.; Hohenstein, E. G.; Martínez, T. J. Hole-hole Tamm-Dancoff-approximated density functional theory: A highly efficient electronic structure method incorporating dynamic and static correlation. _J. Chem. Phys._ 2020, _153_ , 024110
  * Ottochian et al. 2020 Ottochian, A.; Morgillo, C.; Ciofini, I.; Frisch, M. J.; Scalmani, G.; Adamo, C. Double hybrids and time-dependent density functional theory: An implementation and benchmark on charge transfer excited states. _J. Comput. Chem._ 2020, _41_ , 1242–1251
  * Brémond et al. 2021 Brémond, E.; Ottochian, A.; Pérez-Jiménez, A. J.; Ciofini, I.; Scalmani, G.; Frisch, M. J.; Sancho-García, J. C.; Adamo, C. Assessing challenging intra- and inter-molecular charge-transfer excitations energies with double-hybrid density functionals. _J. Comput. Chem._ 2021, _42_ , 970–981
  * Stein et al. 2009 Stein, T.; Kronik, L.; Baer, R. Reliable prediction of charge transfer excitations in molecular complexes using time-dependent density functional theory. _J. Am. Chem. Soc._ 2009, _131_ , 2818–2820
  * Kronik et al. 2012 Kronik, L.; Stein, T.; Refaely-Abramson, S.; Baer, R. Excitation gaps of finite-sized systems from optimally tuned range-separated hybrid functionals. _J. Chem. Theory and Comput._ 2012, _8_ , 1515–1531
  * Körzdörfer and Brédas 2014 Körzdörfer, T.; Brédas, J. L. Organic electronic materials: Recent advances in the DFT description of the ground and excited states using tuned range-separated hybrid functionals. _Acc. Chem. Res._ 2014, _47_ , 3284–3291
  * Vandaele et al. 2022 Vandaele, E.; Mališ, M.; Luber, S. The $\Delta$SCF method for non-adiabatic dynamics of systems in the liquid phase. _J. Chem. Phys_ 2022, _156_ , 130901
  * Levi et al. 2020 Levi, G.; Ivanov, A. V.; Jónsson, H. Variational Density Functional Calculations of Excited States via Direct Optimization. _J. Chem. Theory Comput._ 2020, _16_ , 6968–6982
  * Carter-Fenk and Herbert 2020 Carter-Fenk, K.; Herbert, J. M. State-Targeted Energy Projection: A Simple and Robust Approach to Orbital Relaxation of Non-Aufbau Self-Consistent Field Solutions. _J. Chem. Theory Comput._ 2020, _16_ , 5067–5082
  * Ayers et al. 2015 Ayers, P. W.; Levy, M.; Nagy, Communication: Kohn-Sham theory for excited states of Coulomb systems. _J. Chem. Phys._ 2015, _143_ , 191101
  * Zhekova et al. 2014 Zhekova, H. R.; Seth, M.; Ziegler, T. A perspective on the relative merits of time-dependent and time-independent density functional theory in studies of the electron spectra due to transition metal complexes. An illustration through applications to copper tetrachloride and plastocyanin. _Int. J. Quantum Chem._ 2014, _114_ , 1019–1029
  * Seidu et al. 2015 Seidu, I.; Krykunov, M.; Ziegler, T. Applications of time-dependent and time-independent density functional theory to Rydberg transitions. _J. Phys. Chem._ 2015, _119_ , 5107–5116
  * Cheng et al. 2008 Cheng, C. L.; Wu, Q.; Voorhis, T. V. Rydberg energies using excited state density functional theory. _J. Chem. Phys._ 2008, _129_ , 124112
  * Besley 2021 Besley, N. A. Modeling of the spectroscopy of core electrons with density functional theory. _Wiley Interdiscip. Rev. Comput. Mol. Sci._ 2021, _11_ , e1527
  * Besley et al. 2009 Besley, N. A.; Gilbert, A. T.; Gill, P. M. Self-consistent-field calculations of core excited states. _J. Chem. Phys._ 2009, _130_ , 124308
  * Hait and Head-Gordon 2020 Hait, D.; Head-Gordon, M. Excited State Orbital Optimization via Minimizing the Square of the Gradient: General Approach and Application to Singly and Doubly Excited States via Density Functional Theory. _J. Chem. Theory Comput._ 2020, _16_ , 1699–1710
  * Vandaele et al. 2022 Vandaele, E.; Mališ, M.; Luber, S. The photodissociation of solvated cyclopropanone and its hydrate explored via non-adiabatic molecular dynamics using $\Delta$SCF. _Phys. Chem. Chem. Phys._ 2022, _24_ , 5669–5679
  * Mališ and Luber 2020 Mališ, M.; Luber, S. Trajectory Surface Hopping Nonadiabatic Molecular Dynamics with Kohn-Sham $\Delta$SCF for Condensed-Phase Systems. _J. Chem. Theory Comput_ 2020, _16_ , 4071–4086
  * Pradhan et al. 2018 Pradhan, E.; Sato, K.; Akimov, A. V. Non-adiabatic molecular dynamics with $\Delta$SCF excited states. _J. Phys. Condens. Matter_ 2018, _30_ , 484001
  * Levi et al. 2018 Levi, G.; Pápai, M.; Henriksen, N. E.; Dohn, A. O.; Møller, K. B. Solution Structure and Ultrafast Vibrational Relaxation of the PtPOP Complex Revealed by $\Delta$SCF-QM/MM Direct Dynamics Simulations. _J. Phys. Chem. C_ 2018, _122_ , 7100–7119
  * Burton 2022 Burton, H. G. Energy Landscape of State-Specific Electronic Structure Theory. _J. Chem. Theory Comput._ 2022, _18_ , 1512–1526
  * Helgaker et al. 2014 Helgaker, T.; Jørgensen, P.; Olsen, J. _Molecular Electronic‐Structure Theory_ ; John Wiley & Sons, Ltd, 2014; Chapter 4, pp 107–141
  * Jensen and Jørgensen 1984 Jensen, H. J. A.; Jørgensen, P. A direct approach to second-order MCSCF calculations using a norm extended optimization scheme. _J. Chem. Phys._ 1984, _80_ , 1204–1214
  * Olsen et al. 1983 Olsen, J.; Yeager, D. L.; Jørgensen, P. Optimization and Characterization of a Multiconfigurational Self-Consistent Field (MCSCF) State. _Adv. Chem. Phys._ 1983, _54_ , 1–176
  * Golab et al. 1983 Golab, J. T.; Yeager, D. L.; Jørgensen, P. Proper characterization of MC SCF stationary points. _Chem. Phys._ 1983, _78_ , 175–199
  * Perdew and Levy 1985 Perdew, J. P.; Levy, M. Extrema of the density functional for the energy: Excited states from the ground-state theory. _Phys. Rev. B, Condens. Matter_ 1985, _31_ , 6264
  * Kowalczyk et al. 2011 Kowalczyk, T.; Yost, S. R.; Voorhis, T. V. Assessment of the $\Delta$SCF density functional theory approach for electronic excitations in organic dyes. _J. Chem. Phys._ 2011, _134_ , 054128
  * Davidson 1975 Davidson, E. R. The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real-symmetric matrices. _J. Comput. Phys._ 1975, _17_ , 87–94
  * Pulay 1980 Pulay, P. Convergence acceleration of iterative sequences. The case of SCF iteration. _Chem. Phys. Lett._ 1980, _73_ , 393–398
  * Pulay 1982 Pulay, P. Improved SCF Convergence Acceleration. _J. Comput. Chem._ 1982, _3_ , 556–560
  * Taka et al. 2022 Taka, A. A.; Lu, S.-Y.; Gowland, D.; Zuehlsdorff, T. J.; Corzo, H. H.; Pribram-Jones, A.; Shi, L.; Hratchian, H. P.; Isborn, C. M. Comparison of Linear Response Theory, Projected Initial Maximum Overlap Method, and Molecular Dynamics-Based Vibronic Spectra: The Case of Methylene Blue. _J. Chem. Theory Comput._ 2022, _18_ , 3039–3051
  * Gilbert et al. 2008 Gilbert, A. T.; Besley, N. A.; Gill, P. M. Self-consistent field calculations of excited states using the maximum overlap method (MOM). _J. Phys. Chem._ 2008, _112_ , 13164–13171
  * Ivanov et al. 2021 Ivanov, A. V.; Levi, G.; Jónsson, E. O.; Jónsson, H. Method for Calculating Excited Electronic States Using Density Functionals and Direct Orbital Optimization with Real Space Grid or Plane-Wave Basis Set. _J. Chem. Theory Comput._ 2021, _17_ , 5034–5049
  * Ivanov et al. 2021 Ivanov, A. V.; Jónsson, E. O.; Vegge, T.; Jónsson, H. Direct Energy Minimization Based on Exponential Transformation in Density Functional Calculations of Finite and Extended Systems. _Comput. Phys. Commun._ 2021, _267_ , 108047
  * Levi et al. 2020 Levi, G.; Ivanov, A. V.; Jónsson, H. Variational calculations of excited states: Via direct optimization of the orbitals in DFT. _Faraday Discuss._ 2020, _224_ , 448–466
  * Voorhis and Head-Gordon 2002 Voorhis, T. V.; Head-Gordon, M. A geometric approach to direct minimization. _Mol. Phys._ 2002, _100_ , 1713–1721
  * Schmerwitz et al. 2022 Schmerwitz, Y. L. A.; Ivanov, A. V.; Jónsson, E. O.; Jónsson, H.; Levi, G. Variational Density Functional Calculations of Excited States: Conical Intersection and Avoided Crossing in Ethylene Bond Twisting. _J. Phys. Chem. Lett._ 2022, _13_ , 3990–3999
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. _Phys. Rev. Lett._ 1996, _77_ , 3865
  * Perdew et al. 1997 Perdew, J. P.; Burke, K.; Ernzerhof, M. Erratum: Generalized Gradient Approximation Made Simple. _Phys. Rev. Lett._ 1997, _78_ , 1396
  * Jake et al. 2018 Jake, L. C.; Henderson, T. M.; Scuseria, G. E. Hartree-Fock symmetry breaking around conical intersections. _J. Chem. Phys._ 2018, _148_ , 024109
  * Tóth and Pulay 2016 Tóth, Z.; Pulay, P. Finding symmetry breaking Hartree-Fock solutions: The case of triplet instability. _J. Chem. Phys._ 2016, _145_ , 164102
  * Jiménez-Hoyos et al. 2011 Jiménez-Hoyos, C. A.; Henderson, T. M.; Scuseria, G. E. Generalized Hartree-Fock description of molecular dissociation. _J. Chem. Theory Comput._ 2011, _7_ , 2667–2674
  * Li and Paldus 2009 Li, X.; Paldus, J. Do independent-particle-model broken-symmetry solutions contain more physics than the symmetry-adapted ones? The case of homonuclear diatomics. _J. Chem. Phys._ 2009, _130_ , 084110
  * Coulson and Fischer 1949 Coulson, C. A.; Fischer, I. XXXIV. Notes on the molecular orbital treatment of the hydrogen molecule. _London Edinburgh Philos. Mag. J. Sci._ 1949, _40_ , 386–393
  * Perdew et al. 2021 Perdew, J. P.; Ruzsinszky, A.; Sun, J.; Nepal, N. K.; Kaplan, A. D. Interpretations of ground-state symmetry breaking and strong correlation in wavefunction and density functional theories. _Proc. Natl. Acad. Sci._ 2021, _118_ , e2017850118
  * Yu et al. 2016 Yu, H. S.; Li, S. L.; Truhlar, D. G. Perspective: Kohn-Sham density functional theory descending a staircase. _J. Chem. Phys._ 2016, _145_ , 130901
  * Cohen et al. 2008 Cohen, A. J.; Mori-Sánchez, P.; Yang, W. Insights into current limitations of density functional theory. _Science_ 2008, _321_ , 792–794
  * Gräfenstein et al. 2002 Gräfenstein, J.; Kraka, E.; Filatov, M.; Cremer, D. Can Unrestricted Density-Functional Theory Describe Open Shell Singlet Biradicals? _Int. J. Mol. Sci_ 2002, _3_ , 360–394
  * Cremer et al. 2002 Cremer, D.; Filatov, M.; Polo, V.; Kraka, E.; Shaik, S. Implicit and Explicit Coverage of Multi-reference Effects by Density Functional Theory. _Int. J. Mol. Sci_ 2002, _3_ , 604–638
  * Cremer 2001 Cremer, D. Density functional theory: Coverage of dynamic and non-dynamic electron correlation effects. _Mol. Phys._ 2001, _99_ , 1899–1940
  * Gräfenstein et al. 2000 Gräfenstein, J.; Hjerpe, A. M.; Kraka, E.; Cremer, D. An accurate description of the Bergman reaction using restricted and unrestricted DFT: Stability test, spin density, and on-top pair density. _J. Phys. Chem._ 2000, _104_ , 1748–1761
  * Wittbrodt and Schlegel 1996 Wittbrodt, J. M.; Schlegel, H. B. Some reasons not to use spin projected density functional theory. _J. Chem. Phys._ 1996, _105_ , 6574–6577
  * Mališ et al. 2022 Mališ, M.; Vandaele, E.; Luber, S. Spin-Orbit Couplings for Nonadiabatic Molecular Dynamics at the $\Delta$SCF Level. _J. Chem. Theory Comput._ 2022, _18_ , 4082–4094
  * Vaucher and Reiher 2017 Vaucher, A. C.; Reiher, M. Steering Orbital Optimization out of Local Minima and Saddle Points Toward Lower Energy. _J. Chem. Theory Comput._ 2017, _13_ , 1219–1228
  * Pelzer and Wigner 1932 Pelzer, H.; Wigner, E. P. Über die Geschwindigkeitskonstante von Austauschreaktionen. _Z. Phys. Chem. B_ 1932, _15_ , 445–471
  * Eyring 1935 Eyring, H. The activated complex in chemical reactions. _J. Chem. Phys._ 1935, _3_ , 63–71
  * Wigner 1938 Wigner, E. The transition state method. _Trans. Faraday Soc._ 1938, _34_ , 29–41
  * Cerjan and Miller 1981 Cerjan, C. J.; Miller, W. H. On finding transition states. _J. Chem. Phys._ 1981, _75_ , 2800–2801
  * Simons et al. 1983 Simons, J.; Jorgensen, P.; Taylor, H.; Orment, J. Walking on Potential Energy Surfaces. _J. Phys. Chem_ 1983, _87_ , 2745–2753
  * Banerjee et al. 1985 Banerjee, A.; Adams, N.; Simons, J.; Shepard, R. Search for Stationary Points on Surfaces. _J. Phys. Chem._ 1985, _89_ , 52–57
  * Marie and Burton 2023 Marie, A.; Burton, H. G. A. Excited states, symmetry breaking, and unphysical solutions in state-specific CASSCF theory. 2023; https://arxiv.org/abs/2301.11731
  * Hoffmann et al. 2002 Hoffmann, M. R.; Sherrill, C. D.; Leininger, M. L.; Schaefer, H. F. Optimization of MCSCF excited states using directions of negative curvature. _Chem. Phys. Lett._ 2002, _355_ , 183–192
  * Murtagh and Sargent 1970 Murtagh, B. A.; Sargent, R. W. H. Computational experience with quadratically convergent minimisation methods. _Comput. J._ 1970, _13_ , 185–194
  * Powell 1973 Powell, M. J. D. A New Algorithm for Unconstrained Optimization. _Nonlinear Programming_ 1973, 31–65
  * Bofill 1994 Bofill, J. M. Updated Hessian Matrix and the Restricted Step Method for Locating Transition Structures. _J. Comput. Chem._ 1994, _15_ , 1–11
  * Henkelman and Jónsson 1999 Henkelman, G.; Jónsson, H. A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives. _J. Chem. Phys._ 1999, _111_ , 7010–7022
  * Olsen et al. 2004 Olsen, R. A.; Kroes, G. J.; Henkelman, G.; Arnaldsson, A.; Jónsson, H. Comparison of methods for finding saddle points without knowledge of the final states. _J. Chem. Phys._ 2004, _121_ , 9776–9792
  * Kästner and Sherwood 2008 Kästner, J.; Sherwood, P. Superlinearly converging dimer method for transition state search. _J. Chem. Phys._ 2008, _128_ , 014106
  * Manuel Plasencia Gutiérrez and Jónsson 2016 Manuel Plasencia Gutiérrez, C. A.; Jónsson, H. Improved Minimum Mode Following Method for Finding First Order Saddle Points. _J. Chem. Theory Comput._ 2016, _13_ , 125–134
  * Lánczos 1950 Lánczos, C. An Iteration Method for the Solution of the Eigenvalue Problem of Linear Differential and Integral Operators. _J. Res. Natl. Bur. Stand._ 1950, _45_ , 255–282
  * Crouzeix et al. 1994 Crouzeix, M.; Philippe, B.; Sadkane, M. The Davidson Method. _SIAM J. Sci. Comput._ 1994, _15_ , 62–76
  * Lehtola et al. 2020 Lehtola, S.; Blockhuys, F.; Van Alsenoy, C. An overview of self-consistent field calculations within finite basis sets. _Molecules_ 2020, _25_ , 1–23
  * Perdew and Zunger 1981 Perdew, J. P.; Zunger, A. Self-interaction correction to density-functional approximations for many-electron systems. _Phys. Rev. B, Condens. Matter_ 1981, _23_ , 5048–5079
  * Lehtola et al. 2016 Lehtola, S.; Head-Gordon, M.; Jónsson, H. Complex orbitals, multiple local minima and symmetry breaking in Perdew-Zunger self-interaction corrected density-functional theory calculations. _J. Chem. Theory Comput._ 2016, _12_ , 3195
  * Baker 1905 Baker, H. F. Alternants and Continuous Groups. _Proc. London Math. Soc._ 1905, _2_ , 24–47
  * Campbell 1896 Campbell, J. E. On a Law of Combination of Operators bearing on the Theory of Continuous Transformation Groups. _Proc. London Math. Soc._ 1896, _s1-28_ , 381–390
  * Campbell 1897 Campbell, J. E. On a Law of Combination of Operators (Second Paper). _Proc. London Math. Soc._ 1897, _s1-29_ , 14–32
  * Hausdorff 1906 Hausdorff, F. Die symbolische Exponentialformel in der Gruppentheorie. _Ber. Verh. Kgl. Sächs. Ges. Wiss. Leipzig., Math.-phys. Kl._ 1906, _58_ , 19–48
  * Salem and Rowland 1972 Salem, L.; Rowland, C. The electronic properties of diradicals. _Angew. Chem. Int. Ed._ 1972, _11_ , 92–111
  * Gould et al. 2014 Gould, N.; Ortner, C.; Packwood, D. An Efficient Dimer Method With Preconditioning And Linesearch. 2014; http://arxiv.org/abs/1407.2817
  * Sharada et al. 2015 Sharada, S. M.; Stück, D.; Sundstrom, E. J.; Bell, A. T.; Head-Gordon, M. Wavefunction stability analysis without analytical electronic Hessians: Application to orbital-optimised second-order Møller-Plesset theory and VV10-containing density functionals. _Mol. Phys._ 2015, _113_ , 1802–1808
  * Mortensen et al. 2005 Mortensen, J. J.; Hansen, L. B.; Jacobsen, K. W. Real-space grid implementation of the projector augmented wave method. _Phys. Rev. B, Condens. Matter_ 2005, _71_ , 035109
  * Enkovaara et al. 2010 Enkovaara, J.; Rostgaard, C.; Mortensen, J. J.; Chen, J.; Dułak, M.; Ferrighi, L.; Gavnholt, J.; Glinsvad, C.; Haikola, V.; Hansen, H. A.; Kristoffersen, H. H.; Kuisma, M.; Larsen, A. H.; Lehtovaara, L.; Ljungberg, M.; Lopez-Acevedo, O.; Moses, P. G.; Ojanen, J.; Olsen, T.; Petzold, V.; Romero, N. A.; Stausholm-Møller, J.; Strange, M.; Tritsaris, G. A.; Vanin, M.; Walter, M.; Hammer, B.; Häkkinen, H.; Madsen, G. K.; Nieminen, R. M.; Nørskov, J. K.; Puska, M.; Rantala, T. T.; Schiøtz, J.; Thygesen, K. S.; Jacobsen, K. W. Electronic structure calculations with GPAW: A real-space implementation of the projector augmented-wave method. _J. Phys. Condens. Matter_ 2010, _22_ , 253202
  * Larsen et al. 2009 Larsen, A. H.; Vanin, M.; Mortensen, J. J.; Thygesen, K. S.; Jacobsen, K. W. Localized atomic basis set in the projector augmented wave method. _Phys. Rev. B, Condens. Matter_ 2009, _80_ , 195112
  * Dunning 1989 Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. _J. Chem. Phys._ 1989, _90_ , 1007–1023
  * Kendall et al. 1992 Kendall, R. A.; Dunning, T. H.; Harrison, R. J. Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions. _J. Chem. Phys._ 1992, _96_ , 6796–6806
  * Woon and Dunning 1994 Woon, D. E.; Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. IV. Calculation of static electrical response properties. _J. Chem. Phys._ 1994, _100_ , 2975–2988
  * Weigend and Ahlrichs 2005 Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. _Phys. Chem. Chem. Phys._ 2005, _7_ , 3297–3305
  * Rossi et al. 2015 Rossi, T. P.; Lehtola, S.; Sakko, A.; Puska, M. J.; Nieminen, R. M. Nanoplasmonics simulations at the basis set limit through completeness-optimized, local numerical basis sets. _J. Chem. Phys._ 2015, _142_ , 094114
  * Blöchl 1994 Blöchl, P. E. Projector augmented-wave method. _Phys. Rev. B, Condens. Matter_ 1994, _50_ , 17953
  * Loos et al. 2021 Loos, P. F.; Comin, M.; Blase, X.; Jacquemin, D. Reference Energies for Intramolecular Charge-Transfer Excitations. _J. Chem. Theory Comput._ 2021, _17_ , 3666–3686
  * Lehtola et al. 2018 Lehtola, S.; Steigemann, C.; Oliveira, M. J.; Marques, M. A. Recent developments in LIBXC — A comprehensive library of functionals for density functional theory. _SoftwareX_ 2018, _7_ , 1–5
  * Barbatti and Crespo-Otero 2016 Barbatti, M.; Crespo-Otero, R. In _Density-Functional Methods for Excited States_ ; Ferré, N., Filatov, M., Huix-Rotllant, M., Eds.; Springer International Publishing: Cham, 2016; pp 415–444
  * Mewes et al. 2014 Mewes, J. M.; Jovanović, V.; Marian, C. M.; Dreuw, A. On the molecular mechanism of non-radiative decay of nitrobenzene and the unforeseen challenges this simple molecule holds for electronic structure theory. _Phys. Chem. Chem. Phys._ 2014, _16_ , 12393–12406
  * Gray et al. 2017 Gray, H. B.; Záliš, S.; Vlček, A. Electronic structures and photophysics of d8-d8 complexes. _Coord. Chem. Rev._ 2017, _345_ , 297–317
  * Haldrup et al. 2019 Haldrup, K.; Levi, G.; Biasin, E.; Vester, P.; Laursen, M. G.; Beyer, F.; Kjær, K. S.; Brandt Van Driel, T.; Harlang, T.; Dohn, A. O.; Hartsock, R. J.; Nelson, S.; Glownia, J. M.; Lemke, H. T.; Christensen, M.; Gaffney, K. J.; Henriksen, N. E.; Møller, K. B.; Nielsen, M. M. Ultrafast X-Ray Scattering Measurements of Coherent Structural Dynamics on the Ground-State Potential Energy Surface of a Diplatinum Molecule. _Phys. Rev. Lett._ 2019, _122_ , 63001
  * Monni et al. 2018 Monni, R.; Capano, G.; Auböck, G.; Gray, H. B.; Vlček, A.; Chergui, M. Vibrational coherence transfer in the ultrafast intersystem crossing of a diplatinum complex in solution. _Proc. Natl. Acad. Sci._ 2018, _115_ , E6396–E6403
  * Van Der Veen et al. 2011 Van Der Veen, R. M.; Cannizzo, A.; Mourik, F. V.; Vlček, A.; Chergui, M. Vibrational Relaxation and Intersystem Crossing of Binuclear Metal Complexes in Solution. _J. Am. Chem. Soc._ 2011, _133_ , 305–315
  * Clodfelter et al. 1994 Clodfelter, S. A.; Doede, T. M.; Brennan, B. A.; Nagle, J. K.; Bender, D. P.; Turner, W. A.; Lapunzina, P. M. Luminescent Metal-Metal-Bonded Exciplexes Involving Tetrakis@-diphosphito)diplatinate(II) and Thallium(1). _J. Am. Chem. Soc._ 1994, _116_ , 11379–11386
  * Christensen et al. 2010 Christensen, M.; Haldrup, K.; Kjær, K. S.; Cammarata, M.; Wulff, M.; Bechgaard, K.; Weihe, H.; Harrit, N. H.; Nielsen, M. M. Structure of a short-lived excited state trinuclear Ag-Pt-Pt complex in aqueous solution by time resolved X-ray scattering. _Phys. Chem. Chem. Phys._ 2010, _12_ , 6921–6923