Description of Conical Intersections with Density Functional Methods.

Top Curr Chem (2015) DOI: 10.1007/128_2015_631 # Springer International Publishing Switzerland 2015

Description of Conical Intersections with Density Functional Methods Miquel Huix-Rotllant, Alexander Nikiforov, Walter Thiel, and Michael Filatov

Abstract Conical intersections are perhaps the most significant mechanistic features of chemical reactions occurring through excited states. By providing funnels for efficient non-adiabatic population transfer, conical intersections govern the branching ratio of products of such reactions, similar to what the transition states do for ground-state reactivity. In this regard, intersections between the ground and the lowest excited states play a special role, and the correct description of the potential energy surfaces in their vicinity is crucial for understanding the mechanism and dynamics of excited-state reactions. The methods of density functional theory, such as time-dependent density functional theory, are widely used to describe the excited states of large molecules. However, are these methods suitable for describing the conical intersections or do they lead to artifacts and, consequently, to erroneous description of reaction dynamics? Here we address the first part of this question and analyze the ability of several density functional approaches, including the linear-response time-dependent approach as well as the spin-flip and ensemble formalisms, to provide the correct description of conical intersections and the potential energy surfaces in their vicinity. It is demonstrated that the commonly used linear-response time-dependent theory does not yield a proper description of these features and that one should instead use alternative computational approaches. M. Huix-Rotllant Institute of Physical and Theoretical Chemistry, Goethe University Frankfurt, Max-von-LaueStr. 7, 60438 Frankfurt am Main, Germany A. Nikiforov and W. Thiel Max-Planck-Institut fu¨r Kohlenforschung, Kaiser-Wilhelm-Platz 1, 45470 Mu¨lheim an der Ruhr, Germany M. Filatov (*) Institut fu¨r Physikalische und Theoretische Chemie, Universita¨t Bonn, Beringstr. 4, 53115 Bonn, Germany e-mail: [email protected]

M. Huix-Rotllant et al.

Keywords conical intersections • ensemble DFT • excited states • spin-flip TD-DFT • TD-DFT

Contents 1 Introduction 2 Conical Intersections 2.1 Branching Plane Vectors 2.2 Double Bond Torsion and Conical Intersections 3 Computational Methods 3.1 REKS Method 3.2 Linear-Response Methods 4 Application of DFT Methods to Conical Intersections 4.1 Ethylene 4.2 Penta-2,4-Dieniminium Cation, PSB3 4.3 Ketene 5 Conclusions and Outlook References

1 Introduction The molecular potential energy surface (PES) is perhaps one of the most significant concepts in molecular physics which enables one to model molecular structure in terms of specific spatial arrangements of atoms (molecular geometries) [1]. Defining (meta)stable molecular conformations (local minima) and transient species (saddle points or transition states) in terms of the molecular electronic energy–molecular geometry relationship [2], the PES concept is a basis for theories of molecular structure and reactivity and serves as a starting point for even more approximate models, such as force-field molecular mechanics [3]. Rooted in the Born–Oppenheimer approximation [4], the PES concept is naturally valid in situations where the coupling between the motion of the nuclei and electrons is negligibly weak, that is, when the separation between the electronic states is much greater than the characteristic energy of nuclear motion [5]. This assumption breaks down in the vicinity of points where two (or more) PESs corresponding to distinct electronic states become degenerate, the surface crossing points [5–8]. Near these points, the non-adiabatic coupling between the electronic and nuclear degrees of freedom becomes decisive for the dynamics of transformations occurring in the excited states [9–12] as well as in the ground state [13, 14] of molecules. Especially interesting and at the same time especially challenging for their accurate theoretical description are situations of accidental (that is, not symmetry-imposed) crossing of PESs of states of the same spatial and spin symmetry, the conical intersections (CIs) [6, 15–19]. Such crossings occur in the subspace of M2 internal molecular degrees of freedom (for a nonlinear molecule, M ¼ 3 N 6, N – the number of atoms). The degeneracy between the electronic states is lifted along the two remaining directions, which span the branching plane

Description of Conical Intersections with Density Functional Methods

(BP) of the CI [6]. At a CI, the non-adiabatic coupling between the nuclear and electronic degrees of freedom becomes divergent [6, 9]; hence, the manifold of CI points, the CI seam, plays a decisive role for the dynamics of the electronic states by providing funnels for the efficient population transfer [9–12, 20]. As the BP defines all possible directions of exiting the strong non-adiabatic coupling region [6], it becomes a very important descriptor of the CI which, very similar to the transition vector in transition state theory [21, 22], determines the branching of possible products of the reaction occurring through the CI [23]. The correct theoretical description of CIs requires the use of computational methods capable of describing various electron correlation effects in the ground and the excited states with high and unbiased accuracy [10, 12, 24]. Naturally, since the first realization of the importance of CIs for the dynamics of excited states [9, 10], the computational investigation of CIs has been the realm of multi-reference (MR) methods of wavefunction theory (WFT), such as the complete active space self-consistent field (CASSCF) [25–28], the CAS-based second-order perturbation theory (CASPT2) [29], and the MR configuration interaction (MRCI) [30] methods. Although capable of delivering highly accurate results, these methods are computationally demanding and can be used routinely to calculate small molecular systems limited to just a handful of atoms. To go beyond this limitation and to match the challenges presented by the rapidly expanding use of ultrafast spectroscopy in photochemistry and photobiology [31, 32] one needs to employ computational methods capable of realistically describing large molecules. As the methods of density functional theory (DFT) [33, 34] offer a reasonably accurate description of electron correlation effects (which are vital for molecular bonding and structure) at a typical mean-field computational cost, these methods seem to represent a viable alternative to MR-WFT approaches. However, there is a caveat. Originally [33], DFT was formulated for ground states only, and it was deemed inappropriate to apply it to variational calculation of individual excited states [35–37]. A way around this problem is offered by the response formalism implemented in the currently widely used linear-response (LR) time-dependent DFT (TD-DFT) approach [38, 39]. LR-TD-DFT (or TD-DFT for brevity) enables one to obtain the excitation energies from the poles of the ground-state density-density response function [38, 39], thus enabling the use of ground-state theory for obtaining excited states. Although capable of describing crossings between excited states correctly, TD-DFT may experience difficulties with the S0/S1 Cls as the ground electronic state (S0) is the variationally optimized reference state for the response calculation and is thus described on a different footing than the excited state (S1), which is a response state. In particular, there is no coupling between the ground and excited states in most approximate LR-TD-DFT methods, which results in a wrong dimensionality of the S0/S1 crossing seam, M1 instead of M2 [40–43]. An interesting modification of the original LR-TD-DFT formalism is implemented in the spin-flip (SF) TD-DFT method which employs a variationally optimized high-spin open-shell reference state (e.g., triplet state) to access lower spin states (e.g., singlets) by allowing one-electron transitions with simultaneous


inversion of the spin (spin-flip transitions) [44–47]. In SF-TD-DFT, the ground singlet state and the excited singlet states are obtained on the same footing as response states. This enables some interaction between these states and reinstates the correct dimensionality of the S0/S1 conical intersection [40, 43, 48]. However, the downside of SF-TD-DFT is a substantial spin-contamination of the so-obtained singlet state, which may lead to the occurrence of mixed spin-symmetry states, thus complicating the identification of physically meaningful excited states [48]. An alternative to the response formalism, for obtaining excited states in the context of DFT, is offered by the ensemble DFT approaches [39, 49–52]. Ensemble DFT, formulated in the seminal works of Lieb [53] (ground-state ensembles) and of Gross, Oliveira, and Kohn [49] (ensembles of ground and excited states) enables one to obtain excitation energies from the variational calculation, as contrasted with the response states of LR-TD-DFT and SF-TD-DFT. A practical implementation of ensemble DFT in the form of a generally applicable computational scheme was achieved in the spin-restricted ensemble-referenced Kohn–Sham (REKS) method [54–56] and its state-averaged (SA) extensions, SA-REKS [57] and SI-SA-REKS [42, 58] (SI: state-interaction). In particular, the SI-SA-REKS method has proved its ability to describe properly the S0/S1 CIs by yielding the correct dimensionality of the CI, M2 [43, 58, 59]. Being a spin-restricted method which employs the same spatial orbitals for electrons with opposite spins, REKS is free from the spincontamination which infests the SF-TD-DFT description, and yields a spin and space (if present) symmetry adapted description of ground and excited states. The aforementioned computational approaches (LR-TD-DFT, SF-TD-DFT, and REKS) cover all the practically accessible implementations of DFT for describing the ground and excited states of molecules. In this chapter, the ability of these computational methods to describe CIs, and particularly S0/S1 CIs, of molecules is assessed in relation to probably the best description that the MR-WFT world can offer, namely MRCISD (MRCI with single and double excitations). When using approximate density functionals, LR-TD-DFT often yields good results for the energies of one-electron valence transitions and it seems tempting to believe that TD-DFT should be capable of yielding PES crossings as well. Although its inability to do so was clearly spelled out in the past [40, 41], the use of TD-DFT for modeling photodynamics and the underlying mechanistic features, such as CIs, seems to be gaining ground [60–73]. The consequences of using TD-DFT for photodynamics are reviewed in another chapter of this book;1 here, we focus on the ability (or, more precisely, the inability) of this formalism to provide a proper description of CI topography while yielding seemingly reasonable geometries and relative energies of the minimum-energy crossing points. As an alternative to TD-DFT, we investigate the performance of SF-TD-DFT and REKS methodologies which are known to yield the proper dimensionality of the CI seam [43]. Besides the molecular geometry at the minimum-energy CI

