Comparison of Spin-Flip TDDFT-Based Conical Intersection Approaches with XMS-CASPT2

Determining conical intersection geometries is of key importance to understanding the photochemical reactivity of molecules. While many small- to medium-sized molecules can be treated accurately using multireference approaches, larger molecules require a less computationally demanding approach. In this work, minimum energy crossing point conical intersection geometries for a series of molecules have been studied using spin-flip TDDFT (SF-TDDFT), within the Tamm-Dancoff Approximation, both with and without explicit calculation of nonadiabatic coupling terms, and compared with both XMS-CASPT2 and CASSCF calculated geometries. The less computationally demanding algorithms, which do not require explicit calculation of the nonadiabatic coupling terms, generally fare well with the XMS-CASPT2 reference structures, while the relative energetics are only reasonably replicated with the MECP structure as calculated with the BHHLYP functional and full nonadiabatic coupling terms. We also demonstrate that, occasionally, CASSCF structures deviate quantitatively from the XMS-CASPT2 structures, showing the importance of including dynamical correlation.


INTRODUCTION
Photochemical processes are ubiquitous in nature and form the backbone of many important chemical processes. In nature, the photoinduced isomerization of retinal forms the basis of vision, while absorption of light by chlorophyll is important in photosynthesis. In terms of man-made processes, dyesensitized solar cells, 1 fluorescent molecular probes (a full literature review is beyond the scope of the current work; see, e.g., ref 2), chemosensors, 3 energy transfer cassettes, 4 photodynamic therapy agents, 5 and tunable laser dyes 6−8 are just a handful of the many applications of photochemistry. Key to the correct computational description of absorption of radiation by such molecules is the transition dipole moment. 9 Upon absorption, the molecule may relax to the ground electronic state via different routes: emission, phosphoresence, and radiationless decay (via a conical intersection). The last of these relaxation methods provides an extreme test of the robustness of computational methods, since at the conical intersection there are (at least) two degenerate electronic states, and the Born−Oppenheimer approximation breaks down (see, e.g., refs 10 and 11).
Many studies have employed the CASSCF method to explore the excited state pathways (including conical intersections); there are far too many to include here. There are numerous examples also where TDDFT-based methods have been applied to such problems. In general, TDDFT results have been compared to CASSCF (where possible), with mixed accuracy. Minezawa and Gordon 12 compared spin-flip (SF) TDDFT minimum energy crossing point (MECP) conical intersection geometries of ethene with those determined at the MRCI and MS-CASPT2 levels of theory, with SF-TDDFT correctly predicting the three conical intersection geometries determined by the multireference methods. Filatov 13 studied the dependence of the choice of the density functional upon MECP geometry compared with various multireference approaches (CASSCF, CASPT2, MRCI), concluding that the BHHLYP hybrid functional performed the best, while popular contemporary functionals, such as M06-2X, perform relatively poorly. Nikiforov et al. 14 studied a group of small organic molecules, using the restricted ensemble-referenced Kohn−Sham (REKS) approach, comparing their results with MR-CISD calculated geometries. They determined average RMSD differences from the MR-CISD geometries of ∼0.1 Å, although the underlying MCSCF wave functions for some of the molecules had reduced active spaces due to technical limitations. Zhang and Herbert 15 compared SF-TDDFT calculated MECP conical intersection geometries of 9H-adenine to MR-CIS results, 16 noting that the difference between the two methods was nearly indistinguishable. Recently, Segarra-Marti et al. 17 studied the excited state decay of uracil and thymine cations, while including dynamical correlation using extended multistate CASPT2 (XMS-CASPT2). 18 They found that inclusion of dynamical electron correlation resulted in the separation of the energy levels of a "3-state" conical intersection, giving a different geometry and energy.
While it is desirable to use the highest-level theory possible to determine MECP geometries, for systems of interest in both chemistry and biology, it is not always possible to use multireference approaches. In the current study, we compare MECP conical intersection geometries calculated using SF-TDDFT, both with and without explicit calculation of nonadiabatic coupling terms, with both CASSCF and XMS-CASPT2. Our primary motivation is to determine the accuracy of SF-TDDFT S 1 /S 0 MECPs against XMS-CASPT2 for a set of medium-sized molecules, with particular emphasis on methods that do not require explicit computation of the nonadiabatic coupling terms. It is therefore useful to establish a protocol for the optimization of MECP geometries of larger molecules that treat electron correlation sufficiently. Boggio-Pasqua and Bearpark have recently investigated a similar approach to radical polycyclic aromatic hydrocarbons. 19 The rest of the paper is organized as follows: in section 2, we present an overview of the background theory of the approaches used; in section 3, we give the computational details; in section 4, we present the results of our comparisons, and in section 5 we give our conclusions.

