Assessment of Approximate Computational Methods for Conical Intersections and Branching Plane Vectors in Organic Molecules the Predissociation Mechanism of the B̃ 2 a ′ State of Hco via the Conical Intersection with the X̃ 2 a ′ State Assessment of Approximate Computational Methods for Conical Inter

Articles you may be interested in Interpretation of the photoelectron, ultraviolet, and vacuum ultraviolet photoabsorption spectra of bromobenzene by ab initio configuration interaction and DFT computations An intraline of conical intersections for methylamine Theoretical study of the two-photon absorption properties of several asymmetrically substituted stilbenoid molecules J. Quantum-chemical computational methods are benchmarked for their ability to describe conical intersections in a series of organic molecules and models of biological chromophores. Reference results for the geometries, relative energies, and branching planes of conical intersections are obtained using ab initio multireference configuration interaction with single and double excitations (MRCISD). They are compared with the results from more approximate methods, namely, the state-interaction state-averaged restricted ensemble-referenced Kohn-Sham method, spin-flip time-dependent density functional theory, and a semiempirical MRCISD approach using an orthogonalization-corrected model. It is demonstrated that these approximate methods reproduce the ab initio reference data very well, with root-mean-square deviations in the optimized geometries of the order of 0.1 Å or less and with reasonable agreement in the computed relative energies. A detailed analysis of the branching plane vectors shows that all currently applied methods yield similar nuclear displacements for escaping the strong non-adiabatic coupling region near the conical intersections. Our comparisons support the use of the tested quantum-chemical methods for modeling the photochemistry of large organic and biological systems.


I. INTRODUCTION
8][9] The most important descriptors in this regard are the so-called branching plane vectors, 8 i.e., the atomic displacements along which the degeneracy of the electronic states at a CI is lifted.These vectors span the branching plane (BP) which comprises the possible directions for exiting the strong non-adiabatic coupling regime at a CI. 6,7 nceptually, multireference quantum-chemical methods of ab initio wavefunction theory are best suited to investigate CIs and non-adiabatic dynamics in molecules. 10Methods such as complete-active-space second-order perturbation theory (CASPT2) 11 and multireference configuration interaction with single and double excitations (MRCISD or MRCI, for brevity) 12 are capable of providing accurate computational results for small and mid-size molecules.With the growing interest in photobiology and in the photochemistry of large organic and inorganic compounds, for example, light-driven a) Electronic mail: mike.filatov@gmail.comartificial molecular machines, 13 there is a growing demand for simple yet sufficiently accurate computational methods for describing CIs and their BPs.
In recent years, there has been considerable progress in the development of such computational methods.In the domain of density functional theory (DFT), the use of the spinflip (SF) ansatz in the context of time-dependent DFT (SF-TDDFT) [14][15][16][17] allows the treatment of intersections between ground-and excited-state PESs. 18,19 n alternative approach based on ensemble DFT is the state-interaction state-averaged restricted ensemble-referenced Kohn-Sham (SI-SA-REKS) method, which has also been shown to be suited for describing CIs in organic molecules. 20][20][21][22][23] Semiempirical quantum-chemical methods 24,25 have become increasingly popular for investigating the non-adiabatic dynamics of excited states, in particular, the OM2/MRCI approach with an orthogonalization-corrected model Hamiltonian (OM2). 25,26 ][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44] Apart from a case study on adenine, 45 we are not aware of any systematic investigation on the ability of approximate quantum-chemical methods to properly describe the CIs and BPs of larger organic molecules (despite the growing number of such applications).This is the topic of the present article: we compare the results obtained from SF-TDDFT, SI-SA-REKS, and OM2/MRCI calculations with those from high-level ab initio MRCISD calculations covering a representative set of eight organic compounds with 12 CIs.In Sec.II, we briefly describe some key properties of CIs as well as the computational methods employed in this study.In Sec.III, we specify the corresponding computational details.In Sec.IV, we discuss the results obtained for the CIs, focusing on their geometries, relative energies, and branching planes.Our motivation is to provide guidance in selecting suitable computational approaches for non-adiabatic simulations of excited-state dynamics in photochemistry and photobiology.

II. COMPUTATIONAL METHODOLOGY
Here we review the salient features of conical intersections as well as the computational methods used in this work.For more detail, we refer the readers to the original publications cited below.