1 See the chapter “Surface Hopping Dynamics with DFT Excited States” by M. Barbatti and R. Crespo-Otero.


(MECI) and its relative energy, we focus on the BP of the intersection and compare the BP vectors obtained in the DFT calculations with the corresponding MRCISDderived vectors [59]. When benchmarking the DFT methods, we address a set of CIs of diverse topography occurring in unsaturated molecules during double-bond torsion; however, the conclusions drawn from these benchmarks should remain valid for other CIs as well, such as those in cyclic molecules or in nucleobases. The reason why we focus on relatively small molecules is that it would be extremely difficult to carry out MRCISD calculations with good-quality basis sets for larger molecules.

2 Conical Intersections CIs are manifolds of points at which there is a real crossing (degeneracy) between two (or more) adiabatic Born–Oppenheimer PESs of electronic states of the same spatial and spin symmetry. The fundamental conditions for the occurrence of such intersections have been known since the early days of quantum mechanics [15, 17]. For two adiabatic electronic states Ψm and Ψn, which can be represented by the solutions of a 2 2 secular problem

H mm H nm

H mn H nn

Cmm Cnm

Cmn Cnn

¼

Em 0

0 En

Cmm Cnm

Cmn ; Cnn

ð1Þ

formulated in terms of (in general, arbitrary) diabatic orthogonal states Φm and Φn as Ψm ¼ CmmΦm + CnmΦn and Ψn ¼ CmnΦm + CnnΦn, the crossing occurs whenever the two conditions H mm Hnn ¼ Em En ¼ 0;

ð2aÞ

H mn ¼ H nm ¼ 0;

ð2bÞ

are fulfilled [6, 9, 11, 17–19, 74]. In the above equations, Hmn are the matrix ^ |Φni, of the Hamiltonian H ^ in Born–Oppenheimer approxelements, Hmn ¼ hΦm| H imation evaluated with respect to the diabatic states. As was pointed out in the Introduction, these conditions can be fulfilled in the space of M2 internal molecular coordinates; thus, the CIs may occur in molecules with three or more atoms.

2.1

Branching Plane Vectors

The degeneracy of the two electronic states is lifted in the subspace of the two coordinates which can be found expanding the Hamiltonian in (1) around the point of degeneracy through the first order in internal nuclear coordinates Q,


H mm Hnm

Hmn H nn

I2 ∇Q H:δQ 0

1 1 ∇ ð H H Þ ∇ H Q mm nn Q nm B C þ @2 A δQ; 1 ∇Q H mn ∇Q ðH mm H nn Þ 2 g h δQ; ¼ I2 s δQ þ h g

ð3aÞ

ð3bÞ

where H ¼ ðH mm þ H nn Þ=2 is the average diagonal element, δQ is the nuclear displacement vector, ∇Q denotes differentiation with respect to nuclear coordinates, and I2 is a 2 2 unit matrix. The normalized vectors x1 ¼

g ; kgk

ð4aÞ

x2 ¼

h ; khk

ð4bÞ

span the branching plane of the CI which contains all the nuclear displacements lifting the degeneracy of the electronic states Ψm and Ψn. In the vicinity of a CI, the adiabatic PESs depend linearly on the nuclear displacements and have the topography of a double cone [6, 17]; hence the name. The definition of the BP vectors of a CI given in (3) and (4) is not unique [74]; the vectors spanning the same plane can be defined via the adiabatic states Ψm and Ψn, which leads to the following definition: 0

0

x1 ¼

g ; kg0 k 0

h x2 ¼ 0 ; h 0

0

g ¼ ∇Q ðEm En Þ;

ð5aÞ

0 ^ Ψn : h ¼ Ψm ∇Q H

ð5bÞ

0

0

Although the planes spanned by the (x1, x2) and (x1 , x2 ) vector pairs are the same, the individual vectors are not necessarily aligned with each other. Furthermore, as 0 0 was first demonstrated by Ruedenberg et al. [6], the vectors (x1 , x2 ) are rigidly rotated within the plane when traveling along a loop around a CI. This leads to a 0 0 certain arbitrariness in the directions of the (x1 , x2 ) vectors obtained from CI optimization carried out with finite numerical accuracy. As the crossing point at which the optimization stops may be infinitesimally close to the CI, but not exactly 0 0 at the CI, the BP vectors (x1 , x2 ) obtained in a series of CI optimizations utilizing


the same computational method but started from different initial conditions may not coincide with each other.2 0 0 When comparing BP vectors (x1 , x2 ) obtained using different computational schemes, the arbitrariness in their orientation within the BP can be bypassed by aligning the vector pairs such as to maximize the projection of one of the vectors, 0 0 e.g., x1 , obtained using computational method A onto the matching vector (x1 or x1) obtained using method B [59, 76]. As proposed recently [59], such an alignment can be carried out by a similarity transformation SRS1 that spans a 2D orthogonal rotation R and a shear transformation S. The latter transformation is necessary because the BP vectors are not, in general, orthogonal with respect to one another, 0 0 i.e., the inner products x1 x2 and (x1 · x2) are not zero. The shear transformation leaves the inner product invariant upon rotation within the plane. The same publication [59] also introduced a number of useful measures to compare the BPs obtained using different computational methods. Thus, projections pB xkA , k ¼ 1, 2 of the 0 vectors xk (or xk ) obtained using method A onto the BP obtained by method B and the projection rAB of a rectangle spanned by the vector pair (xA1 , xA2 ) onto the (xB1 , xB2 ) rectangle were introduced to quantify the similarity (or discrepancy) between the BP vectors produced by different computational methods [59]. These measures are used in the following, when comparing the computational methods addressed in this chapter. The use of the alignment procedure and the numerical measures described above for comparing BP vectors become especially useful when addressing computational schemes for which the BP vectors cannot be obtained using definitions (4) or (5), 0 such as the TD-DFT or SF-TD-DFT methods. For these methods, only the x1 vector can be determined explicitly by differentiation of the respective electronic energies. 0 The x2 vector is obtained during optimization for a CI using the branching space 0 update method of Maeda et al. [77] as an orthogonal complement to x1 in a plane 0 0 that iteratively converges to the BP of the CI. The (x1 , x2 ) vector pairs so obtained have been shown to approximate sufficiently accurately the true BPs of the optimized CIs [77]. However, the resulting vectors are strictly orthogonal, a property not shared by the explicitly calculated vectors. In the context of the REKS method, the BP vectors are defined by (4), which uses diabatic states and, therefore, yields a unique orientation of the vectors. This feature becomes especially important when analyzing the Cls and their BPs in chemical terms.

2

As shown by Yarkony [74, 75], a rotation of the crossing states which orthogonalizes the BP vectors brings the vectors to a unique orientation, especially when symmetry is present. Application of such a prescription, however, modifies the BP vectors, while leaving the BP unchanged.


2.2

Double Bond Torsion and Conical Intersections

In this chapter, we focus on the Cls arising during double-bond torsion in unsaturated molecules; thus, let us look more closely at this situation. The origin of Cls for double-bond torsion can be traced back to a crossing between the electronic states caused by homolytic and heterolytic breaking of the π-component of the double bond [78]. The homolytic bond breaking results in a diradicaloid electronic configuration, whereas the heterolytic bond breaking leads to an ionic (or zwitterionic) electronic configuration. According to the sign-change theorem of Longuet-Higgins [19], a CI should be present inside a loop connecting the conformations that correspond to the two bond-breaking mechanisms [12, 79–81]; see Fig. 1. Let us now assume a nuclear movement in the direction of the x1 vector while keeping the interstate coupling element Hmn at zero. When passing through the CI, the S0 wavefunction experiences a sudden switch from ionic to diradical (or vice versa) and the S1 wavefunction does precisely the opposite. It is therefore natural to associate with the x1 vector a direction that corresponds to the transition between the uncoupled ionic and diradical states; it is these two electronic states that are included into the ensemble averaging in the SI-SA-REKS method; see Sect. 3.1. Conversely, a displacement along the x2 vector while keeping the energy difference Hmm–Hnn at zero should correspond to increasing (decreasing) the coupling between the states; hence, this motion should contain the torsion about the double bond axis. Indeed, when the two fragments connected by the double bond attain an