BACKGROUND THEORY
We here present a brief overview of the background theory to the methods involved, in order to fully understand the key differences between each approach.
2.1. Branching Planes. In the vicinity of a conical intersection, two (or more) electronic states become degenerate, and the Born−Oppenheimer approximation breaks down. Consider the case where two electronic states, J and K, become degenerate; given N int internal degrees of freedom, the (N int −2)-dimensional space is known as the seam space, in which the two electronic states are degenerate, and the remaining two degrees of freedom are called the branching space. Within the branching space, the degeneracy of the Born−Oppenheimer surfaces is lifted by an infinitesimal shift in nuclear coordinates. The branching space is spanned by two vectors, g JK and h JK . The first is evaluated simply as the difference in gradient vectors of the two Born−Oppenheimer electronic states: The second vector, h JK , is known as the nonadiabatic coupling vector and is defined as The derivative coupling vector, d JK , can then be calculated as The topology of the seam space is determined by the relative orientation and magnitude of the two vectors g JK and h JK .

Multireference Approaches.
Assuming an appropriately chosen active space, the CASSCF and/or XMS-CASPT2 20−23 approach can be used to calculate both g JK and h JK analytically, using state-averaged wave functions (over the Born−Oppenheimer electronic states of interest then, for a CIS calculation, the off-diagonal term H JK (R) must be zero when an S 1 /S 0 conical intersection is calculated, due to Brillouin's theorem (this is not strictly true for TDDFT but is, in general, observed for functionals typically used with exact Hartree−Fock exchange). As a result, the coupling matrix elements of h JK vanish, and the degeneracy is only over one degree of freedom, resulting in an incorrect topology of the conical intersection. No such condition arises for two excited states becoming degenerate. One approach used to tackle this problem is SF-TDDFT, 24 in which the reference state (equivalent to the "ground-state") has a different spin multiplicity to the target states, hence the "ground-state" is also treated as an excited state, thus both g JK and h JK determine the topology around the conical intersection. 25 SF-TDDFT approaches can be used to calculate analytic derivative couplings. 25

Penalty Function Optimization.
Where analytic derivative couplings cannot be calculated, or where their calculation is expensive, more approximate methods can be used. The first considered here is the penalty-constrained optimization algorithm of Martínez et al. 26 In this approach, minimization of the objective function is performed, where α is a parameter employed to avoid singularities, and σ is a Langrange multiplier. The minimization of the penalty function is performed in an iterative manner for increasing values of σ. 2.3.3. Branching Plane Update. The second approximate approach considered in this study is the branching-plane update algorithm of Morokuma et al. 27 The mean energy gradient, G mean , is defined as 3. COMPUTATIONAL DETAILS XMS-CASPT2 calculations were performed for the molecules given in Table 1 (including active spaces). These molecules were chosen as representative of molecular structures found in molecular probes of biological function, including a model of long unsaturated lipid tails (2,4,6-octatriene). The basis sets were chosen to match those used in the references provided ( Table 1). The S 0 , S 1 , and S 1 /S 0 MECP geometries were optimized, using an average over the first two singlet excited states (in C 1 symmetry). The S 1 excited state geometry was used as an initial guess for the conical intersection geometry. In the absence of convergence, a small "kink" was added to the molecular structure to aid convergence by avoiding any limiting symmetry constraint (in the case of a ring structure, the ring was puckered by moving one atom 0.1 Å out of the plane of the ring). In all cases, the MECP geometries to be found were those that match the literature references given in Table 1. A real shift of 0.2 au was used in all XMS-CASPT2 computations. Density fitting was used in all XMS-CASPT2 calculations, employing the TZVPP-JKFIT density fitting basis set, except for fulvene, where the cc-pVDZ-JKFIT set was employed. MECP geometries were determined using the same gradient projection algorithm as for SF-TDDFT. 28 All XMS-CASPT2 geometry optimizations were performed using the BAGEL software. 29,30 CASSCF conical intersection geometry calculations were performed with the Molpro 2015.1 software 28,31 without the use of density fitting. In all cases, we were comparing SF-TDDFT approaches to XMS-CASPT2 for well-known MECPs, not trying to identify novel MECPs.
SF-TDDFT calculations were performed using the BHHLYP 40,41 and ωB97X 42 functionals and the basis sets given in Table 1. The BHHLYP functional has 50% Hartree− Fock exchange; such functionals are noted as successful within the spin-flip methodology. 24,43 The ωB97X functional was chosen as an example of a contemporary range-separated hybrid functional that performs well for a variety of applications. 44 MECP geometries were determined by analytical calculation of the nonadiabatic coupling matrix elements as discussed in section 2.3.1, 25 using the gradient projection algorithm of Bearpark et al. 28 We denote this as NAC for brevity. The penalty-constrained optimization approach (discussed in section 2.3.2; here defined as PC) 26 and branching-plane (discussed in section 2.3.3; here denoted as BP) 27 algorithms were used to determine MECP conical intersection geometries using SF-TDDFT without the explicit calculation of the nonadiabatic coupling terms. The spin-flip approach was used to determine the "reference" TDDFT state, as discussed in section 2. Each of the approaches considered here was performed within the Tamm-Dancoff Approximation, 45 due to restrictions in the implementation of the SF-TDDFT methodology. For the 5FC and azomethane molecules, convergence to the MECP structures using each of the SF-TDDFT approaches was poor. In these cases, increasing the basis set to 6-31G(d,p) and optimizing to the MECP geometry using NAC-BHHLYP gave a good starting geometry for each of the SF-TDDFT approaches with the 6-31G(d) basis set. The DFT and SF-TDDFT calculations were performed with the Q-Chem 5.0 software suite. 46

RESULTS AND DISCUSSION
In this section, we briefly compare each molecule individually with the calculated XMS-CASPT2 MECP geometries, before giving a more general discussion of the performance of the SF-TDDFT-based methods.
4.1. Fulvene. Fulvene has been widely studied as a benchmark for calculations of the MECP conical intersection between the S 0 and S 1 electronic states. 32 Given in Table 2 are selected geometrical parameters for the stable MECP. 32 Qualitatively, each structure calculated using CASSCF and SF-TDDFT methods matches the XMS-CASPT2 reference structure well. The calculated bond lengths of the CASSCF structure are within 0.01 Å of the XMS-CASPT2 structure, while the largest deviation for the SF-TDDFT methods is 0.04 Å. The relative energetics are given in Table S2 (Supporting  Information). Each of the methods using BHHLYP is very similar, with a deviation from the XMS-CASPT2 energies of ∼0.5 eV for the difference between the S 0 minimum and MECP energies.
4.2. 4ABN. Groups studying the photochemistry of 4ABN (and related amino-substituted benzonitriles) are primarily interested in the S 2 /S 1 conical intersection, as this plays a pivotal role in the presence or absence of dual fluorescence  (10,9) 6-31+G(d) 5FC 34 (8,7) 6-31G(d) 9H-adenine 35 (12,10) 6-31G(d,p) 2,4,6-octatriene 36 (6,6) 6-31+G(d) azomethane 37 (6,4) 6-31G(d) azoxymethane 37 (6,4) 6-31G(d) phenol 38 ( Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article bands in polar/nonpolar solvents. 33 Shown in Figure 2 is the XMS-CASPT2 S 1 /S 0 MECP geometry for 4ABN, with selected geometrical parameters given in Table 3. The XMS-CASPT2 and CASSCF geometries both exhibit a "boat-like" conformation, with the CN and −NH 2 groups both pointing away from the plane of the ring, while each of the SF-TDDFT methods has the −NH 2 group nearly planar with the ring. In particular, the −NH 2 group is completely planar for BP-ωB97X. The relative energetics are given in Table S3. Each of the SF-TDDFT approaches, except PC-BHHLYP, has a deviation from the XMS-CASPT2 relative energies of ∼0.8 eV (for the MECP). 4.3. 5FC. Blancafort et al. studied the photophysics of both cytosine and 5-fluorocytosine, characterizing both the S 1 /S 0 and S 2 /S 1 conical intersections to determine the nonradiative decay pathway. 34 Shown in Figure 3 is the XMS-CASPT2 S 1 / S 0 MECP geometry, along with geometrical parameters ( Table  4). The CASSCF geometry shows reasonable agreement with the XMS-CASPT2 geometry, although the C−H and N−H bonds remain planar with the ring, unlike the XMS-CASPT2 geometry. The SF-TDDFT methods have qualitative deficiencies in comparison to XMS-CASPT2; the nonplanar nature of the ring and functional groups is quite different than that observed for XMS-CASPT2 (see Figure 3). In particular, the C−O bond is too long for each of the SF-TDDFT approaches. Qualitatively, the NAC approach most closely resembles the XMS-CASPT2 geometry; both the penalty constrained and branching plane update algorithms deviate further from the XMS-CASPT2 geometry. The relative energetics are given in Table S4. Where the MECP geometries are qualitatively correct, the relative energetics are within 0.4 eV of the XMS-CASPT2 energies.
4.4. 9H-Adenine. Perun et al. discovered conical intersections for the radiationless decay mechanisms of the lowest energy 1 nπ* and 1 ππ* ( 1 L b ) electronically excited states of 9H-adenine, using CASSCF. 35 Single-point CASPT2 energies calculated at the S 0 equilibrium geometry predicted the 1 ππ* ( 1 L b ) as the lowest singlet excited state and the 1 nπ* as the S 3 state, separated by the 1 ππ* ( 1 L a ) state, in contrast to TDDFT and CIPSI results. 47 Given in Tables 5 and 6 are the calculated conical intersection geometrical parameters, for the conical intersections between the 1 ππ* ( 1 L b ) state and the ground state and the 1 nπ* state and the ground state,  respectively. The corresponding XMS-CASPT2 conical intersection geometries are shown in Figures S3(a) and S3(b), respectively. The SF-TDDFT approaches are qualitatively similar to the XMS-CASPT2 geometry, with the exception of the orientation of the C−H bond vector shown in Figure 4. There are some more significant quantitative differences (Tables 4 and 5); in particular, both CASSCF (for the 1 ππ* state) and PC-ωB97X (for the 1 nπ* state) exhibit large differences in the bond lengths compared to XMS-CASPT2.
4.5. 2,4,6-Octatriene. Chattopadhyay et al. characterized the S 1 /S 0 conical intersection on the photoisomerization pathway of 2,4,6-octatriene at the CASSCF/6-31G(d) level, along with S 0 and S 1 equilibrium geometries. 36 Shown in Figure 5 is the XMS-CASPT2 geometry, along with selected geometrical parameters in Table 7. Most of the methods qualitatively match the XMS-CASPT2 geometry, with reasonable quantitative accuracy. The exception is PC-ωB97X, which has a qualitatively different geometry ( Figure  5). Despite numerous efforts, including starting from the XMS-CASPT2 geometry and selecting different electronic states in the spin-flip procedure for the MECP optimization, the geometry shown in Figure 5 was consistently obtained for this method. This suggests that the penalty constrained algorithm with the ωB97X functional is a poor approach to find such MECPs; indeed, when the full nonadiabatic coupling terms are included, with either ωB97X or BHHLYP, then good qualitative and quantitative agreement with XMS-CASPT2 is observed.
4.6. Azomethane and Azoxymethane. Ghosh et al. studied the photoisomerization pathways of azomethane and azoxymethane, determining the S 1 /S 0 conical intersection geometries for each, respectively. 37 Shown in Figure S5 is the XMS-CASPT2 S 1 /S 0 conical intersection geometry for azomethane, along with selected geometrical parameters in Table 8. Each of the SF-TDDFT-based methods qualitatively matches the XMS-CASPT2 geometry, with good quantitative agreement, including the relative energetics (Table S7).
Given in Table 9 are selected geometrical parameters for azoxymethane (the XMS-CASPT2 S 1 /S 0 conical intersection is shown in Figure S6). As for azomethane, the SF-TDDFTbased methods show generally qualitatively correct conical intersection geometries, apart from the N1−O1 bond. The C− N−N−C backbone exhibits more of a kink in comparison to the XMS-CASPT2 geometry, which has a dihedral angle close to 172°, while the poorest performing SF-TDDFT methods show an angle of ∼155°(BP-BHHLYP and BP-ωB97X). Atom numbering taken from Figure 1.  Atom numbering taken from Figure 1.

Journal of Chemical Theory and Computation
pubs.acs.org/JCTC Article 4.7. Phenol. The photodissociation of the O−H bond of phenol has previously been studied, identifying the S 2 /S 0 conical intersection as a critical point on the pathway (as well as an S 2 /S 1 conical intersection). 38 Given in Table 10 are selected geometrical parameters, while the XMS-CASPT2 geometry is shown in Figure S7. Most methods show good quantitative agreement with the XMS-CASPT2 geometry, except for BP-BHHLYP, which, while appearing qualitatively correct, exhibits significant quantitative differences from the XMS-CASPT2 geometry.
4.8. SMAC. Zhao et al. investigated the photoinduced isomerization mechanism of SMAC. 39 In particular, they identified five MECP conical intersections using CASSCF. One of these is related to the excited-state intramolecular proton transfer (ESIPT) process, while the other four involve rotation (denoted as TW, for twist) around the CN bond (see Figure 1). We retain the authors' original naming convention for each of the conical intersection geometries here for convenience. For the ESIPT conical intersection, the XMS-CASPT2 geometry has the C−N−C(methyl) plane perpendicular to the plane of the aromatic ring ( Figure 6). The PC-BHHLYP, PC-ωB97X, and NAC-BHHLYP methods all qualitatively match the XMS-CASPT2 geometry, albeit with some significant differences quantitatively, especially the O−H distance (Table 11). The best quantitative correlation occurs     Atom numbering taken from Figure 1.

Journal of Chemical Theory and Computation
pubs.acs.org/JCTC Article for the PC-ωB97X functional. The BP-BHHLYP geometry shows an angle of ∼45°between the two planes formed by the ring and the substituent, while the BP-ωB97X and NAC-ωB97X geometries exhibit near-planarity ( Figure 6). The two XMS-CASPT2 conical intersection geometries, labeled "TWin1" and "TWin2", show the N−CH 3 perpendicular to the plane of the aromatic ring, with the methyl group pointing down or up, respectively (see Figures S8(b) and S8(c)). In all cases, the SF-TDDFT methods are qualitatively similar to the XMS-CASPT2 geometries, with fairly good quantitative agreement (Tables 12 and 13) for most geometrical parameters.
The XMS-CASPT2 geometries for the "TWout1" and "TWout2" conical intersections are similar to above, but the −OH proton points away from the nitrogen atom (while for "TWin1" and "TWin2" the proton points toward the nitrogen atom; see Figures S8(d) and S8(e)). Each of the SF-TDDFT approaches shows good qualitative agreement with XMS-CASPT2 for "TWout1", although some of the bond lengths are different by as much as 0.1 Å (Table 14). The ωB97X   Atom numbering taken from Figure 1. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article functional is qualitatively correct for "TWout2" compared to XMS-CASPT2 for each of the methods considered, but PC-BHHLYP and NAC-BHHLYP show very different features (Table 15 and Figure 7), with a significant kink in the ring at position 1 (see Figure 1 for atom numbering). 4.9. Discussion. The convergence of a conical intersection geometry optimization with SF-TDDFT was often challenging. However, for some of the SF-TDDFT approaches (and basis sets) considered here, the optimization problem seems almost pathological. The two molecules, 4ABN and azomethane, hint at such problems. In both cases, several of the SF-TDDFT methods failed to achieve convergence using the 6-31G(d) basis set (as used in refs 33 and 37), but the geometry converged smoothly with the larger 6-31G(d,p) basis set. Using this converged geometry as a guess, the 6-31G(d) calculations then swiftly converged, suggesting that the choice of method is relatively insensitive in the immediate region of a conical intersection; the problem is getting to such a region!     17 In some cases, there were significant differences between CASSCF and XMS-CASPT2. We expect the geometrical parameters coincident with the g and h vectors to be accurately described by CASSCF, but other bond lengths may be expected to be longer. Indeed, we found in a few cases differences between the calculated geometries from XMS-CASPT2 and CASSCF (e.g., 9H-adenine; see Table 5). In addition, the successful application of CASSCF (and, by extension, XMS-CASPT2) is limited by the choice of active space, which itself may be limited by (a) an inexperienced user and/or (b) technical limitations within a given piece of software. The second of these limitations is noted for 9H-adenine, where the authors were restricted to (6,6) for the conical intersection search, despite identifying the "correct" active space as 12 electrons in 10 orbitals. 35 In some cases in the current work, where the SF-TDDFTbased method failed to give a qualitatively correct geometry for a conical intersection, we tried starting from the CASSCF and/ or XMS-CASPT2 geometries and reoptimizing, e.g., 4ABN ( Figure 2) and 2,4,6-octatriene ( Figure 5). In each of these cases, the combination of algorithm and functional (plus basis set) led the geometry away from that determined by XMS-CASPT2 to the qualitatively incorrect structures observed. This was mainly observed when nonadiabatic coupling terms were neglected (i.e., the PC or BP approaches) or in the case of the ESIPT MECP of SMAC, with NAC-ωB97X. This agrees well with the work of Herbert et al., 48 who recommend the use of BHHLYP when identifying MECP structures with SF-TDDFT.
Given in Table 16 are the mean deviations, mean unsigned deviations, and maximum errors for the geometrical parameters. While the mean deviations look very encouraging, the mean unsigned error is a more realistic measure for each of the geometrical parameters. The maximum errors are due to the few cases noted above, in particular the failure of PC-ωB97X to correctly describe the 2,4,6-octatriene MECP. While spincontamination can be an issue with SF-TDDFT, in the case of 2,4,6-octatriene with ⟨S 2 ⟩ values stay well below the default thresholds (1.20), and the failure can thus be attributed to the penalty-constrained projection algorithm, rather than SF-TDDFT. While NAC-ωB97X fails to correctly describe the ESIPT MECP of SMAC, NAC-BHHLYP matches the XMS-CASPT2 geometry. Once again, spin-contamination is not an issue here; Herbert et al. recommend the use of BHHLYP when using SF-TDDFT, and our results confirm this. We also calculated s x and s y tilt parameters, 48 which confirmed each of   the MECP geometries obtained in this study had a peaked topology. The BP-BHHLYP and PC-BHHLYP approaches show generally good agreement with the XMS-CASPT2 computed MECP geometries; they also closely match the NAC-BHHLYP geometries, with less computational effort required. The BP-ωB97X approach also shows reasonable agreement, but the PC-ωB97X shows some significant deviations. The relative energetics for each molecule considered are given in Tables S2−S20 and Figure S9 in the Supporting Information. The qualitative picture here is generally good in comparison to the XMS-CASPT2 energies, with a few exceptions (related to the qualitative disagreement of the MECP geometries). For 9Hadenine, only CASSCF qualitatively matches the gap relative to the vertical excitation energy for the MECP energy; this is lower than the vertical excitation energy, while all of the SF-TDDFT approaches give a relative energy higher than the vertical excitation energy. The deviation from the XMS-CASPT2 energies in most cases is within 1 eV with largely the correct qualitative trend but with relatively low quantitative accuracy. We note that NAC-BHHLYP most closely follows the XMS-CASPT2 relative energy values ( Figure S9), while PC-BHHLYP and PC-ωB97X both outperform the branching plane update method.

CONCLUSIONS
We have studied different SF-TDDFT-based approaches for optimizing conical intersection geometries and compared them to the XMS-CASPT2 method. NAC-BHHLYP is the most reliable method for calculating the MECP geometries, but BP-BHHLYP and PC-BHHLYP also demonstrate good agreement, while having a substantially reduced computational cost in comparison to NAC-BHHLYP. Keal et al. 49 concluded that the penalty-constrained approach should only be used where full nonadiabatic coupling terms cannot be calculated; while we have demonstrated that PC-BHHLYP appears to be reliable in most situations, we would also agree that NAC-BHHLYP should be employed where possible. For reasonable relative energetics, only the NAC-BHHLYP approach can be recommended; thus, initial MECP optimization could be performed by PC-BHHLYP and refined by NAC-BHHLYP. The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.