A. Conical intersections and branching plane vectors
][48] Two adiabatic electronic states representing the solutions of the secular problem (1) become degenerate if the following conditions are fulfilled: 3,8,49 H mm − H nn = E m − E n = 0, (2a) where Ĥ is the electronic Hamiltonian in Born-Oppenheimer approximation, and its matrix elements H mn are formulated in terms of (in general, arbitrary) diabatic orthogonal states m and n .As these conditions can be satisfied in the space of N−2 internal molecular coordinates, [46][47][48] CIs occur in polyatomic molecules with three or more atoms.In the space of all nuclear coordinates, the degeneracy between the adiabatic electronic states m and n is lifted along two directions defined by the gradients of the conditions in Eqs.(2a) and (2b), 3,8,[46][47][48][49][50][51] x 1 = ∇ Q (H mm − H nn ), (3a) In Eq. ( 3), ∇ Q denotes to differentiation with respect to the nuclear coordinates.In the vicinity of a crossing point, the adiabatic electronic states depend linearly on the displacement along the x 1 and x 2 vectors, and therefore the PESs near a CI have the topography of a double cone. 8,46 [51] For the mechanism of photochemical reactions, the branching plane spanned by x 1 and x 2 defines all possible directions of exiting the strong non-adiabatic coupling region and then propagating on one of the intersecting PES, 4,5,50 thus playing a role similar to the transition vector in the transition state theory of ground-state chemical reactions.
In Eq. ( 3), the branching plane vectors are defined via the diabatic states m and n .Alternatively, the BP vectors can be also defined via the adiabatic states, see Eq. ( 4), As first demonstrated by Atchity et al., 8 the orientation of the BP vectors (4) may not coincide with the x 1 and x 2 vectors of Eq. ( 3), even when obtained with the same computational method.The BP vectors (4) can be changed by a rigid rotation through an arbitrary angle within the branching plane without changing the plane itself.Therefore, although the BP remains invariant under such a transformation, it is difficult to directly compare the individual BP vectors obtained from different computational methods.
To simplify a visual comparison between the BP vectors produced by the different methods employed presently, the normalized BP vectors can be subjected to a similarity transformation SRS −1 , where R is a two-dimensional (2D) orthogonal rotation matrix and S is a 2D shear transformation; the latter is needed because the BP vectors are not necessarily orthogonal.The angle of the 2D rotation is determined from the condition to maximize the projection of the vector a (see Figure 1) onto the vector a, when the two sets of BP vectors {a,b} and {a ,b } are obtained using two different computational methods labeled I and J.The shear transformation is set up by the constraint of keeping the inner product (a • b) between the BP vectors the same before and after the transformation.By applying the above transformation, the BP vector pairs obtained from different methods are reoriented such that the x 1 vectors are aligned.
To facilitate a quantitative comparison between the BP vectors from different methods, we use visual comparisons and the following measures of consistency between the branching planes: (1) projections of the BP vectors obtained using method I onto the whole branching plane obtained by method J, and (2) projection of the rectangle enveloped by the BP vectors from method I onto that from method J, see Figure 1.
For an arbitrary unit vector a , the norm of a projection p(a ) onto a plane spanned by (generally non-orthogonal) unit vectors a and b, see Figure 1, is given by For orthogonal vectors a and b, Eq. ( 6) simplifies to the usual scalar product.Note, however, that the BP vectors defined in Eqs. ( 3) and ( 4) are not required to be orthogonal.For a plane BP I (defined by two unit vectors, a and b ), the projection onto the plane BP J (defined by the vectors a and b) can be defined via the area of a rectangle enveloped by the projections p(a ) and p(b ) of the individual vectors.The area of a rectangle enveloped by a and b on the plane BP I is given by and the area of its projection onto the plane BP J is given by As the unit vectors a and b are not, in general, mutually orthogonal, the area s I is not equal to unity.Therefore, instead of p J (s I ), we use the ratio r IJ = p J (s I )/s I as a measure of parallelity between the planes BP I and BP J .When the two vector pairs (a,b) and (a ,b ) are equivalent (in the sense that the two planes are parallel and the vector pairs may only be rigidly rotated within the plane with respect to one another), the ratio r IJ is equal to unity, which is its maximum attainable value.The same is also true for r JI defined in analogy to r IJ .Obviously, the minimum possible value of r IJ (and r JI ) is zero, for two mutually orthogonal planes.As the BP-onto-BP projection can be characterized by a single number, it is the measure that we use in the remainder of the paper to compare the branching planes produced by different methods.The projections of the individual BP vectors are documented in the supplementary material. 52ptimizations of minimum-energy conical intersections (MECIs) may result in different molecular orientations when using different quantum-chemical methods.Therefore, the optimized MECI geometries were aligned using the Kabsch algorithm as implemented in the VMD suite of programs. 53,54  discussed below in more detail, the MECI geometries obtained from different quantum-chemical methods are quite similar to each other.Hence, it can be expected that the effect of geometrical differences on the BP projections r IJ is minimal and that differences between the BPs produced by different quantum-chemical methods will thus mainly originate from the different shape of the intersecting PESs near the MECI points.

B. REKS and SA-REKS
A complete description of the REKS method can be found in the original literature, see Refs.55-57, and thus only its most relevant features will be repeated here.1][62][63] The ensemble representation for the density leads to fractional occupation numbers (FONs) of a few frontier KS orbitals. 60,61,63,64 Th REKS ground-state energy is minimized with respect to the KS orbitals and the FONs of the frontier active orbitals.The currently employed version of the REKS method treats two fractionally occupied frontier orbitals, as, for example, in a diradical state resulting from (near) degeneracy of the bonding π and the anti-bonding π * frontier orbitals of an alkene near a ca.90 • twist around the double bond. 65By analogy with the CASSCF method, the current version of REKS is dubbed REKS(2,2) to indicate that fractional occupations are considered for two active orbitals which contain a total of two electrons.
Drawing an analogy with wavefunction theory, REKS(2,2) describes a strongly correlated system for which a model two-configuration wavefunction can be introduced, 65 where ϕ a and ϕ b are the frontier (active) orbitals with FONs n a and n b , and the unbarred and barred orbitals are occupied with spin-up and spin-down electrons, respectively.For such a system, the lowest singlet excited state can be approximated by an open-shell singlet wavefunction (10), 65 7][68] However, direct application of the variational principle to an individual excited state (which possesses the same spin and space symmetry as the ground state) is not formally permitted in DFT 69 and may lead to artifacts in practical calculations. 70ith the use of the ensemble formalism, 71 this limitation is bypassed by applying the variational principle to an ensemble of the ground state (described by REKS(2,2)) and the excited state (described by ROKS) which leads to the state-averaged (SA) REKS method. 72In SA-REKS, the energy of an ensemble of states 11) is minimized with respect to the KS orbitals, which are common for both states, and with respect to the FONs of the active orbitals in the REKS(2,2) treatment.In Eq. (11), equal weighting factors C 0 and C 1 are employed.Having completed the KS orbital optimization, the individual energies, E REKS(2,2) 0 and E ROKS 1 , are calculated using these KS orbitals.
The SA-REKS method describes the ground state (S 0 ) and the lowest singlet excited (S 1 ) state of a homosymmetric biradical. 65In a heterosymmetric biradical, generated either by chemical substitution or by an asymmetric geometric distortion of a homosymmetric biradical, the states given in Eqs. ( 9) and ( 10) can mix with each other.This mixing becomes important in the vicinity of conical intersections. 73It can be taken into account in the SI-SA-REKS method 20,74 via a simple 2 × 2 secular equation, as in Eq. (1), in which the diagonal matrix elements are given by the energies E REKS(2,2) 0 and E ROKS 1 , while the off-diagonal matrix element is calculated using the Lagrange multiplier W ab between the SA-REKS active orbitals, Equation ( 12) has been derived 20,74 using the Slater-Condon rules and the variational condition for the SA-REKS orbitals. 55,67,68 Povided that equal weighting factors are employed in Eq. ( 11), the SA-REKS orbitals remain unaffected by the application of the state interaction procedure. 20,74 t is also worth noting that the SI-SA-REKS energy expression can be derived using the adiabatic connection argument.The interested reader may refer to the recent work of Fromager and Franck, 75 in which an energy expression similar to the one used in SI-SA-REKS was obtained by constructing a generalized adiabatic connection path between non-interacting and fully interacting ensembles of the ground and excited states of the type ( 9) and (10); cf.Sec.III in Ref. 75.
In the framework of the SI-SA-REKS method, a CI occurs when the energies E REKS(2,2) 0 and E ROKS 1 become equal to one another and the coupling matrix element (12) vanishes.Thus, the BP vectors, x 1 and x 2 , can be formally defined by Eq. (3), where H mm and H nn are replaced by E REKS(2,2) 0 and E ROKS 1 , respectively, and H mn is replaced by H 01 from Eq. (12).Currently, analytic derivatives of the SI-SA-REKS energies with respect to the nuclear coordinates are not yet implemented and hence differentiation has to be carried out numerically.The SI-SA-REKS method has been applied to the calculation of conical intersections in a number of organic molecules and models of biological chromophores and has been found to be accurate and reliable in reproducing the energetic and geometric characteristics of conical intersections. 20