Fig. 1 Schematic representation of conical intersection for double bond isomerization. Reproduced with permission from [81]. Copyright © (2014) American Chemical Society


approximately orthogonal orientation during the torsion, the interaction between the fragment wavefunctions vanishes and this results in a vanishing Hmn element. Thus, the gradient of the interstate coupling element in the SI-SA-REKS method (see Sect. 3.1) can be conveniently associated with the x2 vector. From this brief discussion, the advantage of the specific orientation of the BP vectors as in Fig. 1 is that it offers a simple chemical interpretation of the CI and the factors influencing its occurrence. In particular, the relative electronegativity of the fragments connected by the double bond and, consequently, the relative preference for one of the bond breaking mechanisms define the molecular geometry at a CI [58, 81]. For example, in alkenes, the homolytic π-bond breaking is energetically preferred over the heterolytic bond breaking, and this results in a pronounced pyramidalization at a CI which is needed to stabilize the ionic structure and (along with the torsion) to reach the point of surface crossing. By contrast, in organic molecules with strongly electron-withdrawing (or electron-donating) functional groups, the two bond breaking mechanisms may become nearly isoenergetic and a CI can be reached without requiring pyramidalization; in this case, the molecular geometry at the CI corresponds to double bond torsion combined with stretching/compression of the other single and double bonds, hence a bond length alternation (BLA) distortion [81]. Because the BP vectors obtained in the SI-SAREKS calculations have a unique orientation and a transparent chemical interpretation, the vectors obtained by other methods, TD-DFT, SF-TD-DFT, and MRCISD, are aligned with them by using the recently introduced similarity transformation [59].

3 Computational Methods In the following, we give a brief overview of the computational methods used to study the CIs in this chapter. A complete description of these computational schemes can be found in other chapters of this book, and so we recapitulate only those features most relevant for PES crossings. The interested reader is advised to consult the other chapters in this book and the literature references cited therein for more detail on the computational methodologies.

3.1

REKS Method

The REKS method [54–56] employs the ideas behind ensemble DFT to describe the non-dynamic correlation in the ground and excited states of molecules and to access the excited states through the application of the variational principle [42, 57, 58]. Perhaps the most important realization in ensemble DFT that was rigorously proved by Lieb [53] and computationally verified by Schipper et al. [82, 83] and by Morrison [84] is that the density and the ground-state energy of a strongly


correlated fermionic system is to be exactly represented by a weighted sum (ensemble) of the densities (and energies) of several electronic configurations (ensemble components). In the REKS method, the ensemble representation of the non-interacting Kohn–Sham reference system is used to describe the non-dynamic correlation arising because of near degeneracy of several electronic configurations; such as in situations with dissociating chemical bonds, near transition states of symmetry-forbidden reactions, biradicaloid species, etc. The ensemble representation leads to fractional occupation numbers (FONs) of several frontier Kohn–Sham orbitals, which are obtained simultaneously with the orbitals from the variational optimization of the total ground-state energy. The excited states are accessed within the REKS methodology by applying the ensemble variational principle proved by Gross et al. [49]: M M M X X X ^ ΦK λK EK ; 0 λK 1; λK ΦK H λK ¼ 1 ; K¼1

K¼1

ð6Þ

K¼1

which states that the energy of an ensemble of trial wavefunctions ΦK representing several lowest states of the many-body Hamiltonian Hˆ is always bounded from below by the weighted sum of the exact eigenenergies of this Hamiltonian taken with the same (positive definite) weighting factors λK. Restricting the ensemble averaging to two states, S0 and S1, the ground and the lowest excited states of an atom or a molecule, one obtains the excitation energy from the variational optimization of an ensemble Eω ¼ ð1 ωÞE0 þ ωE1 ;

ð7Þ

of the two states by taking the energy difference ΔE ¼ E1 E0 ¼

E ω E0 ; ω

ð8Þ

where 0 ω 1 is a fixed (that is, not variationally optimized) weighting factor. For a system with two strongly correlated electrons in two orbitals (e.g., a biradical or a dissociating bond), the ground state as described by the REKS(2,2) method corresponds to a two-configurational model wavefunction [85], rffiffiffiffiffi rffiffiffiffiffi na nb Φ0 ¼ . . . ϕa ϕa i . . . ϕb ϕb i; 2 2

ð9Þ

where ϕa and ϕb are the fractionally occupied frontier orbitals, na and nb are the respective FONs, and the unbarred and barred orbitals are occupied with α- and βspin electrons, respectively. When the two active orbitals belong to two different irreducible representations of the molecular symmetry group (e.g., a homosymmetric biradical, H2 with stretched bond, double bond torsion in C2H4,


etc.), the lowest singlet excited state arising from a one-electron excitation in the space of the two orbitals ϕa and ϕb can be approximated by an open-shell singlet (OSS) wavefunction [85], 1 1 Φ1 ¼ pffiffiffi . . . ϕa ϕb i þ pffiffiffi . . . ϕb ϕa i:: 2 2

ð10Þ

The OSS state can be described by the spin-restricted open-shell Kohn–Sham (ROKS) method [86, 87]. Combining the two energies, EREKS(2,2) and EROKS in an ensemble as in (7) and variationally optimizing the KS orbitals and the FONs of the active orbitals in REKS, one obtains the state-averaged energy Eω from which the excitation energy is obtained by (8). In the SA-REKS method described, the same set of KS orbitals is used to construct the two energies, EREKS(2,2) and EROKS. Typically, the weighting factor ω is set to 1/2, which corresponds to an equi-ensemble of the S0 and S1 states. The SA-REKS method [57] treats S0 and S1 as uncoupled states, e.g., as in a homosymmetric biradical where the interaction between the states is prevented by symmetry. When the two states belong to the same symmetry species, e.g., as in a heterosymmetric biradical, the states approximated by the wavefunctions in (9) and (10) interact with each other and this interaction should be taken into account when calculating the energies of the individual ensemble components in (7) and (8). Within the REKS formalism, the uncoupled S0 and S1 states can be obtained by diagonalizing a 2 2 secular matrix

EREKSð2;2Þ H 01

H01 ; EROKS

ð11Þ

where the coupling matrix element H01: H 01 ¼

pffiffiffiffiffi ^ pffiffiffiffiffi ^ pffiffiffiffiffi pffiffiffiffiffi na ϕa na F a ϕa nb ϕa nb F b ϕb ¼ ð na nb Þεab ;

ð12Þ

is obtained by the application of Slater–Condon rules and the variational condition for the SA-REKS orbitals. In (12), F^ a and F^ b are the Fock operators for the openshell orbitals and εab is the off-diagonal Lagrange multiplier in the open-shell Lagrangian [54]. Equations (11) and (12) constitute the state-interaction SI-SAREKS method [42, 58]. Provided that the weighting factor ω in (7) is set to 1/2, the state-averaged energy Eω in the SA-REKS and SI-SA-REKS methods remains the same and the SA-REKS orbitals can be used in the SI-SA-REKS method. Note that (11) and (12) can also be obtained by using the adiabatic connection argument for an ensemble of two states [52].3

3 See also the chapter “Ensemble DFT approach to excited states of strongly correlated molecular systems” by M. Filatov.


In the framework of the SI-SA-REKS method, a CI occurs when the two energies EREKS(2,2) and EROKS become equal and the coupling matrix element H01 vanishes. Hence, the BP vectors x1 and x2 of the CI can be obtained by the differentiation of the respective energy differences and matrix elements as in (4), where the g and h vectors are given by (13) [59, 81]: g ¼ ∇Q EREKSð2;2Þ EROKS ;

ð13aÞ

h ¼ ∇Q H 01 :

ð13bÞ

The BP vectors as defined in the SI-SA-REKS method are in accord with the chemical interpretation presented in Sect. 2 [81]. Indeed, the x1 (or g) vector points in the direction of the difference of the gradients of the energies of the noninteracting states S0 and S1, exactly as it should be when moving along the x1 vector in Fig. 1. The same is true for the x2 vector as well. As the SI-SA-REKS energy formula can be obtained by the adiabatic connection argument, and hence within the domain of DFT, the BP vectors given in (13) can be regarded as approximations to the true DFT BP vectors, had these vectors been known from the exact theory. Thus, when using the SI-SA-REKS method, the excited states are obtained from a variational calculation in strict correspondence with the ensemble variational principle of DFT [49]. The variational calculation of excited states is also implemented in constricted variational DFT (CV-DFT) [88, 89] reviewed in another chapter of this book. The use of the SI-SA-REKS method offers markedly improved accuracy when describing certain types of excitations, particularly excitations of extended π-conjugated molecular systems [90], charge-transfer excitations in donor-acceptor systems [91], and excitations of strongly correlated groundstate systems (biradicals, molecules with broken bonds, etc.) [90]. With regard to CIs, the SI-SA-REKS method yields the correct dimensionality (M 2) of the CI seam and the correct shape of the S0 and S1 PESs in its vicinity [43, 59]. This follows not only from the numerical tests presented below in this chapter but also from the general theoretical argument that the coupling between the S0 and S1 states is consistently taken into account in the SI-SA-REKS method.

