Conical intersections are degeneracies between Born-Oppenheimer potential
energy surfaces that facilitate non-adiabatic transitions between excited
states, *i.e.*, internal conversion and intersystem crossing processes, both of
which represent a breakdown of the Born-Oppenheimer
approximation.^{Matsika:2011, Herbert:2016b} Although simultaneous
intersections between more than two electronic states are
possible,^{Matsika:2011} consider for convenience the two-state case, and let

$$\mathbf{H}=\left(\begin{array}{cc}\hfill {H}_{JJ}(\mathbf{R})\hfill & \hfill {H}_{JK}(\mathbf{R})\hfill \\ \hfill {H}_{JK}^{\ast}(\mathbf{R})\hfill & \hfill {H}_{KK}(\mathbf{R})\hfill \end{array}\right).$$ | (9.3) |

denote the matrix representation of the vibronic (vibrational + electronic)
Hamiltonian [Eq. (4.2)] in a basis of two electronic states,
$J$ and $K$. (Electronic degrees of freedom have been integrated out of this
expression, and $\mathbf{R}$ represents the remaining, nuclear coordinates.)
By definition, the Born-Oppenheimer states are the ones that diagonalize
$\mathbf{H}$ at a particular molecular geometry $\mathbf{R}$, and thus two
conditions must be satisfied in order to obtain degeneracy in the
Born-Oppenheimer representation: ${H}_{JJ}={H}_{KK}$ and ${H}_{JK}=0$. As
such, degeneracies between two Born-Oppenheimer potential energy surfaces exist
in subspaces of dimension ${N}_{\mathrm{int}}-2$, where ${N}_{\mathrm{int}}=3{N}_{\mathrm{atoms}}-6$ is the number of internal (vibrational) degrees of freedom (assuming
the molecule is non-linear). This (${N}_{\mathrm{int}}-2$)-dimensional subspace is
known as the *seam space* because the two states are degenerate everywhere
within this space. In the remaining two degrees of freedom, known as the
*branching space*, the degeneracy between Born-Oppenheimer surfaces is
lifted by an infinitesimal displacement, which in a three-dimensional plot
resembles a double cone about the point of intersection, hence the name
*conical intersection*.

The branching space is defined by the span of a pair of vectors ${\mathbf{g}}_{JK}$ and ${\mathbf{h}}_{JK}$. The former is simply the difference in the gradient vectors of the two states in question,

$${\mathbf{g}}_{JK}=\frac{\partial {E}_{J}}{\partial \text{\mathbf{R}}}-\frac{\partial {E}_{K}}{\partial \text{\mathbf{R}}},$$ | (9.4) |

and is readily evaluated at any level of theory for which analytic energy
gradients are available (or less-readily, via finite difference, if they are
not!). The definition of the *non-adiabatic coupling vector*
${\mathbf{h}}_{JK}$, on the other hand, is more involved and not directly
amenable to finite-difference calculations:

$${\mathbf{h}}_{JK}=\u27e8{\mathrm{\Psi}}_{J}\left|\left(\partial \widehat{H}/\partial \mathbf{R}\right)\right|{\mathrm{\Psi}}_{K}\u27e9.$$ | (9.5) |

This is closely related to the *derivative coupling* vector

$${\mathbf{d}}_{JK}=\u27e8{\mathrm{\Psi}}_{J}\left|\left(\partial /\partial \mathbf{R}\right)\right|{\mathrm{\Psi}}_{K}\u27e9=\frac{{\mathbf{h}}_{JK}}{{E}_{J}-{E}_{K}}.$$ | (9.6) |

The latter expression for ${\mathbf{d}}_{JK}$ demonstrates that the coupling
between states becomes large in regions of the potential surface where the two
states are nearly degenerate. The relative orientation and magnitudes of the
vectors ${\mathbf{g}}_{JK}$ and ${\mathbf{h}}_{JK}$ determined the topography
around the intersection, *i.e.*, whether the intersection is “peaked” or
“sloped”; see Ref. Herbert:2016b for a pedagogical overview.

Algorithms to compute the non-adiabatic couplings ${\mathbf{d}}_{JK}$ are not
widely available in quantum chemistry codes, but thanks to the efforts of the
Herbert and Subotnik groups, they are available in Q-Chem when the wave
functions ${\mathrm{\Psi}}_{J}$ and ${\mathrm{\Psi}}_{K}$, and corresponding electronic energies ${E}_{J}$
and ${E}_{K}$, are computed at the CIS or TDDFT
level,^{Fatehi:2011, Zhang:2014b, Zhang:2015a, Ou:2015} or at the
corresponding spin-flip (SF) levels of theory (SF-CIS or SF-TDDFT). The
spin-flip implementation^{Zhang:2014b} is particularly significant,
because only that approach—and not traditional spin-conserving CIS or
TDDFT—affords correct topology around conical intersections that involve the
ground state. To understand why, suppose that $J$ in Eq. (9.3)
represents the ground state; call it $J=0$ for definiteness. In linear
response theory (TDDFT) or in CIS (by virtue of Brillouin’s theorem), the
coupling matrix elements between the reference (ground) state and all of the
excited states vanish identically, hence ${H}_{0K}(\mathbf{R})\equiv 0$. This
means that there is only one condition to satisfy in order to obtain
degeneracy, hence the branching space is one- rather than two-dimensional, for
any conical intersection that involves the ground state.^{Levine:2006}
(For intersections between two excited states, the topology should be correct.)
In the spin-flip approach, however, the reference state has a different spin
multiplicity than the target states; if the latter have spin quantum number
$S$, then the reference state has spin $S+1$. This has the effect that the
ground state of interest (spin $S$) is treated as an excitation, and thus on a
more equal footing with excited states of the same spin, and it rigorously
fixes the topology problem around conical intersections.^{Zhang:2014b}

Nonadiabatic (derivative) couplings are available for both CIS and TDDFT. The
CIS non-adiabatic couplings can be obtained from direct differentiation of the
wave functions with respect to nuclear
positions.^{Fatehi:2011, Zhang:2014b} For TDDFT, the same procedure can be
carried out to calculate the approximate non-adiabatic couplings, in what has
been termed the “pseudo-wave function” approach.^{Ou:2014, Zhang:2014b}
Formally more rigorous TDDFT non-adiabatic couplings derived from quadratic
response theory are also available, although they are subject to certain
undesirable, accidental singularities if for the two states $J$ and $K$ in
Eq. (9.5), the energy difference $|{E}_{J}-{E}_{K}|$ is quasi-degenerate
with the excitation energy ${\omega}_{I}={E}_{I}-{E}_{0}$ for some third state,
$I$.^{Zhang:2015a, Ou:2015} As such, the pseudo-wave function method is
the recommended approach for computing non-adiabatic couplings with TDDFT,
although in the spin flip case the pseudo-wave function approach is rigorously
equivalent to the pseudo-wave function approach, and is free of singularities.^{Zhang:2015a}

Finally, we note that there is some evidence that SF-TDDFT calculations are
most accurate when used with functionals containing $\sim $50% Hartree-Fock
exchange,^{Shao:2003, Huix-Rotllant:2010} and many studies with this method
(see Ref. Herbert:2016b for a survey) have used the BH&HLYP
functional, in which LYP correlation is combined with Becke’s “half and half”
(BH&H) exchange functional, consisting of 50% Hartree-Fock exchange and
50% Becke88 exchange (EXCHANGE = BHHLYP in Q-Chem.)