C. SF-TDDFT
Spin-flip time-dependent density-functional theory couples electronic configurations arising from excitation operators with M s = ±1, unlike the usual formulation of TDDFT, in which only spin-preserving excitation operators with M s = 0 are allowed.5][16][17] The SF-TDDFT scheme of Ziegler and Wang relies on the non-collinear spin DFT framework, which operates with KS spinors, rather than spin-free orbitals, and where both α and β spin components are present such that the exchange-correlation density functional E xc [ρ αα , ρ ββ ] depends on the diagonal part of a 2 × 2 density matrix By performing a change of variables of the type one can define the exchange correlation functional as where the total density is ρ( r) = ρ αα ( r) + ρ ββ ( r) and the square of magnetization is The collinear limit can then be recovered by setting the magnetization to the spin-density limit, s → ρ αα − ρ ββ .In this limit, ρ +/ − → ρ α/β , thus recovering the usual functional form.
In the non-collinear framework, the exchange-correlation kernel of TDDFT is generalized to include both spinpreserving and spin-flip excitations, By taking the collinear limit, the surviving terms are In this block-diagonal kernel, one can easily recognize the 2 × 2 spin-preserving excitation block, corresponding to the usual TDDFT excitations, and the spin-flip block.
The reference in a two-electron two-orbital model, this gives rise to two closed-shell configurations (double excitation and ground state) and two open-shell configurations which can be coupled to singlet and triplet configurations (see Fig. 2).Outside this model, open-shell spin-uncompensated configurations are present so that spin-contaminated excitation energies are obtained, which can be partially purified by an a posteriori correction scheme. 76he SF-TDDFT scheme has several advantages over regular TDDFT, especially when describing a conical intersection region involving the ground state.First, the reference state is a triplet, which satisfies more easily the non-interacting pure state v-representability requirement of KS DFT.Second, SF-TDDFT includes some extra doubleexcitation character in the excited states due to the configuration arising from the a † β i α excitation.The coupling of this configuration with the ground-state configuration (i † β a α excitation) can give the correct dimensionality of the conical intersection region.This is a clear advantage over the usual spin-preserving TDDFT, in which the ground state and the excited states are decoupled. 77Recently, Truhlar and co-workers have introduced an ad hoc coupling term between the ground and the singly excited configurations in the context of spinpreserving TDDFT, thus recovering the correct dimensionality of the CI seam. 78However, the presence of this coupling term violates Brillouin's theorem for the KS orbitals used in the configuration interaction with single excitations, which is employed on top of the usual TDDFT formalism. 78The coupling introduced in this way is in general too small, because it does not recover the effects of doubly excited configurations that are responsible for the interaction between the ground and excited configurations in the context of the exact theory. 79SF-TDDFT is able to correctly introduce the ground and excited state coupling, albeit only within the restricted configuration space represented in Fig. 2. In general, most CIs are well represented by SF-TDDFT, but more complicated intersections are described only approximately. 18One of the most important problems is that SF-TDDFT states are frequently spinmixed.In order to correct this problem, an a posteriori spin purification scheme has been proposed. 76However, when us-ing this scheme, the correct dimensionality of the CI seam is lost. 23ithin the SF-TDDFT approach, conical intersections can be located using the branching space update method of Maeda, Ohno, and Morokuma, 80 which employs a projected gradient algorithm for optimizing the branching plane.The x 1 -vector is calculated by Eq. (4a) and the x 2 -vector is approximated by an iterative update scheme.Note that in this approach the x 1 and x 2 vectors are defined for the adiabatic states and are kept orthogonal, although the exact BP vectors need not necessarily have this property. 8

D. MRCISD
A detailed description of the MRCI approach can be found, e.g., in Ref. 81; here we present a brief overview.In MRCI, the reference wavefunction is represented as a combination of several configuration state functions (CSFs).CSFs are linear combinations of Slater determinants with the same spatial-orbital occupations; they are eigenfunctions of the S 2 and S z operators, where the total spin of the MRCI wavefunction defines the total spin of the CSFs.The MRCI wavefunction is expanded in terms of Nth-order excitations of each CSF, where N electrons are promoted from occupied to vacant orbitals, within the spin and symmetry constraints of the system.In the presently employed MRCI approach, we allow all single and double excitations from the multiconfigurational reference wavefunction (MRCISD).
The CI expansion space is generated using the Graphical Unitary Group Approach (GUGA) of Shavitt. 82,83 n GUGA, like in any graphical representation of the CI wavefunction, 82 the orbitals are distributed into groups and the occupation of each group is restricted within the set of expansion determinants.This occupation-based approach to construct the CI wavefunction is quite convenient in multireference calculations, as the latter typically divide the orbital space into several subsets of orbitals with restricted occupations.Moreover, it eliminates redundancy problems that may occur in MRCI when constructing the CI expansion space by exciting the reference CSFs.
The MRCISD treatment uses an adiabatic basis.The energy gradients, and thus the vector x 1 , are calculated analytically as described in Ref. 84.The vector x 2 is split into two components 85,86 Here accounts for the change of the CI vectors C m upon geometry changes, and arises from the geometry dependence of the underlying CSF expansions.The subscript CSF in Eq. ( 21) denotes a scalar product in CSF space, and subscript r in Eq. ( 22) implies integration over electronic coordinates.

