9.7.3 Minimum-Energy Crossing Points

The seam space of a conical intersection is really a (hyper)surface of dimension $N_{\mathrm{int}}-2$, and while the two electronic states in question are degenerate at every point within this space, the electronic energy varies from one point to the next. To provide a simple picture of photochemical reaction pathways, it is often convenient to locate the minimum-energy crossing point (MECP) within this $(N_{\mathrm{int}}-2)$-dimensional seam. Two separate minimum-energy pathway searches, one on the excited state starting from the ground-state geometry and terminating at the MECP, and the other on the ground state starting from the MECP and terminating at the ground-state geometry, then affords a photochemical mechanism. (See Ref. 1099 for a simple example.) In some sense, then, the MECP is to photochemistry what the transition state is to reactions that occur on a single Born-Oppenheimer potential energy surface. One should be wary of pushing this analogy too far, because whereas a transition state reasonably be considered to be a bottleneck point on the reaction pathway, the path through a conical intersection may be downhill and perhaps therefore more likely to proceed from one surface to the other at a point “near" the intersection, and in addition there can be multiple conical intersections between the same pair of states so more than one photochemical mechanism may be at play. Such complexity could be explored, albeit at significantly increased cost, using non-adiabatic “surface hopping" ab initio molecular dynamics, as described in Section 9.8.6. Here we describe the computationally-simpler procedure of locating an MECP along a conical seam.

Recall that the branching space around a conical intersection between electronic states $J$ and $K$ is spanned by two vectors, ${𝐠}_{JK}$ [Eq. (9.4)] and ${𝐡}_{JK}$ [Eq. (9.5)]. While the former is readily available in analytic form for any electronic structure method that has analytic excited-state gradients, the non-adiabatic coupling vector ${𝐡}_{JK}$ is not available for most methods. For this reason, several algorithms have been developed to optimize MECPs without the need to evaluate ${𝐡}_{JK}$, and three such algorithms are available in Q-Chem.

Martínez and coworkers⁵⁷⁰ developed a penalty-constrained MECP optimization algorithm that consists of minimizing the objective function

$$F_{\sigma }(𝐑)=\frac{1}{2}\left[E_I(𝐑)+E_J(𝐑)\right]+\sigma \left(\frac{{\left[E_I(𝐑)-E_J(𝐑)\right]}^2}{E_I(𝐑)-E_J(𝐑)+\alpha }\right),$$

(9.7)

where $\alpha$ is a fixed parameter to avoid singularities and $\sigma$ is a Lagrange multiplier for a penalty function meant to drive the energy gap to zero. Minimization of $F_{\sigma }$ is performed iteratively for increasingly large values $\sigma$.

A second MECP optimization algorithm is a simplification of the penalty-constrained approach that we call the “direct” method. Here, the gradient of the objective function is

$$𝐆={\mathrm{𝐏𝐆}}_{\mathrm{mean}}+2(E_K-E_J){𝐆}_{\mathrm{diff}},$$

(9.8)

where

$${𝐆}_{\mathrm{mean}}=\frac{1}{2}({𝐆}_J+{𝐆}_K)$$

(9.9)

is the mean energy gradient, with ${𝐆}_i=\partial E_i/\partial 𝐑$ being the nuclear gradient for state $i$, and

$${𝐆}_{\mathrm{diff}}=\frac{{𝐆}_K-{𝐆}_J}{||{𝐆}_K-{𝐆}_J||}$$

(9.10)

is the normalized difference gradient. Finally,

$$𝐏=\mathrm{𝟏}-{𝐆}_{\mathrm{diff}}{𝐆}_{\mathrm{diff}}^{\top }$$

(9.11)

projects the gradient difference direction out of the mean energy gradient in Eq. (9.8). The algorithm then consists in minimizing along the gradient $𝐆$, with for the iterative cycle over a Lagrange multiplier, which can sometimes be slow to converge.

The third and final MECP optimization algorithm that is available in Q-Chem is the branching-plane updating method developed by Morokuma and coworkers⁶¹² and implemented in Q-Chem by Zhang and Herbert.¹⁰⁹⁹ This algorithm uses a gradient that is similar to that in Eq. (9.8) but projects out not just ${𝐆}_{\mathrm{diff}}$ in Eq. (9.11) but also a second vector that is orthogonal to it, representing an iteratively-updated approximation to the branching space.

None of these three methods requires evaluation of non-adiabatic couplings, and all three can be used to optimize MECPs at the CIS, SF-CIS, TDDFT, SF-TDDFT, and SOS-CIS(D0) levels. The direct algorithm can also be used for EOM-XX-CCSD methods (XX = EE, IP, or EA). It should be noted that since EOM-XX-CCSD is a linear response method, it suffers from the same topology problem around conical intersections involving the ground state that was described in regards to TDDFT in Section 9.7.1. With spin-flip approaches, correct topology is obtained.¹¹⁰⁰

Analytic derivative couplings are available for (SF-)CIS and (SF-)TDDFT, so for these methods one can alternatively employ an optimization algorithm that makes use of both ${𝐠}_{JK}$ and ${𝐡}_{JK}$. Such an algorithm, due to Schlegel and coworkers,⁴⁵ is available in Q-Chem and consists of optimization along the gradient in Eq. (9.8) but with a projector

$$𝐏=\mathrm{𝟏}-{𝐆}_{\mathrm{diff}}{𝐆}_{\mathrm{diff}}^{\top }-{\mathrm{𝐲𝐲}}^{\top }$$

(9.12)

where

$$𝐲=\frac{(\mathrm{𝟏}-{\mathrm{𝐱𝐱}}^{\top }){𝐡}_{JK}}{||(\mathrm{𝟏}-{\mathrm{𝐱𝐱}}^{\top }){𝐡}_{JK}||},$$

(9.13)

in place of the projector in Eq. (9.11). Equation (9.12) has the effect of projecting the span of ${𝐠}_{JK}$ and ${𝐡}_{JK}$ (i.e., the branching space) out of state-averaged gradient in Eq. (9.8). The tends to reduce the number of iterations necessary to converge the MECP, and since calculation of the (optional) ${𝐡}_{JK}$ vector represents only a slight amount of overhead on top of the (required) ${𝐠}_{JK}$ vector, this last algorithm tends to yield significant speed-ups relative to the other three.¹¹⁰⁰ As such, it is the recommended choice for (SF-)CIS and (SF-)TDDFT.

It should be noted that while the spin-flip methods cure the topology problem around conical intersections that involve the ground state, this method tends to exacerbate spin contamination relative to the corresponding spin-conserving approaches.¹¹⁰² While spin contamination is certainly present in traditional, spin-conserving CIS and TDDFT, it presents the following unique challenge in spin-flip methods. Suppose, for definiteness, that one is interested in singlet excited states. Then the reference state for the spin-flip methods should be the high-spin triplet. A spin-flipping excitation will then generate S${}_0$, S${}_1$, S${}_2,\mathrm{\dots }$ but will also generate the $M_S=0$ component of the triplet reference state, which therefore appears in what is ostensibly the singlet manifold. Q-Chem attempts to identify this automatically, based on a threshold for $⟨{\widehat{S}}^2⟩$, but severe spin contamination can sometimes defeat this algorithm,¹⁰⁹⁹ hampering Q-Chem’s ability to distinguish singlets from triplets (in this particular example). An alternative might be the state-tracking procedure that is described in Section 9.7.5.