3.2

Linear-Response Methods

In linear-response methods with local potentials, the excitation energies are obtained from the residues of the exact Fourier-transformed density-density response function χ [38], ð 0 0 0 δρðr; ωÞ ¼ χ r; r ; ω δvext r ; ω dr ; ð14Þ


where vext(r, ω) corresponds to the Fourier-transformed time-dependent external potential. Separating the time-dependent Hamiltonian into time-independent (zeroth-order) and time-dependent (first-order) parts, one arrives, after a few steps, at the linear-response equation 0 0 χ LR r; r ; ω ¼ χ s r; r ; ω ðð 0 þ χ s ðr; r1 ; ωÞ f Hxc ðr1 ; r2 ; ωÞχ LR r2 ; r ; ω dr1 dr2 ;

ð15Þ

where χ LR(r,r0 ,ω) is the LR function, χ s(r,r0 ,ω) is the zeroth-order response function, and fHxc(r1,r2,ω) ¼ |r1 – r2|1 + fxc(r1,r2,ω) is the Hartree plus the exchangecorrelation (xc) kernel. All LR approaches extract the poles of some form of (15), differing mainly in the choice of the zeroth-order density used to construct the non-interacting response function and the xc kernel.

3.2.1

LR-TD-DFT

In LR-TD-DFT, the poles of the LR function (15) are cast in matrix form

AðωÞ B * ð ωÞ

BðωÞ A* ðωÞ

1 X ¼ω Y 0

0 1

X ; Y

ð16Þ

where the matrices A(ω) and B(ω) are defined as KS ½AðωÞai, b j ¼ εKS δi j δab þ ai f Hxc ðωÞb j ; a εi ½BðωÞai, b j ¼ ia f Hxc ðωÞb j :

ð17aÞ ð17bÞ

In (17), i, j, . . ., a, b,. . . and p, q,. . . are the indices for occupied, virtual, and general (occupied or virtual) spin-orbitals, and εKS are the eigenvalues of the KS p Hamiltonian. Adiabatic LR-TD-DFT based on a KS initial density employs the XC kernel taken in the adiabatic approximation (i.e., implying locality of the exchangecorrelation (XC) kernel in the time domain), 0 δ2 Exc ½ρ AA f xc : r; r ¼ δρðrÞδρðr0 Þ

ð18Þ

The fact that the xc kernel is frequency independent makes it impossible to account properly for doubly excited states within the adiabatic approximation [92, 93]. These reasons are believed to be the major causes of the failure of LRTD-DFT methods to describe correctly the low-lying excited states in


highly-conjugated systems [94]. Recent advances in LR-TD-DFT beyond the adiabatic approximation are surveyed in another chapter of this book.4 When using the conventional closed-shell KS reference state, the LR-TD-DFT description in the adiabatic approximation breaks down near an S0/S1 conical intersection. Indeed, the assumption of pure-state v-representability employed in the traditional KS theory fails near the surface crossing points where the two electronic states become (nearly) degenerate. This is usually reflected in severe convergence problems of the KS self-consistent field iterations, which can be alleviated by either allowing holes below the Fermi level [82] or using a more general density Ansatz, such as in ensemble DFT [53, 82]. Furthermore, the interaction between the ground and the excited states is missing in adiabatic LR-TD-DFT, which results in a linear crossing instead of a conical intersection [40]. Surprisingly, it was observed by Levine et al. [40] and by Cordova et al. [95] that, even though the PESs are defective, the geometries and branching planes of minimal energy crossing points may look reasonable in comparison with other, theoretically justified, methods. The topography of conical intersections between different excited states is expected to be correct even at the adiabatic LR-TD-DFT level; however because of the missing effect of doubly-excited configurations its geometry and energy level may be incorrect [40].

3.2.2

Spin-Flip TD-DFT

Several schemes have been developed that attempt to ameliorate some of the problems arising in the adiabatic LR-TD-DFT description of S0/S1 conical intersections. Among them, one of the most successful is the spin-flip time-dependent DFT, SF-TD-DFT. In SF-TD-DFT, electronic configurations arising from the excitation operators with ΔMs ¼ 1 are coupled, unlike the usual formulation of TD-DFT in which only spin-preserving excitation operators with ΔMs ¼ 0 are allowed. Several formulations of SF-TD-DFT appear in the literature [44– 47]. Ziegler and Wang’s formulation of SF-TD-DFT [45] relies on the non-collinear spin DFT framework, which operates with the KS spinors, rather than spin-free orbitals. In this formulation, the xc kernel is given by 0 ν α ðrÞ ν β ðrÞ xc f SF ; ¼ xc xc r; r ρα ðr0 Þ ρβ ðr0 Þ

ð19Þ

where νσxc and ρσ (r0 ) are the xc potential and the electronic density of spin σ, respectively. The configurations arising from the spin-flip αβ and βα excitations are decoupled. Starting from an open-shell triplet reference in a two-electron

4 See the chapter “Current status and recent developments in linear response time-dependent density-functional theory” by Mark E. Casida and Miquel Huix-Rotllant.