E. OM2/MRCI
We give only a brief outline of the OM2/MRCI method here, more details are found in Refs.24, 25, 87, and 88.The orthogonalization-corrected OMx methods are semiempirical methods based on the neglect of diatomic differential overlap (NDDO).In the variational determination of the molecular orbital (MO) coefficients C and MO energies E from the Fock matrix F, the NDDO integral approximation directly leads to the standard eigenvalue problem, FC = CE, in contrast to the ab initio MO secular equations, FC = SCE (involving the overlap matrix S), which need to be transformed to the eigenvalue problem by an explicit orthogonalization step.The Pauli exchange repulsion and the asymmetric splitting of bonding and antibonding orbitals arise from these orthogonalization corrections in the ab initio framework.Standard MNDO-type semiempirical methods formally neglect these and related effects, which is believed to be responsible for many deficiencies of these standard methods. 24he OMx methods include orthogonalization corrections explicitly in the electronic calculation.The OM2 variant is the most complete model of the OMx family and offers a good compromise between accuracy and computational efficiency, especially for excited states. 89It includes valence-shell orthogonalization corrections both in the one-center part 87 and in the two-center part of the core Hamiltonian. 88,90 e MOs obtained from the OM2 calculation are used as input in a subsequent GUGA-MRCI calculation. 91Conceptually, the semiempirical OM2 MO calculation already accounts for a large fraction of dynamical correlation, so that the MRCI treatment mainly needs to recover non-dynamical correlation effects.This can normally be achieved by using rather small active spaces, a small number of reference configurations, and only single and double excitations in the CI expansion (MR-CISD).
The OM2/MRCI formalism uses an adiabatic basis, and hence the BP vectors are defined by Eq. (4).The vector x 1 between two adiabatic states I and J is calculated using the analytic expression where q α represents an individual nuclear coordinate, C m are the CI coefficients of state m, and H is the CI electronic Hamiltonian matrix.The x 2 vector is calculated using a previously implemented semi-analytic gradient module, 92 which utilizes transition density matrices according to the formulation of Lengsfield III and Yarkony. 93

III. COMPUTATIONAL DETAILS
The present DFT calculations employed the 6-31+G** basis set 94 and the BH&HLYP 95 hybrid density functional.In a recent study of conical intersections in organic molecules, 20 it was found that the BH&HLYP functional reproduces the geometries and relative energies of conical intersections obtained from high-level ab initio MRCI calculations more accurately than some other density functionals (including longrange corrected and meta-GGA functionals).
The REKS(2,2), SA-REKS, and SI-SA-REKS calculations were performed using the COLOGNE2012 program. 96he ground-state equilibrium geometries were optimized with the REKS(2,2) method.The geometries of conical intersections were obtained with the use of the CIOpt program, 97 which employs a penalty function minimization in connection with numerically calculated gradients for the ground and excited states.In the REKS calculations, the S 0 /S 1 energy gap at the optimized CI geometry was less than 0.01 kcal/mol in each case.The numerical integrations employed grids comprising 75 radial points and 302 angular points per atom, and a SCF convergence criterion of 10 −8 for the density matrix was applied in all REKS calculations.When computing the BP vectors by numerical differentiation as described in Sec.II B, the SCF convergence threshold was tightened to 10 −10 and the increment in the atomic coordinates was chosen to be 0.015 Å.
The SF-TDDFT calculations were done with the GAMESS-US program 98,99 using a grid of 96 radial points and 302 angular points for the ground-state calculations and 48 radial points and 110 angular points for the TDDFT calculations.The density convergence threshold for the density matrix was set to 10 −6 .The ground-state equilibrium geometries were optimized using the SF-TDDFT analytic gradients.The conical intersections were located with the branching space update method of Maeda, Ohno, and Morokuma. 80In this approach, the x 1 vector is calculated analytically, while the x 2 vector is approximated through updates.The S 0 /S 1 energy gap at the optimized CI geometries was always less than 0.08 kcal/mol.
1][102][103][104] The initial stateaveraged MCSCF calculations (with the ground and the relevant excited states included into the averaging) employed molecule-specific active spaces, namely: 2 electrons in 2 orbitals (2,2) for methyliminium; (4,4) for ethylene, butadiene, PSB3, and styrene; and (6,5) for ketene.The large size of stilbene and HBI forced us to restrict the MCSCF active space to (2,2).Where available, symmetry constraints were employed when calculating the FC conformations; for most of the molecules, C s symmetry was used, except for ethylene (D 2h ), trans-butadiene (C 2h ), cis-butadiene (C 2v ), and cis-stilbene (C 1 ).The MRCISD calculations employed the same active spaces as the preceding MCSCF treatments.All the results were obtained with the 6-31+G** basis, except for stilbene and HBI where we could not afford diffuse functions and thus used the 6-31G** basis.In the MRCISD optimization of conical intersections, the convergence criterion was an energy gap of less than 10 −4 E h (ca.0.06 kcal/mol) between the two intersecting states.Both the x 1 and x 2 vectors were calculated analytically, and the geometries were updated using a modified version of the direct inversion in the iterative subspace (GDIIS) algorithm. 85he OM2/MRCISD calculations were performed with the MNDO program. 105The GUGA-MRCI calculations employed three reference configurations (closed-shell, single HOMO-LUMO excitation, and double HOMO-LUMO excitation) and molecule-specific active spaces: (2,2) for ethylene and methyliminium; (4,4) for butadiene and PSB3; (6,5) for ketene; (6,6) for styrene; (8,8) for stilbene; and (10,9) for the anionic HBI chromophore.Where available, symmetry constraints were employed when calculating the FC conformations; for most of the molecules, C s symmetry was used, except for trans-butadiene (C 2h ), cis-butadiene (C 2v ), and cisstilbene (C 1 ).Conical intersections were optimized using a modified version 106 of the Lagrange-Newton algorithm proposed by Manaa and Yarkony, 107 with a diagonal initial Hessian containing empirical values for the corresponding internal coordinates. 106The x 1 and x 2 vectors were calculated analytically at each optimization step, while the Hessian was updated according to the BFGS procedure. 106The trust radius for each geometry update was set to 0.1 Å.The final energy gap at the optimized CI geometry was always less than 10 −5 kcal/mol.
The BP vectors obtained by the methods outlined above were normalized against their respective Frobenius norms.The following units are adopted throughout the remainder of this article: geometric parameters in Å, total energies in Hartree atomic units (a.u.), and relative energies in electron volts (eV).

