9.8.1 Nonadiabatic Couplings

Conical intersections are degeneracies between Born-Oppenheimer potential energy surfaces that facilitate nonadiabatic transitions between excited states, i.e., internal conversion and intersystem crossing processes, both of which represent a breakdown of the Born-Oppenheimer approximation.^{758, 459, Herbert:2022} Although simultaneous intersections between more than two electronic states are possible,⁷⁵⁸ consider for convenience the two-state case, and let

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 $𝐑$ represents the remaining, nuclear coordinates.) By definition, the Born-Oppenheimer states are the ones that diagonalize $𝐇$ at a particular molecular geometry $𝐑$, 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}}=3N_{\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 ${𝐠}_{JK}$ and ${𝐡}_{JK}$. The former is simply the difference in the gradient vectors of the two states in question,

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 nonadiabatic coupling vector ${𝐡}_{JK}$, on the other hand, is more involved and not directly amenable to finite-difference calculations:

This is closely related to the derivative coupling vector

The latter expression for ${𝐝}_{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 ${𝐠}_{JK}$ and ${𝐡}_{JK}$ determined the topography around the intersection, i.e., whether the intersection is “peaked” or “sloped”;^{45, 1256} see Ref. 459 for a pedagogical overview.

Algorithms to compute the nonadiabatic couplings ${𝐝}_{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,^{314, 1271, 1272, 847} or at the corresponding spin-flip (SF) levels of theory (SF-CIS or SF-TDDFT). The spin-flip implementation¹²⁷¹ 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.^Herbert:2022

To understand why, suppose that $J$ in Eq. (9.5) 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}(𝐑)\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.⁶⁷² (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.^{1271, 459, Herbert:2022} This can have consequences for simulation of internal conversion to the ground state, where for example S${}_1$ lifetimes may be significantly affected by warping of the potential energy surfaces around the S${}_1$/S${}_0$ intersection in conventional TDDFT.¹²⁷⁴

Nonadiabatic (derivative) couplings are available for both CIS and TDDFT. The CIS nonadiabatic couplings can be obtained from direct differentiation of the wave functions with respect to nuclear positions.^{314, 1271} For TDDFT, the same procedure can be carried out to calculate the approximate nonadiabatic couplings, in what has been termed the “pseudo-wave function” approach.^{848, 1271} Formally more rigorous TDDFT nonadiabatic 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.7), 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$.^{1272, 847} As such, the pseudo-wave function method is the recommended approach for computing nonadiabatic 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.¹²⁷²

Finally, we note that there is mounting evidence that SF-TDDFT calculations are most accurate when used with functionals containing $\sim$50% Hartree-Fock exchange,^{1022, 493, Herbert:2022} and many studies with this method have used the BH&HLYP functional; see Refs. 459 and Herbert:2022 for reviews. The BH&HLYP functional combines 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.)