10.6 Nonadiabatic Couplings and Optimization of Minimum-Energy Crossing Points

10.6.3 Minimum-Energy Crossing Points

The seam space of a conical intersection is really a (hyper)surface of dimension Nint-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 (Nint-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. 1037 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 10.7.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. (10.3)] and 𝐡JK [Eq. (10.4)]. 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 coworkers557 developed a penalty-constrained MECP optimization algorithm that consists of minimizing the objective function

Fσ(𝐑)=12[EI(𝐑)+EJ(𝐑)]+σ([EI(𝐑)-EJ(𝐑)]2EI(𝐑)-EJ(𝐑)+α), (10.6)

where α is a fixed parameter to avoid singularities and σ is a Lagrange multiplier for a penalty function meant to drive the energy gap to zero. Minimization of Fσ is performed iteratively for increasingly large values σ.

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

𝐆=𝐏𝐆mean+2(EK-EJ)𝐆diff, (10.7)

where

𝐆mean=12(𝐆J+𝐆K) (10.8)

is the mean energy gradient, with 𝐆i=Ei/𝐑 being the nuclear gradient for state i, and

𝐆diff=𝐆K-𝐆J||𝐆K-𝐆J|| (10.9)

is the normalized difference gradient. Finally,

𝐏=𝟏-𝐆diff𝐆diff (10.10)

projects the gradient difference direction out of the mean energy gradient in Eq. (10.7). 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 coworkers595 and implemented in Q-Chem by Zhang and Herbert.1037 This algorithm uses a gradient that is similar to that in Eq. (10.7) but projects out not just 𝐆diff in Eq. (10.10) 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 10.6.1. With spin-flip approaches, correct topology is obtained.1038

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,60 is available in Q-Chem and consists of optimization along the gradient in Eq. (10.7) but with a projector

𝐏=𝟏-𝐆diff𝐆diff-𝐲𝐲 (10.11)

where

𝐲=(𝟏-𝐱𝐱)𝐡JK||(𝟏-𝐱𝐱)𝐡JK||, (10.12)

in place of the projector in Eq. (10.10). Equation (10.11) has the effect of projecting the span of 𝐠JK and 𝐡JK (i.e., the branching space) out of state-averaged gradient in Eq. (10.7). 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.1038 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.1040 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 S0, S1, S,2 but will also generate the MS=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 S^2, but severe spin contamination can sometimes defeat this algorithm,1037 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 10.6.5.