IV. RESULTS AND DISCUSSION
Most of the conical intersections considered in this work arise from crossings that occur in π -conjugated systems along the internal coordinate that describes double bond torsion.[110] These CIs can be classified as twisted-pyramidalized (tw-pyr, for brevity) or as twisted-bond_length_alternating (tw-BLA) depending on the relative preference for the two bond-breaking mechanisms; the tw-pyr CIs are typical of molecules with dominant homolytic bond breaking, while the tw-BLA CIs occur in molecules for which both π -bond breaking processes are nearly isoenergetic. 22,108 other type of CIs commonly occurring in organic molecules can be vaguely classified as n/π CIs that originate from the crossing between an electron configuration (n 2 π * 0 ) with a doubly occupied lone pair and a singly excited (n 1 π * 1 ) configuration.In the presently studied molecules, such CIs occur in ketene, 117 ethylene (the ethylidene or methylcarbene CI), 118 and methyliminium (the methylimine CI), 108,109 see Figure 3.
The two types of CIs probe different aspects of quantumchemical computations.The correct description of the CIs occurring along the π -bond torsion coordinate requires a balanced description of the relative stability of covalent (diradical) and ionic configurations, while for the n/π 's CIs, the relative stability of configurations involving lone-pair excitation must be treated properly.
The set of molecules shown in Figure 3 is by no means exhaustive; for example, cyclic molecules like nucleobases are not present.However, for such molecules, the ab initio MRCISD reference calculations with sufficiently large active spaces and basis sets would be computationally too demanding.Even for the molecules in Figure 3, the ab initio MRCISD calculations were carried out with the rather small 6-31+G** basis set to make optimization of the conical intersections feasible for the larger molecules.In case of the smaller molecules  (ethylene and methyliminium), we would be able to afford ab initio MRCISD optimizations with larger basis sets, but for the sake of consistency, we prefer to employ the 6-31+G** basis throughout to allow for evenhanded and unbiased comparisons between the methods.Given these limitations, we emphasize that the current ab initio MRCISD/6-31+G** results cannot be regarded as having the status of highly accurate benchmark data; they are used in the following as reference values obtained at a uniform and reasonably reliable ab initio level that reflects the current state of the art.The results of the calculations are summarized in Figure 3 and in Tables I and II; the latter present the relevant relative energies and comparisons between the computed CI geometries and branching planes, respectively.
In an overall assessment, all the methods provide consistent descriptions of the MECI geometries and BP vectors for the test molecules considered presently.On average, the RMSDs of the MECI geometries are on the order of 0.1 Å and no conspicuous discrepancies between these geometries are seen in Figure 3.The same is true for the BP projections which for all pairs of methods vary between 0.7 and 0.9, thus indicating a consistent description of the nuclear motions that lead to an escape from the CI region.This would seem to suggest that the methods tested in this work should show similar behavior in non-adiabatic molecular dynamics (NAMD) simulations.One may also expect that this will hold true for larger molecules, for which the application of first-principles methods becomes prohibitively costly.In the following, we discuss the results for the individual molecules (using the abbreviations SSR and SF for SI-SA-REKS and SF-TDDFT, respectively).

A. Ethylene
The π → π * (HOMO→LUMO) excitation of ethylene results in rapid photoisomerization, [118][119][120] and two distinct MECIs can be identified along the CI seam. 7,118 he twpyr MECI (see Figure 3) arises from a crossing between the PESs of the electronic states for homolytic and heterolytic πbond breaking. 20,22,73,108,109 Theheterolytic π -bond breaking process is much less favorable than the homolytic one. 20,22 cording to a previously proposed rule, 20,22    b BP projections are given as a matrix with the off-diagonal elements r IJ representing the projection of the BP obtained using method I onto the BP obtained using method J.
MECI occurs in the proximity of a ground-state conformation that corresponds to heterolytic π -bond breaking and features strong pyramidalization to stabilize the shift of the π electron pair toward one of the carbon atoms.
The best estimate of the π → π * vertical excitation energy of ethylene is 7.8 eV. 111 The energy of the tw-pyr MECI with respect to the ground-state minimum is much lower: the values currently predicted by the different methods lie in a narrow energy range between 4.8 and 5.0 eV.The BP vectors of this MECI are shown in Figure 4, where the MRCISD, SF, and OM2/MRCI vectors are aligned with the SSR vectors as described in Sec.II.The x 1 vector mainly corresponds to the pyramidalization of one of the carbon atoms and points into the direction that connects the transition states of the heterolytic and homolytic π -bond breaking on the S 0 surface. 4,22 e x 2 vector reflects the torsion of the π -bond and is thus associated with the photoisomerization reaction. 4,22 he character of the two BP vectors is consistently reproduced by all the methods, and the overlaps between the BPs (i.e., the r IJ projections) are rather large, on the order of 0.8-0.9,see Table II.The tw-pyr MECI geometries produced by the different methods agree well with each other, as evidenced by the RMSDs in Table II; compared with the reference ab initio MRCISD geometry, the RMSDs for SSR, SF, and OM2 are 0.0312, 0.0766, and 0.0750 Å, respectively.
The n/π -type MECI along the CI seam of ethylene occurs near the ethylidene (or methylcarbene) geometry that results from a H-atom transfer from one carbon to the other one. 118Under the constraint of C s symmetry, this MECI arises from a crossing between the 1 A and 1 A states of ethylidene, which would be linear in C s symmetry but becomes conical in the unsymmetrical C 1 case. 117The x 1 and x 2 vectors of this MECI correspond to an a mode (bending of the H-C-C valence angle) and an a mode (wagging the terminal H atom), respectively. 117All the currently applied methods consistently reproduce the ethylidene MECI geometry: the RMSDs are typically below 0.1 Å, and the BP vectors are qualitatively consistent with each other.The numerical values of the BP projections r IJ are however somewhat lower than in the case of the tw-pyr MECI.As shown in Table II of the supplementary material, 52 it is the projections of the x 2 vector that are responsible for the decreased overlap between the BPs produced by different methods.Most likely, this is caused by rather flat S 0 and S 1 PESs in the direction of this vector.For instance, in the case of the SSR method, the norm of the x 2 vector (0.0010 Hartree/Bohr) is ca.40 times smaller than that for the tw-pyr MECI (0.0443 Hartree/Bohr); this implies that the derivatives, Eq. (3), in this direction are small, and any tiny differences are amplified when normalizing the BP vectors.Nevertheless, in spite of the existing quantitative discrepancies, the BP vectors FIG. 4. BP vectors of twisted-pyramidalized (upper panel) and ethylidene (lower panel) MECIs of ethylene calculated using the methods employed in this work and aligned using the similarity transformation described in Sec.II.J. Chem.Phys.141, 124122 (2014) of the ethylidene MECI produced by the different methods are qualitatively similar.It is also noteworthy that the ethylidene MECI does not play an important mechanistic role in the photoisomerization, as it is only very rarely visited by nuclear trajectories initiated at the FC point. 7Therefore, the discrepancies observed for this MECI should not lead to any significant consequences for the dynamic description of ethylene photoisomerization.
In our calculations, we employ a relatively modest 6-31+G** basis set, and it is thus of interest to check the effects of basis set extension on the computed MECI geometries and relative energies.In Ref. 20, the sensitivity of the SSR method to basis set extension was studied, and only minor changes were found in the geometries and energies of the tw-pyr MECI of ethylene and stilbene.Here, we report the effect of basis extension and active-space variation on the geometry and relative energy of the tw-pyr MECI of ethylene obtained by the ab initio MRCISD method, see Table III.Evidently, enlarging the active space from (2,2) to (4,4) only leads to slight changes of the MECI energies (tw-pyr MECI, decrease by 0.07 eV; eth MECI, increase by 0.15 eV), while basis set extension from 6-31+G** (double-ξ quality) to 6-311+G** (triple-ξ quality) 94 causes a slight decrease by ca.0.1 eV.Overall the variations in the MECI energies remain below 0.2 eV.The optimized MECI geometries show almost no dependence on the chosen basis set and active space for both MECIs; the largest RMSD observed in any of these calculations is as small as 0.0042 Å.These findings suggest that the use of a double-ξ quality basis set in ab initio MR-CISD calculations provides a description of MECIs, which is not substantially altered by moderate basis set extensions that keep the computational effort within reasonable bounds.The use of much larger basis sets for ab initio MRCISD MECI optimizations is currently not feasible (especially not for the larger molecules of our test set).