two-orbital model, this gives rise to two closed-shell configurations (a double excitation and a ground state) and two open-shell configurations which can be coupled to a singlet or a triplet configuration (see Fig. 2). Beyond this model, openshell spin-uncompensated configurations are present so that spin-contaminated excitation energies are obtained, which can be partially purified by an a posteriori correction scheme [96]. The SF-TD-DFT approach has several advantages over adiabatic LR-TD-DFT, especially when describing a conical intersection region involving the ground state. First, the reference state is a triplet, which satisfies more easily the KS condition of non-interacting pure state v-representability. Second, SF-TD-DFT includes some extra double-excitation character in the excited states because of the configuration arising from the a{β iα excitation. Thus, ground and excited states are coupled in SF-TD-DFT, and this yields 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 [40]. SF-TD-DFT is able to introduce the ground and excited state coupling correctly, albeit only within the restricted configuration space represented in Fig. 2. In general, most CIs are well represented by SF-TD-DFT, but more complicated intersections are described only approximately [48]. One of the most important problems is that SF-TD-DFT states are frequently spin-mixed. In order to correct this problem, an a posteriori spin purification scheme has been proposed [96]. However, when using this scheme, the correct dimensionality of the CI seam is lost [43]. Within the two LR approaches, TD-DFT and SF-TD-DFT, CIs can be located using the branching space update method of Maeda et al. [77] which employs a 0 projected gradient algorithm for optimizing the branching plane. The x1 -vector is 0 calculated by (5) and the x2 -vector is approximated by an iterative update scheme. 0 0 Note that in this approach the x1 and x2 vectors are strictly orthogonal, although the exact BP vectors need not necessarily have this property [6].

Fig. 2 Spin-flip excitations from an open-shell triplet reference


4 Application of DFT Methods to Conical Intersections In this section we discuss the application of the computational methods described above to study Cls in a number of organic molecules and organic and biological chromophores [59]. The geometries of the MECIs optimized using the SI-SAREKS (denoted for brevity as SSR), SF-TD-DFT (abbreviated to SF), and LRTD-DFT (abbreviated to TD) methods are superimposed with the MRCISD geometries in Fig. 3. All the calculations reported here employed the 6-31+G** basis set and the BH&HLYP density functional (DFT calculations), with the exception of stilbene and anionic HBI, for which the MRCISD calculations employed a smaller 6-31G** basis set. The results of the MRCISD, SSR, and SF calculations are taken from [59] and the TD calculations were carried out using the same settings as in the SF calculations. In particular, the Tamm–Dancoff approximation (TDA) was used in the SF and TD calculations and both methods applied the algorithm by Maeda et al. [77] to locate the MECI points. The majority of CIs arising in the molecules in Fig. 3 originate from a crossing between diradical and ionic electronic configurations along the reaction coordinate which corresponds to double bond torsion [12, 78, 79]. As discussed in Sect. 2, these CIs can be classified as twist-pyramidalization CIs (tw-pyr, for brevity), for

Fig. 3 Superimposed geometries at the respective MECI points optimized using MRCISD/6-31 + G** (standard colors), SSR-BH&HLYP/6-31 + G** (solid blue), SF-TD-DFT-BH&HLYP/631 + G** (solid green), LR-TD-DFT-BH&HLYP/6-31 + G** (solid red). MRCISD calculations employed the 6-31G** basis set for stilbene and anionic HBI. The MRCISD, SSR, and SF geometries are taken from [59]


situations where the homolytic breaking of the π-bond is energetically preferred, and as twist-bond_length_alternation (tw-BLA) CIs, for cases when the two π-bond breaking mechanisms, homolytic and heterolytic, either become nearly isoenergetic or the latter mechanism is preferred [81]. An accurate description of this type of CI requires a balanced description of the relative stability of diradical and ionic electronic configurations; thus the ability of computational methods to deliver such a description is probed by this type of CI. Another type of CI shown in Fig. 3 can be described as n/π* CIs that originate from the crossing between an electron configuration with a doubly occupied lone pair orbital (n2π*0) and a singly excited configuration (n1π*1) [97]. These CIs occur in ketene, ethylene (ethylidene or methylcarbene CI), and methyliminium cation (methylimine CI) as shown in Fig. 3. This type of CI probes the accuracy of description of lone-pair excitations by the tested theoretical methods. The set of CIs shown in Fig. 3 misses some other types, in particular those found in cyclic molecules, such as nucleobases; however, good-quality reference MRCISD data for these CIs would be extremely difficult to obtain because of the large size of these molecules. We therefore omit these CIs from comparisons and focus instead on the ability of the DFT computational schemes to reproduce the molecular geometry and relative energy of the MECI points and the shape of the S0 and S1 PESs in their vicinity. To this end, we compare the BP vectors obtained using the DFT approaches with the respective MRCISD vectors, and we investigate the S0 and S1 PESs around the optimized MECI points by scanning the PESs (1) in the direction of the BP vectors and (2) along a closed path around the MECI point [98]. The latter scans enable us to evaluate the ability of the respective theoretical methods to describe the correct dimensionality of the CI seam [43, 98]. As is evident from Table 1, all the DFT methods are capable of describing the vertical excitations at the Franck–Condon (FC) point sufficiently accurately as compared to the best estimates of these energies. The accuracy of the MRCISD description is sometimes inferior to the DFT methods, probably because of certain restrictions with regard to the primary active space and the basis set employed [59]. The TD method tends to overestimate the vertical excitation energies as compared to the other two DFT methods. The relative energies of the MECI points obtained in the DFT calculations are in semi-quantitative agreement with the MRCISD energies; deviations of ca. 0.5 eV or slightly more (e.g., for ketene or PSB3) are observed in Table 1. However, the MRCISD energies in Table 1 cannot be regarded as highly accurate benchmark reference data because of the aforementioned basis set and active space restrictions; these data illustrate what is currently achievable at the computationally affordable MR-WFT level of calculation [59]. The accuracy of the optimized molecular geometries at the MECI points is characterized in Fig. 3 in graphical form and in Table 2 in terms of root-meansquare-deviations (RMSDs) between the optimized Cartesian coordinates. With the exception of the PSB3 cation, for which the SF and TD methods predict a much too strong pyramidalization of the C3 atom, there are no conspicuous differences between the DFT and MRCISD geometries. This is reflected in the DFT vs

cis-FC pointm tw-pyr MECI

l 0 ΔEScistrans

cis-FC pointm,o transoid MECI Stilbene trans-FC pointk

l 0 ΔEScistrans

Butadiene trans-FC pointn

tw-pyr CI2j

Geometry Ethylene FC pointe tw-pyr MECI eth MECIf Styrene FC point tw-pyr CI1h

(4.6) [103]

(4.1) [103]

(6.18) [99]

5.76 4.68

(4.88) [100]

4.54 4.26

4.18 0.18

(5.35)q (0.13)

(5.74) (4.50)

5.54 5.58

5.76 0.14

4.97

5.10 4.29

7.55 4.96 3.96

SSRb

6.68 4.89

6.89 0.07

5.28

8.11 4.79 4.69

(7.80) [99]

MRCIa

4.49 4.20

4.35 0.35

5.96 5.19

5.94 0.11

4.83

5.29 4.01

7.70 4.82 4.52

SFc

4.86 n.a.i

4.42 0.33

6.05 5.40

6.43 0.14

4.72

5.25 n.a.i

7.53 4.73 4.26

TDd

tw-pyr MECIphr

Anionic HBI FC point tw-pyr MECIImq

l 0 ΔEScistrans cis-FC pointm tw-BLA MECI Ketene FC point MECI

(3.06) [104]

(3.84) [102]

(4.20) [101]

Geometry Methyliminium FC point tw-BLA MECI meth MECIg Penta-2,4-dieniminium cation, PSB3 trans-FC pointk

(3.20)

(3.88)p (2.87)

4.01 2.43

4.24 2.38

4.47 0.16

9.52 3.82 5.63

MRCI

3.24

3.00 2.97

3.93 3.00

4.08 3.00

4.18 0.15

8.62 3.83 4.54

SSR

3.20

3.14 2.84

4.01 2.82

4.28 3.08

4.43 0.13

9.22 4.27 5.46

SF

2.81

3.56 2.60

4.06 3.27

4.93 2.82

4.97 0.14

8.32 4.14 5.22

TD

Table 1 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 S1 and S0 states. The MRCISD, SSR, and SF data are taken from [59]


b

MRCISD/6-31+G** method SI-SA-REKS-BH&HLYP/6-31+G** method c SF-TD-DFT-BH&HLYP/6-31+G** method d LR-TD-DFT-BH&NLYP/6-31+G** method e Franck–Condon point corresponding to the So minimum f Ethylidene MECI g Methylimine MECI h Twisted-pyramidalized MECI1 i LR-TD-DFT energy is not available because of lack of convergence j Twisted-pyramidalized MECI2 k trans-conformation l Energy difference between cis- and trans-conformations in the ground state m cis-conformation n trans-conformation, transition to the bright Bu state o Transition to the B1 state p 6-31G** basis set q Twisted-pyramidalized MECI with twisted imidazole ring r Twisted-pyramidalized MECI with twisted phenyl ring

a


MRCI 0.0000 0.0609 0.0698 0.0789

SSR 0.0609 0.0000 0.0712 0.0504

SF 0.0698 0.0712 0.0000 0.0365

TD 0.0658 0.0420 0.0304 0.0000

BP projectionb MRCI SSR MRCI 1.0000 SSR 0.7845 SF 0.8996 TD 0.7378

SF 0.8059 1.0000 0.8064 0.7677

TD 0.9027 0.8089 1.0000 0.8111

0.6307 0.6459 0.6721 1.0000

b

RMSDs are given as a matrix with the off-diagonal elements corresponding to the values for methods A and B BP projections are given as a matrix with the off-diagonal elements rAB representing the projection of the BP obtained using method A onto the BP obtained using method B

a

MRCI SSR SF TD

RMSD (Å)a

Table 2 Root-mean-square deviations (RMSDs) between the MECI geometries and branching plane projections averaged over all molecules in Fig. 3. The MRCISD, SSR, and SF data are taken from [59]



MRCISD RMSDs in Table 2, which vary in the range between 0.061 Å (SSR) and 0.079 Å (TD). An even better agreement is seen between the geometries obtained by the different DFT methods, especially for the SF and TD methods (0.030 Å). The BP projections in Table 2 indicate that all methods produce nearly the same BP vectors at the MECI geometries. The very good agreement between the DFT methods and MRCISD seems to indicate that these much simpler methods can be confidently used in lieu of MRCISD for locating CIs and analyzing the S0 and S1 PESs in their vicinity. However, there is a caveat: despite the seemingly good agreement with the high-level MR-WFT approach, not all the DFT methods produce the correct dimensionality of the CI seam. In the following, we focus on three molecules and their respective CIs, ethylene and its tw-pyr MECI, cationic PSB3 and its tw-BLA MECI, as well as ketene molecule and its n/π* MECI, as representatives of the specific types of CIs discussed above.

4.1

Ethylene

On the S0/S1 CI seam in ethylene there are two distinct MECI points that correspond to a tw-pyr MECI (the major photorelaxation funnel) and to a MECI point with partial transfer of a hydrogen atom; see the tw-pyr and ethylidene MECIs of C2H4 in Fig. 3 [105–108]. The tw-pyr MECI originates from a crossing between the PESs for homolytic and heterolytic π-bond breaking along the double bond torsion mode. As the heterolytic π-bond breaking is strongly disfavored in ethylene, the tw-pyr MECI features a strong pyramidalization of one of the carbon atoms to stabilize the ionic electronic configuration. All the DFT methods employed here predict the vertical excitation energy at the FC point in a good agreement with the best theoretical estimate of 7.8 eV. The energy level of the tw-pyr MECI is much lower than the FC point. The DFT methods yield tw-pyr MECI energies varying between 4.96 eV (SSR) and 4.73 eV (TD), in good agreement with the MRCISD value 4.79 eV; see Table 1. The BP vectors predicted by the DFT methods are in very good agreement with the MRCISD BP vectors, as can be seen from a visual comparison of the vectors in Fig. 4 and from the BP projections (see Sect. 2) which lie in the range 0.99–0.97 in each case. From these results it seems that all the DFT methods are capable of describing the tw-pyr MECI in ethylene correctly. However, a closer look at the S0 and S1 PESs near the MECI point reveals that the surface crossing predicted by adiabatic LRTD-DFT is not actually a conical intersection but is rather a linear crossing. As seen in Fig. 5, the S1 PES produced by adiabatic LR-TD-DFT falls below the S0 PES when moving along the x1 direction.5 Furthermore, the shape of the S0 and S1 PESs

5

The conventional KS DFT/LR-TD-DFT calculations experienced severe convergence problems in the vicinity of the MECI point, which are reflected in the shape of the S1 PES in the upper right panel of Fig. 5.

M. Huix-Rotllant et al. Fig. 4 BP vectors of twisted-pyramidalized MECI of ethylene calculated using MRCISD, SSR-BH&HLYP, SF-TDBH&HLYP, and LR-TDBH&HLYP methods

along the x2 direction as predicted by adiabatic LR-TD-DFT is inconsistent with the SSR and SF PESs. It appears as though the S0 state is unaware of the proximity to the S1 state. The SSR and SF methods yield some interaction between these states which is reflected in the shapes of the S0 and S1 PESs; the latter appear to be (approximately) symmetrically split around some median level, whereas no such (quasi) symmetric splitting can be seen in the adiabatic LR-TD-DFT PESs; see the lower right panel of Fig. 5. Thus, adiabatic LR-TD-DFT predicts a linear S0/S1 PES crossing, as expected from the vanishing off-diagonal matrix element in (1); hence, only one degree of freedom is sufficient to make the PESs cross, which leads to the M 1 dimensionality of the crossing seam in adiabatic LR-TD-DFT [40, 41, 43]. These observations are corroborated by the plot in Fig. 6, which shows the S0–S1 energy difference along a loop with a radius of 0.002 Å around the MECI point. The SSR and SF energy differences always remain finite which means that there is a finite gap at all geometries which lift the degeneracy of the S0 and S1 states at the MECI point. By contrast, the adiabatic LR-TD-DFT energy difference curve crosses the zero level twice when going around the loop. It is also noteworthy that the shape of the energy difference curve as produced by adiabatic LR-TD-DFT is markedly different from those of the other two methods. The angle θ that parameterizes the loop is taken from the direction of the x1 vector. This implies that the coupling between the S0 and S1 states should vanish at the points θ ¼ 0 and


Fig. 5 S0 (blue) and S1 (red) PES profiles near the tw-pyr MECI of ethylene. Upper panels – profiles along the x1 BP vector, lower panels – along the x2 BP vector. Leftmost panels – profiles obtained using the SSR-BH&HLYP/6-31 + G** method, middle panels – obtained by the SF-TDBH&HLYP/6-31 + G** method, rightmost panels – using the TD-BH&HLYP/6-31 + G** method. For each computational method, the MECI energy level is taken as the origin of the energy scale (in kcal/mol). The variable x corresponds to displacement (in Å) from the MECI geometry in the direction of the x1 or x2 vector, respectively

Fig. 6 The S0–S1 energy difference (in kcal/mol) along a loop around the MECI point of ethylene. The solid line corresponds to the SSR-BH&HLYP/6-31 + G** method, the dashed line to the SF-TD-BH&HLYP/6-31 + G** method, and the dotted line to the TD-BH&HLYP/6-31 + G** method. The loop radius is 0.002 Å. For each method, the angle θ (in deg) is taken from the direction of the x1 vector


θ ¼ 180 ; note that the off-diagonal element Hmn in (1) is supposed to zero out when moving along the x1 direction. The SSR and SF methods predict the maximum S0– S1 splitting in the x1 direction, whereas adiabatic LR-TD-DFT yields the maximum splitting approximately in the x2 direction (θ 0 ). As found by Gozem et al. [43, 98] for the example of PSB3, the multi-state CASPT2 and MRCISD methods yield energy difference curves around a MECI which are similar in shape to those predicted by SSR and SF and dissimilar to TD-DFT. Therefore, in spite of its seemingly good predictions for the geometry and energy of the tw-pyr MECI of ethylene, TD-DFT fails to give the correct conical topography of the crossing and yields the wrong shape of the two PESs in the vicinity of the crossing point. As discussed by Barbatti and Crespo-Otero (see Footnote 1) this has severe consequences for the description of the dynamics near the crossing points, rendering TD-DFT practically useless for modeling non-adiabatic population transfer through PES crossings.

4.2

Penta-2,4-Dieniminium Cation, PSB3

In the protonated Schiff base PSB3 (penta-2,4-dieniminium), which is often used as the simplest model of the Schiff base retinal chromophore, the presence of a strong electron-withdrawing substituent, the cationic nitrogen, leads to the occurrence of a tw-BLA MECI which originates from a crossing between the diradical and ionic electronic configurations obtained by homolytic and heterolytic breaking of the central π-bond along the torsion mode [42, 43, 81, 109, 110]. The ionic configuration is stabilized by the presence of the electron-withdrawing group and no pyramidalization is needed to reach the crossing point with the diradical configuration. The SSR and SF methods predict the vertical excitation energy at the FC point in good agreement with the available theoretical best estimate, the quantum Monte– Carlo excitation energy of 4.2 eV for the cis-conformation. The TD-DFT FC excitation energy is in error by more than 0.7 eV, which is typical of this method when applied to extended π-conjugated systems, such as, e.g., cyanine dyes [90]. Despite the discrepancy at the FC point, the tw-BLA MECI of PSB3 is predicted by all the DFT methods to lie at approximately the same energy level of ca. 3 eV; see Table 1. The BP vectors of PSB3 correspond to the BLA distortion (the x1 vector) and to the torsion mode (the x2 vector) in good agreement with the MRCISD vectors; see Fig. 7. The shape of the S0 and S1 PESs in the vicinity of the MECI point (see Fig. 8) exhibits the same trend as for ethylene; the SSR and SF methods yield a sole crossing point at the MECI and the surfaces are (approximately) symmetrically split when moving along the x2 direction within the branching plane. The adiabatic LRTD-DFT surfaces do not show any interaction between the S0 and S1 states along the x2 direction (see right lower panel of Fig. 8) and there are several crossing points between the states when moving along the x1 direction, which indicates that the


Fig. 7 BP vectors of the twisted-BLA MECI of PSB3 calculated using MRCISD, SSR-BH&HLYP, SF-TD-BH&HLYP, and LR-TD-BH&HLYP methods

Fig. 8 S0 (blue) and S1 (red) PES profiles near the tw-BLA MECI of PSB3. See caption of Fig. 5 for the legend


surface crossing has a linear rather than conical character. Furthermore, the adiabatic LR-TD-DFT method experiences severe convergence problems in the vicinity of the crossing point which prevented us from obtaining a reasonable S0–S1 energy difference plot along a loop around the MECI. Nevertheless, the curves in Fig. 8 corroborate the fact that the adiabatic LR-TD-DFT method is incapable of yielding the correct topography of the S0/S1 crossing point and that the other two DFT methods, SSR and SF, are thus preferred for studying the non-adiabatic dynamics of the excited state of PSB3 and related molecules.

4.3

Ketene

Ketene displays an S0/S1 CI along the C–C–O bending mode which, as first analyzed by Yarkony [97], originates from a crossing between the 1A’ and 1A” electronic states. This crossing is linear under the Cs symmetry constraint; however, it becomes conical when the symmetry restriction is lifted, as there emerges another direction which lifts the degeneracy between these states. Hence, the BP vectors of this CI correspond to the a0 -symmetric C–C–O bending mode (x1) and the out-ofplane H–C–C–O torsional mode (x2). As shown in Fig. 9, the DFT methods yield BP vectors in very good agreement with the MRCISD vectors. However, all the DFT methods predict the MECI point Fig. 9 BP vectors of twisted MECI of ketene calculated using MRCISD, SSR-BH&HLYP, SF-TDBH&HLYP, and LR-TDBH&HLYP methods


Fig. 10 S0 (blue) and S1 (red) PES profiles near the MECI of ketene. See caption of Fig. 5 for the legend

at a somewhat higher, by ca. 0.4–0.8 eV, energy level than MRCISD. This difference may likely be caused by a somewhat stiffer C–C bond predicted by the DFT calculations. This can be asserted from the x1 vectors in Fig. 9, which describes a C–C bond stretching combined with the C–C–O bending motion: there is a notably greater degree of bond stretching in the MRCISD vector. As can be seen from Fig. 10, which shows the S0 and S1 PES profiles along the BP vectors, this motion noticeably destabilizes the S0 (1A0 ) state in which the C–C π-bonding orbital is doubly occupied. It is therefore plausible that MRCISD somewhat underestimates the strength of the C–C bond in ketene because of the rather small basis set and the limited active space employed in the calculations [59]. Although the DFT methods agree with one another on the strength of the C–C bond in ketene, the TD-DFT method differs from the other two methods as to the character of the S0/S1 crossing; whereas SSR and SF predict it to be conical, TD yields a linear crossing. This is confirmed by the plots of the S0–S1 energy difference along a loop around the minimal energy crossing point in Fig. 11. In this plot, a loop with a greater radius, 0.02 Å, was chosen to bypass certain SCF convergence problems in the DFT calculations with a smaller radius (0.002 Å was used in Fig. 6). The adiabatic LR-TD-DFT energy difference in Fig. 11 (dotted curve) again crosses the zero line twice which indicates a linear crossing. The same linear character of the S0/S1 PES crossing can be seen in the rightmost plots in Fig. 10, which show the PES profiles along the BP vectors. Similar to the cases of ethylene and PSB3, there is no coupling between the crossing states in the adiabatic LR-TD-DFT curves. By contrast, the SSR and SF methods correctly predict (nearly) the symmetric splitting of the states which indicates a proper description of the coupling and therefore a correct conical character of the PES crossing point.

M. Huix-Rotllant et al. Fig. 11 The S0–S1 energy difference (in kcal/mol) along a loop around the MECI point of ketene. See text for details of the loop. See caption of Fig. 6 for the legend

5 Conclusions and Outlook In this chapter we have analyzed the ability of practically available computational schemes based on density functional theory to describe conical intersections and the ground- and excited-state potential energy surfaces in their vicinity. We have investigated three DFT-based methods, namely the popular LR-TD-DFT approach as well as the less common SF-TD-DFT and SI-SA-REKS methods, and have evaluated their performance by comparing with reference results from standard MRCISD-WFT calculations. In accord with previous analysis from the literature [40, 41, 43], we find that the LR-TD-DFT method is incapable of properly describing the topography of the S0/S1 crossing point, yielding a linear instead of a conical intersection seam (i.e., M 1 rather than M 2 dimensionality). In LR-TD-DFT, the excited states are obtained as the response states and thus do not couple with the ground state (the reference state), which leads to the observed linear character of the intersection. In addition, there are generally serious SCF convergence problems near the crossing point in the conventional DFT calculation, which often result in aborted TD-DFT calculations. These shortcomings render the LR-TD-DFT method largely useless for modeling non-adiabatic relaxation processes. The other two methods considered here, SF-TD-DFT and SI-SA-REKS, do not suffer from the erratic SCF convergence near the crossing point and yield the correct dimensionality of the CI seam and the correct shape of the S0 and S1 PESs in its vicinity. Being based on ensemble DFT, the SI-SA-REKS method is capable of correctly treating the multi-reference character of the nearly degenerate crossing states. This method has the further advantage of not suffering from the erroneous spin-contamination which plagues the SF-TD-DFT states. Although the use of the high-spin (triplet) state as the reference state in SF-TD-DFT removes the SCF convergence issues near the surface crossing points, the spin contamination inherent in this approach often complicates the identification of proper excited


configurations and results in the occurrence of unphysical mixed-spin states. The use of ad hoc spin-purification procedures [96] should be discouraged as it destroys the correct dimensionality of the CI seam [43]. In view of the widespread availability of the LR-TD-DFT method in quantum chemical codes, we would like to raise concerns about using this method for investigating the dynamics of excited states. It is our hope that the arguments presented in this chapter can help researchers in this area to make reasonable choices when selecting the computational methodology for applications.

References 1. McNaught AD, Wilkinson A (1997) IUPAC. Compendium of chemical terminology, 2nd edn. (The ”Gold Book”). Blackwell, Oxford 2. Mezey PO (1987) Potential energy hypersurfaces. Elsevier, New York 3. Allinger NL (1976) In: Gold V, Bethell D (eds) Advances in physical organic chemistry, vol 13. Academic, London, pp 1–82 4. Born M, Oppenheimer R (1927) Ann Phys 84:457 5. Baer M (2006) Beyond Born–Oppenheimer: electronic nonadiabatic coupling terms and conical intersections. Wiley, Hoboken 6. Atchity GJ, Xantheas SS, Ruedenberg K (1991) J Chem Phys 95:1862 7. Domcke W, Yarkony DR, K€ oppel H (eds) (2004) Conical intersections. Electronic structure, dynamics and spectroscopy. Advanced series in physical chemistry, vol 15. World Scientific, Singapore 8. Domcke W, Yarkony DR, K€ oppel H (eds) (2011) Conical intersections. Theory, computation and experiment. Advanced series in physical chemistry, vol. 17 World Scientific, Singapore 9. Yarkony DR (1996) Rev Mod Phys 68:985 10. Bernardi F, Olivucci M, Robb MA (1996) Chem Soc Rev 25:321 11. Yarkony DR (2004) In: Domcke W, Yarkony DR, K€ oppel H (eds) Conical intersections. electronic structure, dynamics and spectroscopy. Advanced series in physical chemistry, vol 15. World Scientific, Singapore, pp 41–127 12. Migani A, Olivucci M (2004) In: Domcke W, Yarkony DR, K€ oppel H (eds) Conical intersections. electronic structure, dynamics and spectroscopy. Advanced series in physical chemistry, vol 15. World Scientific, Singapore, pp 271–320 13. Butler LJ (1998) Annu Rev Phys Chem 49:125 14. Soto J, Arenas JF, Otero JC, Pelez D (2006) J Phys Chem A 110:8221 15. Hund F (1927) Z Phys 40:742 16. von Neumann J, Wigner E (1929) Physik Z 30:467 17. Teller E (1937) J Phys Chem 41:109 18. Herzberg G, Longuet-Higgins HC (1963) Discuss Faraday Soc 35:77 19. Longuet-Higgins HC (1975) Proc R Soc Lond Ser A 344:147 20. Truhlar DG, Mead CA (2003) Phys Rev A 68:032501 21. Polanyi JC (1972) Acc Chem Res 5:161 22. Polanyi JC (1987) Science 236:680 23. Sellner B, Barbatti M, Lischka H (2009) J Chem Phys 131:024312 24. Robb MA (2011) In: Domcke W, Yarkony DR, K€oppel H (eds) Conical intersections. Theory, computation and experiment. Advanced series in physical chemistry, vol 17. World Scientific, Singapore, pp 3–50

M. Huix-Rotllant et al. 25. Docken KK, Hinze J (1972) J Chem Phys 57:4928 ¨ hrn J (eds) 26. Ruedenberg K (1976) K.R. Sundberg. In: Calais JL, Goscinski O, Linderberg J, O Quantum science. Plenum, New York, pp 505–515 27. Ruedenberg K (1979) In: Report on the NRCC 1978 workshop on post-Hartree–Fock quantum chemistry. Lawrence Berkeley Laboratory, Univ. of California, Report LBL 8233, UC4, CONF 780883, pp 46–64 28. Roos BO (1987) In: Lawley KP (ed) Ab initio methods in quantum chemistry II. Wiley, New York, pp 399–446 29. Andersson K, Malmqvist P, Roos BO (1992) J Chem Phys 96:1218 30. Shavitt I (1977) In: Schaefer HF III (ed) Modern theoretical chemistry, vol 3. Methods of electronic structure theory. Plenum, New York, pp 189–275 31. Zewail AH (2000) J Phys Chem A 104:5660 32. Zewail AH (2010) Chem Phys 378:1 33. Hohenberg P, Kohn W (1964) Phys Rev 136:B864 34. Kohn W, Sham LJ (1965) Phys Rev 140:A1133 35. Gaudoin R, Burke K (2004) Phys Rev Lett 93:173001 36. Gaudoin R, Burke K (2005) Phys Rev Lett 94:029901 37. Ziegler T, Krykunov M, Autschbach J (2014) J Chem Theory Comput 10:3980 38. Casida ME, Jamorski C, Bohr F, Guan JO, Salahub DR (1994) In: Karna SP, Yeates AT (ed) Nonlinear optical materials: theory and modeling, ACS symposium series, vol 628, Div Comp Chem, 1996. Symposium on nonlinear optical materials – theory and modeling, at the 208th national meeting of the American-Chemical-Society, Washington, DC, Aug 21–25, 1994, pp 145–163 39. Marques MAL, Gross EKU (2003) In: Fiolhais C, Nogueira F, Marques MAL (eds) A primer in density-functional theory. Lecture notes in physics, vol 620. Springer, Berlin, pp 144–184 40. Levine BG, Ko C, Quenneville J, Martıńez TJ (2006) Mol Phys 104:1039 41. Yang S, Martıńez TJ (2011) In: Domcke W, Yarkony DR, K€ oppel H (eds) Conical intersections. Theory, computation and experiment, advanced series in physical chemistry, vol 17. World Scientific, Singapore, pp 347–374 42. Huix-Rotllant M, Filatov M, Gozem S, Schapiro I, Olivucci M, Ferre´ N (2013) J Chem Theory Comput 9:3917 43. Gozem S, Melaccio F, Valentini A, Filatov M, Huix-Rotllant M, Ferre´ N, Frutos LM, Angeli C, Krylov AI, Granovsky AA, Lindh R, Olivucci M (2014) J Chem Theory Comput 10:3074 44. Shao Y, Head-Gordon M, Krylov AI (2003) J Chem Phys 118:4807 45. Wang F, Ziegler T (2004) J Chem Phys 121:12191 46. Rinkevicius Z, Vahtras O, Ågren H (2010) J Chem Phys 133:114104 47. Bernard YA, Shao Y, Krylov AI (2012) J Chem Phys 136:204103 48. Huix-Rotllant M, Natarajan B, Ipatov A, Wawire CM, Deutsch T, Casida ME (2010) Phys Chem Chem Phys 12:12811 49. Gross EKU, Oliveira LN, Kohn W (1988) Phys Rev A 37:2805 50. Gross EKU, Oliveira LN, Kohn W (1988) Phys Rev A 37:2809 51. Oliveira LN, Gross EKU, Kohn W (1988) Phys Rev A 37:2821 52. Franck O, Fromager E (2014) Mol Phys 112:1684 53. Lieb EH (1983) Int J Quant Chem 24:243 54. Filatov M, Shaik S (1999) Chem Phys Lett 304:429 55. Filatov M, Shaik S (2000) J Phys Chem A 104:6628 56. Moreira IDPR, Costa R, Filatov M, Illas F (2007) J Chem Theory Comput 3:764 57. Kazaryan A, Heuver J, Filatov M (2008) J Phys Chem A 112:12980 58. Filatov M (2013) J Chem Theory Comput 9:4526 59. Nikiforov A, Gamez JA, Thiel W, Huix-Rotllant M, Filatov M (2014) J Chem Phys 141:124122

Description of Conical Intersections with Density Functional Methods 60. Send R, Sundholm D (2007) J Phys Chem A 111:8766 61. Tapavicza EE, Tavernelli I, R€ othlisberger U, Filippi C, Casida ME (2008) J Chem Phys 129:124108 62. Baranovskii VI, Sizova OV (2008) J Struct Chem 49:803 63. Delchev VB, Ivanova IP (2012) Monatshefte Chem 143:1141 64. Malisˇ M, Loquais Y, Gloaguen E, Biswal HS, Piuzzi F, Tardivel B, Brenner V, Broquier M, Jouvet C, Mons M, Dosˇlic´ N, Ljubic´ I (2012) J Am Chem Soc 134:20340 65. Baranovskii VI, Mal’tsev DA, Sizova OV (2013) Russ Chem Bull 61:973 66. Kumar A, Sevilla MD (2013) Photochem Photobiol Sci 12:1328 67. Zwijnenburg MA (2013) Phys Chem Chem Phys 15:11119 68. Moreno M, Ortiz-Sanchez JM, Gelabert R, Lluch JM (2013) Phys Chem Chem Phys 15:20236 69. Kobayashi R, Amos RD (2013) Mol Phys 111:1574 70. Talhi O, Lopes GR, Santos SM, Pinto DCGA, Silva AMS (2014) J Phys Org Chem 27:756 71. Wiebeler C, Bader CA, Meier C, Schumacher S (2014) Phys Chem Chem Phys 16:14531 72. Belfon KAA, Gough JD (2014) Chem Phys Lett 593:174 73. Baranovskii VI, Maltsev DA (2014) Comp Theor Chem 1043:71 74. Yarkony DR (2001) J Phys Chem A 105:6277 75. Yarkony DR (2000) J Chem Phys 112:2111 76. Toniolo A, Ben-Nun M, Martıńez TJ (2002) J Phys Chem A 106:4679 77. Maeda S, Ohno K, Morokuma K (2010) J Chem Theory Comput 6:1538 78. Bonacˇic´-Koutecky´ V, Koutecky´ J, Michl J (1987) Angew Chem Int Ed 26:170 79. Haas Y, Cogan S, Zilberg S (2005) Int J Quant Chem 102:961 80. Dick B, Haas Y, Zilberg S (2008) Chem Phys 347:65 81. Filatov M, Olivucci M (2014) J Org Chem 79:3587 82. Schipper PRT, Gritsenko OV, Baerends EJ (1998) Theor Chem Acc 99:329 83. Schipper PRT, Gritsenko OV, Baerends EJ (1999) J Chem Phys 111:4056 84. Morrison RC (2002) J Chem Phys 117:10506 85. Salem L, Rowland C (1972) Angew Chem Int Ed 11:92 86. Ziegler T, Rauk A, Baerends EJ (1977) Theor Chim Acta 43:261 87. Filatov M, Shaik S (1998) Chem Phys Lett 288:689 88. Krykunov M, Ziegler T (2013) J Chem Theory Comput 9:2761 89. Krykunov M, Seth M, Ziegler T (2014) J Chem Phys 140:18A502 90. Filatov M, Huix-Rotllant M (2014) J Chem Phys 141:024112 91. Filatov M (2014) J Chem Phys 141:124123 92. Maitra NT, Zhang F, Cave R, Burke K (2004) J Chem Phys 120:5932 93. Cave RJ, Zhang F, Maitra NT, Burke K (2004) Chem Phys Lett 389:39 94. Huix-Rotllant M, Ipatov A, Rubio A, Casida ME (2011) Chem Phys 391:120 95. Cordova F, Doriol LJ, Ipatov A, Casida ME, Filippi C, Vela A (2007) J Chem Phys 127:164111 96. Xu X, Gozem S, Olivucci M, Truhlar DG (2013) J Phys Chem Lett 4:253 97. Yarkony DR (1999) J Phys Chem A 103:6658 98. Gozem S, Schapiro I, Ferre´ N, Olivucci M (2012) Science 337:1225 99. Schreiber M, Silva-Junior MR, Sauer SPA, Thiel W (2008) J Chem Phys 128:134110 100. Leopold DG, Hemley RJ, Vaida V, Roebber JL (1981) J Chem Phys 75:4758 101. Valsson O, Angeli C, Filippi C (2012) Phys Chem Chem Phys 14:11015 102. Chiang SY, Bahou M, Wu YJ, Lee YP (2002) J Chem Phys 117:4306 103. Rice JK, Baronavski AP (1992) J Phys Chem 96:3359 104. Filippi C, Zaccheddu M, Buda F (2009) J Chem Theory Comput 5:2074 105. Ben-Nun M, Martıńez TJ (2000) Chem Phys 259:237 106. Quenneville J, Martıńez TJ (2003) J Phys Chem A 107:829

M. Huix-Rotllant et al. 107. Barbatti M, Paier J, Lischka H (2004) J Chem Phys 121:11614 108. Virshup AM, Chen J, Martıńez TJ (2012) J Chem Phys 137:22A519 109. Gozem S, Huntress M, Schapiro I, Lindh R, Granovsky AA, Angeli C, Olivucci M (2012) J Chem Theory Comput 8:4069 110. Gozem S, Krylov AI, Olivucci M (2013) J Chem Theory Comput 9:284

Switching Exciton Pulses Through Conical Intersections.

Density-functional description of electrides.

Assessment of approximate computational methods for conical intersections and branching plane vectors in organic molecules.

Conical intersections and low-lying electronic states of tetrafluoroethylene.

Communication: conical intersections between vibrationally adiabatic surfaces in methanol.

Derivative Couplings between Time-Dependent Density Functional Theory Excited States in the Random-Phase Approximation Based on Pseudo-Wavefunctions: Behavior around Conical Intersections.

Dissipative dynamics at conical intersections: simulations with the hierarchy equations of motion method.

Density of convex intersections and applications.

Dynamical, spectroscopic and computational imaging of bond breaking in photodissociation: roaming and role of conical intersections.

Density-functional expansion methods: Grand challenges.

Monitoring conical intersections in the ring opening of furan by attosecond stimulated X-ray Raman spectroscopy.

Conical intersections and diabatic potential energy surfaces for the three lowest electronic singlet states of H3 (+).

Ab initio quantum dynamical analysis of ultrafast nonradiative transitions via conical intersections in pyrazine.

Conical Intersections from Particle-Particle Random Phase and Tamm-Dancoff Approximations.

Non-adiabatic dynamics close to conical intersections and the surface hopping perspective.

Designing conical intersections for light-driven single molecule rotary motors: from precessional to axial motion.

Probing non-adiabatic conical intersections using absorption, spontaneous Raman, and femtosecond stimulated Raman spectroscopy.

Non-radiative relaxation of UV photoexcited phenylalanine residues: probing the role of conical intersections by chemical substitution.

Constricted Variational Density Functional Theory Approach to the Description of Excited States.

Geometric phase effects in low-energy dynamics near conical intersections: a study of the multidimensional linear vibronic coupling model.

Description of electron transfer in the ground and excited states of organic donor-acceptor systems by single-reference and multi-reference density functional methods.

Density Functional Partition Theory with Fractional Occupations.

On the accuracy of density functional theory and wave function methods for calculating vertical ionization energies.

Performance of Several Density Functional Theory Methods on Describing Hydrogen-Bond Interactions.