B. Methyliminium
Methyliminium CH 2 NH + 2 differs from ethylene by the presence of a strong electron-withdrawing substituent, the cationic nitrogen atom.The heterolytic C=N π -bond breaking becomes nearly isoenergetic with the homolytic bond breaking and the tw-pyr MECI changes its character to tw-BLA. 20,22 he tw-BLA MECI occurs at an energy ca. 1 eV lower than the tw-pyr MECI of ethylene (see Table I).Due to the stabilization of the CT configuration (C + -N 0 ), the tw-BLA MECI lies even lower than the methylimine MECI which has very little ionic character and thus does not benefit from the stabilization of the CT configuration.
All methods predict the BLA mode (coupled with some pyramidalization at the carbon end) for the x 1 vector and the twist mode for the x 2 vector, see Figure 5. OM2/MRCI, SSR, and SF yield some pyramidalization of the carbon atom, which is smallest for the OM2/MRCI method and largest for the SF method.The ab initio MRCISD tw-BLA MECI does not feature any pyramidalization.These differences between the methods arise from differences in the predicted relative stabilities of the ionic and covalent electronic configurations. 20,74 he DFT-based methods tend to somewhat overestimate the stability of the CT configurations, and hence they invoke a moderate pyramidalization at the carbon atom to destabilize the CT configuration. 20he methylimine MECI of CH 2 NH + 2 is qualitatively similar to the ethylidene MECI of ethylene and occurs at a geometry with one of the amine hydrogen atoms transferred to the carbon atom, see Figures 3 and 5. Similar to ethylene, the x 1 vector of this MECI corresponds to H-N-C bending (a mode in C s symmetry) and the x 2 vector describes a wagging motion of the NH + group (a mode in C s symmetry).The S 0 and S 1 energy profiles in the direction of the x 2 vector are quite flat (e.g., the norm of the SSR x 2 vector is only 0.0040 Hartree/Bohr as compared to the norm of the x 1 vector of 0.1118 Hartree/Bohr); this leads to a large relative variation of the x 2 vectors across the methods.The MECI geometries obtained from the different methods agree with one another within 0.08 Å, and the branching planes from different methods are qualitatively similar (see lower panel of Figure 5).

C. Styrene
The presence of a phenyl group at one end of the ethylenic moiety in styrene leads to the occurrence of two non-equivalent tw-pyr MECIs, which differ in the pyramidalization of the methylenic carbon atom (MECI 1 ) and of the  benzylic carbon atom (MECI 2 ).Qualitatively, the two MECIs are similar to the tw-pyr MECI of ethylene; however, the weak electron-donating character of the aryl group slightly stabilizes the ionic configurations originating from the heterolytic breaking of the ethylenic π -bond and the two MECIs thus occur at a slightly lower energy with respect to the S 0 minimum than in ethylene.MECI 1 experiences a somewhat stronger stabilization than MECI 2 , because, in the CT configuration, the positive charge can be delocalized into the aryl ring. 20,22 r both tw-pyr MECIs of styrene, the x 1 vector corresponds to pyramidalization of one of the ethylenic carbons, the methylenic carbon (MECI 1 ) or the benzylic carbon (MECI 2 ), and the x 2 vector describes torsion about the ethylenic C=C bond (see Fig. 6).This picture is very similar to ethylene, and all methods yield similar geometries and BPs for the two MECIs, see Table II.The OM2/MRCI, SSR, and SF calculations give similar relative energies for the two MECIs, while the ab initio MRCISD MECIs lie ca.0.4 eV higher in energy (Table I).It is conceivable that the use of a more extended basis set (larger than 6-31+G**) might bring the ab initio MRCISD results into better agreement with the results from the approximate methods.

D. Stilbene
The tw-pyr MECI of stilbene is qualitatively quite similar to the tw-pyr MECIs of ethylene and styrene, see Figures 3 and 7.The presence of the second phenyl ring does not result in a noticeable alteration of the energy level of the tw-pyr MECI, e.g., as compared to MECI 1 of styrene; all methods predict this MECI to lie at ca. 4.2-4.5 eV above the S 0 minimum (trans-conformation).We note that the ab initio MRCISD calculations of stilbene employed the smaller 6-31G** basis set; hence, the MRCISD energies may be less reliable than for the other test molecules.For instance, the MRCISD vertical excitation energies for the trans-and cisisomers of stilbene are 5.3 and 5.7 eV, respectively, and thus much higher than the experimental values (measured in nhexane) of 4.1 and 4.6 eV, respectively. 115By contrast, the OM2/MRCI, SSR, and SF vertical excitation energies for the two isomers agree reasonably well with the experimental figures.
Despite the use of the smaller 6-31G** basis set, the ab initio MRCISD MECI geometry agrees with the SSR and SF geometries to within 0.1 Å, see Table II.The OM2/MRCI method yields a larger RMSD for the tw-pyr MECI geometry, which is due to a somewhat greater torsion of the phenyl ring attached to the pyramidalized C atom (see Figure 3).At any rate, the BPs produced by all the methods are qualitatively similar, with projections r IJ on the order of 0.9 (see Table II).Similar to ethylene and styrene, the x 1 vector of stilbene corresponds to the pyramidalization mode and the x 2 vector to the ethylenic bond torsion mode, see Figure 7.

E. Penta-2,4-dieniminium cation, PSB3
The presence of a very strong electron-withdrawing center such as the cationic nitrogen in PSB3 leads to a tw-BLA MECI that corresponds to a crossing between the diradical and closed-shell CT configurations, which occurs in the course of the breaking of the central π -bond by torsion.PSB3 is perhaps one of the most studied molecular systems, for which the MECI point was obtained by a variety of methods, see Refs. 23, 74, 76, 121 and references cited therein.The tw-BLA character of its MECI has been unambiguously established, 22,23,121 and the methods employed in this work support this assignment.
As shown in Figure 8, the x 1 vector corresponds to the BLA mode, while the the MRCISD, SSR, and OM2/MRCI levels agree with each other to within 0.06-0.09Å (RMSD), whereas the SF geometry deviates from the others by a wider margin, ca.0.2-0.35Å.This discrepancy is caused by an excessive pyramidalization of one of the allylic carbon atoms in the SF calculations, which is caused by an overestimation of the relative stability of the ionic configuration resulting from the heterolytic breaking of the central π -bond.A similar overestimation has also been observed for the methyliminium cation, for which the SF method yields a pronounced pyramidalization of the carbon atom at the MECI geometry (see Figure 5).
Both DFT-based methods, SSR and SF, place the tw-BLA MECI ca.0.6 eV higher than MRCISD and OM2/MRCI (3.0 vs. 2.4 eV).This can again be explained by a greater preference of the DFT-based methods for the ionic configuration relative to the covalent one, which makes a somewhat stronger BLA distortion necessary to reach the S 0 /S 1 crossing point. 20,22,23,74 Aspointed out before, 74 the ionic configuration benefits from a more complete inclusion of dynamic electron correlation, and hence this configuration may become more important upon extending the basis set and the level of correlation treatment in the ab initio MRCI calculations.This may have the effect of increasing the current ab initio MECI energy (2.4 eV) and moving it closer to the current DFT-based values (3.0 eV).

F. Buta-1,3-diene
Butadiene is the simplest conjugated polyene.It has a transoid MECI with a -(CH 3 )-kink of the carbon backbone. 4,122 lectronically, this transoid MECI (see Figures 3 and 9) contains a kinked allylic fragment coupled to a terminal methylene radical.It was originally believed to provide the major non-adiabatic decay channel, 4, 122 but this conjecture has recently been challenged. 123he superimposed optimized geometries of the transoid MECI (see Figure 3 MECI geometries (see Table II) indicate a good agreement between the methods used presently.The largest deviation is observed for the MRCISD-OM2/MRCI pair with a RMSD of 0.10 Å.The BPs from the different methods agree very well with each other as confirmed by their projections r IJ in Table II and the visual representations in Figure 9.The x 1 vector describes a cyclization motion (with a σ -bond between the C 2 and C 4 atoms being formed or broken), while the x 2 vector corresponds to a BLA motion (stretching of the formal double bonds and contraction of the formal single bond and vice versa) coupled with a weak torsion of the C 1 methyl group.The transoid MECI lies between ca.4.8 eV and 5.6 eV above the S 0 minimum (see Table I).The computed vertical excitation energies of the trans-conformer (best estimate: 6.18 eV 111 ) cover an energy range of almost 1 eV; the ab initio MRCISD value is clearly too high, which is likely due to an insufficient treatment of dynamic electron correlation.

G. Ketene
Ketene is a paradigmatic example of an n/π conical intersection, which is associated with a crossing between the 1 A and 1 A states along the O-C-C bending mode. 117This crossing is linear in C s symmetry and becomes conical in C 1 symmetry, since upon symmetry lowering there emerges another direction along which the degeneracy of the electronic states can be lifted. 117As shown in Figure 9 (lower panel), all the methods used presently consistently give x 1 and x 2 vectors that describe an O-C-C bending mode and an outof-plane O-C-C-H mode, respectively.The MECI energy is computed to be close to 3.0 eV by the DFT-based methods, significantly higher than the OM2/MRCI values of 2.43 eV, see Table I.The latter is believed to be too low, in view of the fact that OM2/MRCI underestimates the vertical excitation energy of ketene by ca.0.7 eV (experimental value: 3.84 eV 114 ).

H. p-Hydroxybenzylideneimidazolinone anion, HBI
The lowest points on the CI seam of HBI are characterized by a torsion of the imidazole ring (MECI Im ) or the phenyl ring (MECI Ph ) accompanied by a mild pyramidalization of the methine bridge, 20,124 see Figures 3 and 10.These distortions lead to a breaking of the exocyclic π -bond, and the MECIs arise from a crossing between the respective diradical state (homolytic breaking) featuring a phenoxide anion radical (MECI Im ) or a imidazolinone anion radical (MECI Ph ) and the respective CT state (heterolytic breaking). 20As a result of the greater electron affinity of oxybenzyl as compared to imidazolinone, the MECI Im point lies ca.0.3-0.4eV lower in energy than the MECI Ph point. 20,124 e results for the FC point and the two MECI points of HBI are presented in Tables I and II and in Figures 3  and 10.The ab initio MRCISD calculations of anionic HBI could be carried out only with the smaller 6-31G** basis set, which causes a relatively large deviation of the vertical excitation energy (3.88 eV) at the FC point from the accurate quantum Monte-Carlo result of 3.06 eV. 116It is probably the insufficient treatment of dynamic electron correlation at the MRCISD/6-31G** level that is to blame for the overestimation of the excitation energy.The DFT-based methods, SSR and SF, accurately predict the vertical excitation energy in the range of 3.0-3.1 eV, and the OM2/MRCI value is only slightly lower at 2.90 eV.The MECI Im point is predicted by the two DFT-based methods to lie 0.1-0.3eV below the S 1 FC point, while the ab initio MRCISD and OM2/MRCI calculations give significantly larger gaps of 1.0 and 0.5 eV, respectively; however, all four methods agree that the MECI Ph point lies ca.0.1-0.3eV above the MECI Im point, see Table I.The computed MECI geometries agree fairly well with each other, except that OM2/MRCI gives a much smaller pyramidalization at the methine carbon atom.The agreement between the BPs from the three first-principles methods is good with r IJ projections between 0.8 and 1.0, whereas the projections involving OM2/MRCI are lower (0.5-0.7, Table II).For the two MECIs, the x 1 vector always features pyramidalization of the methine bridge coupled with BLA distortion of the two rings, while the x 2 vector mainly corresponds to torsion of the imidazolinone ring (MECI Im ) or of the phenoxy ring (MECI Ph ), see Figure 10.

V. CONCLUSIONS
The most important result of this work is that the currently applied methods yield a consistent description of the geometries and branching planes of the conical intersections over the whole range of organic molecules studied presently.The ab initio MRCISD method, the DFT-based approaches, SI-SA-REKS and SF-TDDFT, as well as the semiempirical OM2/MRCI method yield very similar geometries of the MECI points, with root-mean-square deviations between the different optimized MECI geometries of the order of 0.1 Å or even less.The branching planes produced by the different methods were compared using the BP-onto-BP projection r IJ (see Sec. II) which indicated very good agreement between the methods, with r IJ values varying in the range of 0.8-0.9.
A straightforward visual comparison between the BP vectors produced by two different methods is not feasible as the vectors may be rigidly rotated within the respective branching plane while leaving the latter unchanged. 8The similarity transformation of the BP vectors introduced in this work enables an alignment of the vectors obtained from dif-ferent calculations in such a way that their mutual spatial coincidence is maximized.The BP vectors obtained from SI-SA-REKS calculations refer to diabatic states, see Eq. ( 3), which allows for a chemically sensible interpretation of the nuclear motions that lead away from the CI region. 22Therefore, the BP vectors from other methods were pegged to the SI-SA-REKS vectors.Visual inspection of the transformed BP vectors does confirm that the methods used here yield very similar branching planes.They should thus all be suitable for studying the non-adiabatic dynamics of a wide variety of organic molecules.
Although the ab initio MRCISD calculations employed in this work cannot be regarded as converged with respect to the basis set size and the treatment of dynamic electron correlation, the good agreement between the results from MRCISD and the other methods gives us confidence that these other more approximate methods should perform well in modeling the non-adiabatic dynamics of photoexcited systems.It is particularly gratifying that a very simple semiempirical method, OM2/MRCI, which is several orders of magnitude less time consuming than the first-principles methods, yields MECI geometries and branching planes in close correspondence with the other methods.
Interestingly, at the DFT level, there is generally good agreement between the MECI geometries and BP vectors obtained from the spin-pure SI-SA-REKS method and the spincontaminated SF-TDDFT approach.Spin-contamination of SF-TDDFT S 0 and S 1 states varies for different molecules from very weak (0.01) to rather strong (0.93), as observed for MECIs in ethylene and styrene, respectively.Table XXVI of the supplementary material 52 lists the S 2 expectation values obtained for the SF-TDDFT S 0 and S 1 states at the MECI points.Although an a posteriori spin purification of the SF-TDDFT energies has been proposed, 76 its use lifts the degeneracy of the S 0 and S 1 energies at the CI point obtained in the spin-contaminated SF-TDDFT optimization (see Table XXVI of the supplementary material 52 for the energy gaps).Furthermore, as was recently observed by Gozem et al., 23 the use of the a posteriori spin purification destroys the correct dimensionality of the intersection seam, leading to linear rather than conical intersections.
Yet another word of caution concerns the widespread use of conventional spin-conserving linear-response TDDFT for modeling the ground-and excited-state PESs of photoactive systems.As conventional TDDFT does not give the correct dimensionality of the crossing seam, 23,77 it leads to certain artifacts in NAMD simulations of photoreactions. 45By contrast, the approximate methods used in this work yield the correct CI dimensionality.Their use is further encouraged by the present comparisons for a variety of organic molecules showing that they provide energies, geometries, and topologies of conical intersections that are in good agreement with each other, with ab initio MRCISD results, and with chemical intuition.

124122- 3
FIG. 1. Projection of a plane defined by a (a ,b ) vector pair onto an (a,b) plane.
FIG. 9. BP vectors of the transoid MECI of butadiene (upper panel) and the MECI of ketene (lower panel).

TABLE I .
Relative energies E (in eV) with respect to the lowest ground-state minimum.Energies at the Franck-Condon (FC) point are vertical excitation energies of the lowest singlet π → π * transition (best estimates given boldface in parentheses).Energies of optimized MECI structures refer to conical intersections between the S 1 and S 0 states (notation see text).

TABLE II .
Comparisons between the methods employed in this work: Root-mean-square deviations (RMSDs) between the MECI geometries and branching plane projections (see text).See legend to TableIfor the labels of the CIs.

TABLE II
a RMSD's are given as a matrix with the off-diagonal elements corresponding to the values for methods I and J.

TABLE III .
Relative energies E (in eV) with respect to the ground-state minimum: ab initio MRCISD results for ethylene with different active spaces and basis sets.