For conical intersections of two states (I,J = I + 1) the vectors defining the branching or g-h plane, the energy difference gradient vector gI,J, and the interstate coupling vector hI,J, can be made orthogonal by a one parameter rotation of the degenerate electronic eigenstates. The representation obtained from this rotation is used to construct the parameters that describe the vicinity of the conical intersection seam, the conical parameters, sI,Jx (R), sI,Jy (R), gI,J(R), and hI,J(R). As a result of the orthogonalization these parameters can be made continuous functions of R, the internuclear coordinates. In this work we generalize this notion to construct continuous parametrizations of conical intersection seams of three or more states. The generalization derives from a recently introduced procedure for using non-degenerate electronic states to construct coupled diabatic states that represent adiabatic states coupled by conical intersections. The procedure is illustrated using the seam of conical intersections of three states in parazolyl as an example.

Recently there has been renewed interest in the role of accidental three state conical intersections in nonadiabatic processes.1–13 Studies dealing with: (i) the location of points of conical intersection;4,5,14 (ii) the seam and branching spaces which are of dimension Nint-5 and 5, respectively;15 (iii) the degeneracy space, a subspace of the branching space in which the degeneracy of two states is preserved;16 (iv) the geometric phase;1,9,17,18 and (v) dynamical studies of actual11,12 and model13 systems, have been reported.

The description of the branching, or g-h plane for intersections of two states of the same symmetry is key to interpreting and determining the dynamics near two state conical intersections. As described in Sec. II, see also Ref. 19, for two state intersections the branching plane is defined by two vectors gI,J, the energy difference gradient vector and hI,J, the interstate coupling vector. Since there are two degenerate states at a point of conical intersection these states are defined up to a one parameter rotation. When this flexibility is used to make the gI,J and hI,J vectors orthogonal,20 the parameters that define the topography of the intersection seam, the conical parameters21 |$s_x^{I,J}$| s x I , J (R), |$s_y^{I,J}$| s y I , J (R), gI,J(R), and hI,J(R) can be made continuous functions of R, the atom centered cartesian coordinates. This continuity is essential for the polynomial representations of diabatic states we use to describe adiabatic states coupled by conical intersection seams.22 The orthogonalization procedure, which defines the arbitrary rotation of the two degenerate states, can be exploited in algorithms for locating conical intersections.23 From a pedagogic perspective the rotated gI,J and hI,J vectors reflect the relevant molecular symmetry an attribute which is usually absent in the original or nascent vectors. In this work we generalize these ideas to conical intersections of three or more states. Section II outlines the theory. Section III demonstrates its application to pyrazolyl, an azolyl for which previous work has demonstrated the existence of low-lying three state intersections.5,24,25 Section IV summarizes and concludes.

Let the cJ(R) satisfy the electronic Schrödinger equation in a configuration state function (CSF) basis,

\begin{equation} [ {{\bf H}^{CSF} ({\bf R}) - {\bf I}E^{a,J,(ab)} ({\bf R})} ]{\bf c}^J ({\bf R}) = {\bf 0}, \end{equation}
[ H C S F ( R ) I E a , J , ( a b ) ( R ) ] c J ( R ) = 0 ,
(1)

where R = R0 + δR. In the vicinity of an N state conical intersection at R0, taken to be the zero of energy, through first order in nuclear coordinates, Eq. (1) becomes6 

\begin{eqnarray} && {\bf H}^d ({\bf R}){\bf d}^I ({\bf R})\nonumber\\ &&\quad \equiv \left( {\begin{array}{c@\quad c@\quad c@\quad c} {{\bf h}^{1,1} \cdot \delta {\bf R}} & & & \\[6pt] {{\bf h}^{2,1} \cdot \delta {\bf R}} & {{\bf h}^{2,2} \cdot \delta {\bf R}} & & \\[6pt] \cdot & \cdot & \cdot & \, \\[6pt] {{\bf h}^{N,1} \cdot \delta {\bf R}} & {{\bf h}^{N,2} \cdot \delta {\bf R}} & \cdot & {{\bf h}^{N,N} \cdot \delta {\bf R}} \end{array}} \right)\left( {\begin{array}{c} {d_1^I ({\bf R})} \\[6pt] {d_2^I ({\bf R})} \\[6pt] \cdot \\[6pt] {d_N^I (({\bf R})} \end{array}} \right)\nonumber\\ &&\quad = E^{a,I,(1)} ({\bf R})\left( {\begin{array}{c} {d_1^I ({\bf R})} \\[6pt] {d_2^I ({\bf R})} \\[6pt] \cdot \\[6pt] {d_N^I ({\bf R})} \end{array}} \right), \end{eqnarray}
H d ( R ) d I ( R ) h 1 , 1 · δ R h 2 , 1 · δ R h 2 , 2 · δ R · · · h N , 1 · δ R h N , 2 · δ R · h N , N · δ R d 1 I ( R ) d 2 I ( R ) · d N I ( ( R ) = E a , I , ( 1 ) ( R ) d 1 I ( R ) d 2 I ( R ) · d N I ( R ) ,
(2a)

where

\begin{equation} {\bf h}^{I,J} ({\bf R}) = {\bf c}^I ({\bf R}^0)^\dag \nabla {\bf H}^{CSF} ({\bf R}){\bf c}^J ({\bf R}^0). \end{equation}
h I , J ( R ) = c I ( R 0 ) H C S F ( R ) c J ( R 0 ) .
(2b)

As described in Ref. 26 (see also Ref. 27) the cI(R0) are used as a geometry independent transformation of the generally applicable HCSF, to Hd an NxN matrix whose spectrum agrees with that of HCSF through first order in δR near R0. For both Hd and HCSF the N degenerate states are determined up to an [N] ≡ N(N − 1)/2 parameter orthogonal geometry dependent transformation, U. This indeterminacy makes the hI,J(R) discontinuous functions of R. We have previously discussed20 the two state conical intersection case where N = 2 and [N] = 1. In this case we introduced the orthogonality requirement,

\begin{equation} {\bf g}^{I,J} ({\bf R}) \cdot {\bf h}^{I,J} ({\bf R}) = 0, \end{equation}
g I , J ( R ) · h I , J ( R ) = 0 ,
(3a)

where

\begin{equation} {\bf g}^{I,J} ({\bf R}) = [{\bf h}^{J,J} ({\bf R}) - {\bf h}^{I,I} ({\bf R})]/2. \end{equation}
g I , J ( R ) = [ h J , J ( R ) h I , I ( R ) ] / 2 .
(3b)

For

\begin{equation} {\bf U}^{(2),I,J} (\theta) = \left( {\begin{array}{c@\quad c} {\cos\theta } & { - \sin\theta } \\[6pt] {\sin\theta } & {\cos\theta } \end{array}} \right), \end{equation}
U ( 2 ) , I , J ( θ ) = cos θ sin θ sin θ cos θ ,
(4a)

Eq. (3a) leads to the requirement

\begin{equation} \tan 4\theta ^{(2),I,J} = \frac{{2{\bf g}^{I,J} \cdot {\bf h}^{I,J} }}{{( {h^{{\rm I,J}} } )^2 - ( {g^{{\rm I,J}} } )^2 }}, \end{equation}
tan 4 θ ( 2 ) , I , J = 2 g I , J · h I , J ( h I , J ) 2 ( g I , J ) 2 ,
(4b)

where gI,J = gI,J∥ and hI,J = hI,J∥. The requirement (3a) is a natural choice since in that case Hd has the appealing form

\begin{eqnarray} {\bf H}^d &=& {\bf s}^{I,J} \cdot \delta {\bf R}\left( {\begin{array}{c@\quad c} 1 & 0 \\[6pt] 0 & 1 \end{array}} \right) - g^{I,J} \delta R^g \left( {\begin{array}{c@\quad c} 1 & 0 \\[6pt] 0 & { - 1} \end{array}} \right)\nonumber\\ && +\, h^{I,J} \delta R^h \left( {\begin{array}{c@\quad c} 0 & 1 \\[6pt] 1 & 0 \end{array}} \right), \end{eqnarray}
H d = s I , J · δ R 1 0 0 1 g I , J δ R g 1 0 0 1 + h I , J δ R h 0 1 1 0 ,
(5)

where sI,J = (hI,I + hJ,J)/2, δRg = gI, J/gI, J · δR, and δRh = hI, J/hI, J · δR.

The generalization to three or more states is not obvious since there are not enough parameters in U(N) (only [N]) to make [N + 1] − 1 mutually orthogonal linear combinations of the hI,J. Here [N + 1] − 1 appears since the trace of Hd is excluded as it does not contribute to the state splitting. In this work we show that for an N state degeneracy of the states J, J(min)J ≤ J(max) = J(min) + (N − 1), the [N] conditions,

\begin{equation} {\bf g}^{I,J} ({\bf R}) \cdot {\bf h}^{I,J} ({\bf R}) = 0\quad J^{(\min)} \le I < J \le J^{(\max)}, \end{equation}
g I , J ( R ) · h I , J ( R ) = 0 J ( min ) I < J J ( max ) ,
(6)

determine U(N) and lead to continuous hI,J(R). A condition similar to Eq. (6) was introduced and carefully analyzed in our recent work describing the representation of adiabatic potential energy surfaces coupled by conical intersections.28 Here we provide a concise recapitulation of that analysis (Eqs. (7a)–(7c), (8), (9a)–(9c), (10a), (10b), (11a)–(11d), (12a)–(12c), and (13)–(16)) sufficient to explain the present procedure.

The desired orthogonal transformation U(N) is constructed as a product of 2 × 2 rotations UN,i,ji, j) (see Sec. III C),

\begin{equation} {\bf U}^{(N)} = \prod\limits_m {{\bf U}^{N,i(m),j(m)} (\theta ^{i(m),j(m)} (m))}, \end{equation}
U ( N ) = m U N , i ( m ) , j ( m ) ( θ i ( m ) , j ( m ) ( m ) ) ,
(7a)

where UN, i, ji, j) is the Givens rotation matrix between basis i and j:

\begin{equation} U_{k,l}^{N,i,j} (\theta ^{i,j}) = \delta _{k,l}, \end{equation}
U k , l N , i , j ( θ i , j ) = δ k , l ,
(7b)

except for

\begin{eqnarray} U_{i,i}^{N,i,j} (\theta ^{i,j}) &=& \cos\theta ^{i,j} = U_{j,j}^{N,i,j} (\theta ^{i,j})\;{\rm and}\; -U_{i,j}^{N,i,j} (\theta ^{i,j})\nonumber\\[4pt] &=& \sin\theta ^{i,j} = U_{j,i}^{N,i,j} (\theta ^{i,j}). \end{eqnarray}
U i , i N , i , j ( θ i , j ) = cos θ i , j = U j , j N , i , j ( θ i , j ) and U i , j N , i , j ( θ i , j ) = sin θ i , j = U j , i N , i , j ( θ i , j ) .
(7c)

Below we will suppress the i,j superscript on θi, j when no confusion will result. Instead of dealing directly with Eq. (6) we showed that those equations are equivalent to minimization of the function,

\begin{equation} S = \sum\limits_{K < L} {{\bf h}^{K,L} } \cdot {\bf h}^{K,L}, \end{equation}
S = K < L h K , L · h K , L ,
(8)

which is comprised of [N] (K,L) pairs.

To establish that equivalence, we proceed as follows. Denote the wave function of electronic state I in any basis, such as the CSF basis, as dI, and the corresponding wave function after the Givens rotation as dI(θ). Then the effect of one individual Givens rotation on the wave functions can be written in vector form as

\begin{equation} \left( {\begin{array}{c} {{\bf d}^1 (\theta)} \\[8pt] {{\bf d}^2 (\theta)} \\[8pt] \vdots \\[8pt] {{\bf d}^N (\theta)} \end{array}} \right) = {\bf U}^{N,I,J} (\theta)\left( {\begin{array}{c} {{\bf d}^1 } \\[8pt] {{\bf d}^2 } \\[8pt] \vdots \\[8pt] {{\bf d}^N } \end{array}} \right), \end{equation}
d 1 ( θ ) d 2 ( θ ) d N ( θ ) = U N , I , J ( θ ) d 1 d 2 d N ,
(9a)

so that the rotated vectors satisfy

\begin{equation} \left\{ {\begin{array}{c} {{\bf d}^I (\theta) = {\bf d}^I \cos\theta - {\bf d}^J \sin\theta } \\[8pt] {{\bf d}^J (\theta) = {\bf d}^I \sin\theta + {\bf d}^J \cos\theta } \\[8pt] {{\bf d}^K (\theta) = {\bf d}^K,{\rm where }K \ne I,J} \end{array}} \right.. \end{equation}
d I ( θ ) = d I cos θ d J sin θ d J ( θ ) = d I sin θ + d J cos θ d K ( θ ) = d K , where K I , J .
(9b)

Therefore,

\begin{equation} \left\{ {\begin{array}{c} {{\bf h}^{I,J} (\theta) = - {\bf g}^{I,J} \sin2\theta + {\bf h}^{I,J} \cos2\theta } \\[8pt] {{\bf g}^{I,J} (\theta) = {\bf h}^{I,J} \sin2\theta + {\bf g}^{I,J} \cos2\theta } \\[8pt] {{\bf h}^{I,K} (\theta) = {\bf h}^{I,K} \cos\theta - {\bf h}^{J,K} \sin\theta } \\[8pt] {{\bf h}^{J,K} (\theta) = {\bf h}^{J,K} \cos\theta + {\bf h}^{I,K} \sin\theta } \end{array}} \right., \end{equation}
h I , J ( θ ) = g I , J sin 2 θ + h I , J cos 2 θ g I , J ( θ ) = h I , J sin 2 θ + g I , J cos 2 θ h I , K ( θ ) = h I , K cos θ h J , K sin θ h J , K ( θ ) = h J , K cos θ + h I , K sin θ ,
(9c)

where KI, J. Other h vectors remain unchanged.

Differentiating Eq. (9c) we showed that28 

\begin{equation} \frac{{\partial S}}{{\partial \theta ^{I,J} }} = 4{\bf g}^{I,J} \cdot {\bf h}^{I,J}. \end{equation}
S θ I , J = 4 g I , J · h I , J .
(10a)

For the current dI to satisfy Eq. (6)gI,J · hI,J must vanish for all pairs. From Eq. (10a), Eq. (6) is equivalent to the condition that S is stationary with respect to any rotation, that is for any pairs of states I < J,

\begin{equation} \frac{{\partial S}}{{\partial \theta ^{I,J} }} = 0. \end{equation}
S θ I , J = 0 .
(10b)

Since S, from Eq. (8), is a sum of norms, it is non-negative. Being a non-negative continuous function, S is guaranteed to have a minimum. The solution to Eq. (6) is therefore also guaranteed.

This solution can be obtained by a sequence of 2 × 2 rotations.28 In order to show this, S(θ) is first divided into 3 parts:

\begin{eqnarray} S( \theta ) &=& \| {{\bf h}^{I,J} (\theta)} \|^2 + \sum\limits_{K \ne I,J} {( {\| {{\bf h}^{I,K} (\theta)} \|^2 + \| {{\bf h}^{J,K} (\theta)} \|^2 } )}\nonumber\\ && + \sum\limits_{K < L,K \ne I,L \ne J} {\| {{\bf h}^{K,L} (\theta)} \|^2 }, \end{eqnarray}
S ( θ ) = h I , J ( θ ) 2 + K I , J ( h I , K ( θ ) 2 + h J , K ( θ ) 2 ) + K < L , K I , L J h K , L ( θ ) 2 ,
(11a)

where using Eq. (9c) the three terms from the left to the right in Eq. (11a) are, respectively,

\begin{eqnarray} \| {{\bf h}^{I,J} (\theta)} \|^2 &=& \frac{1}{2}( {( {h^{I,J} } )^2 + ( {g^{I,J} } )^2 } ) + \frac{1}{2} (( {h^{I,J} } )^2\nonumber\\ && -\, ( {g^{I,J} } )^2 )\cos 4\theta - {\bf h}^{I,J} \cdot {\bf g}^{I,J} \sin4\theta,\qquad \end{eqnarray}
h I , J ( θ ) 2 = 1 2 ( ( h I , J ) 2 + ( g I , J ) 2 ) + 1 2 ( ( h I , J ) 2 ( g I , J ) 2 ) cos 4 θ h I , J · g I , J sin 4 θ ,
(11b)
\begin{eqnarray} &&\hspace*{-5pt}\| {{\bf h}^{I,K} (\theta)} \|^2 + \| {{\bf h}^{J,K} (\theta)} \|^2\nonumber\\ &&\hspace*{-5pt}\quad = \| {{\bf h}^{I,K} \cos\theta - {\bf h}^{J,K} \sin\theta } \|^2 + \| {{\bf h}^{J,K} \cos\theta + {\bf h}^{I,K} \sin\theta } \|^2 \nonumber\\ &&\hspace*{-5pt}\quad= \| {{\bf h}^{I,K} } \|^2 \cos^2 \theta \!-\! 2{\bf h}^{I,K} \cdot {\bf h}^{J,K} \cos\theta \sin\theta \!+\! \| {{\bf h}^{J,K} } \|^2 \sin^2 \theta \nonumber\\ &&\hspace*{-5pt}\qquad+\, \| {{\bf h}^{I,K} } \|^2 \sin^2 \theta + 2{\bf h}^{I,K} \cdot {\bf h}^{J,K} \cos\theta \sin\theta\nonumber\\ &&\hspace*{-5pt}\qquad +\, \| {{\bf h}^{J,K} } \|^2 \cos^2 \theta \nonumber\\ &&\hspace*{-5pt}\quad= \| {{\bf h}^{I,K} } \|^2 + \| {{\bf h}^{J,K} } \|^2, \end{eqnarray}
h I , K ( θ ) 2 + h J , K ( θ ) 2 = h I , K cos θ h J , K sin θ 2 + h J , K cos θ + h I , K sin θ 2 = h I , K 2 cos 2 θ 2 h I , K · h J , K cos θ sin θ + h J , K 2 sin 2 θ + h I , K 2 sin 2 θ + 2 h I , K · h J , K cos θ sin θ + h J , K 2 cos 2 θ = h I , K 2 + h J , K 2 ,
(11c)

and

\begin{equation} \| {{\bf h}^{K,L} (\theta)} \|^2 = \| {{\bf h}^{K,L} } \|^2. \end{equation}
h K , L ( θ ) 2 = h K , L 2 .
(11d)

Since the last two parts are both constant with respect to θ, we only need to find the minimum of the first term ‖hI, J(θ)‖2 with respect to angle θ.

Defining C and α as

\begin{equation} C = \sqrt {\frac{1}{4} [ {( {h^{I,J} } )^2 - ( {g^{I,J} } )^2 } ]^2 + [ {{\bf h}^{I,J} \cdot {\bf g}^{I,J} } ]^2 }, \end{equation}
C = 1 4 [ ( h I , J ) 2 ( g I , J ) 2 ] 2 + [ h I , J · g I , J ] 2 ,
(12a)
\begin{equation} \frac{1}{2}[ {( {h^{I,J} } )^2 - ( {g^{I,J} } )^2 } ] = C\cos\alpha, \end{equation}
1 2 [ ( h I , J ) 2 ( g I , J ) 2 ] = C cos α ,
(12b)
\begin{equation} {\bf h}^{I,J} \cdot {\bf g}^{I,J} = C\sin\alpha, \end{equation}
h I , J · g I , J = C sin α ,
(12c)

Eq. (11b) becomes

\begin{equation} \| {{\bf h}^{I,J} (\theta)} \|^2 = \frac{1}{2} [ {( {h^{I,J} } )^2 + ( {g^{I,J} } )^2 } ] + C\cos( {\alpha + 4\theta } ). \end{equation}
h I , J ( θ ) 2 = 1 2 [ ( h I , J ) 2 + ( g I , J ) 2 ] + C cos ( α + 4 θ ) .
(13)

Note that using Eqs. (12a)–(12c), Eq. (4b) is equivalent to

\begin{equation} C\sin(4\theta + \alpha) = 0, \end{equation}
C sin ( 4 θ + α ) = 0 ,
(14)

which has solutions

\begin{equation} 4\theta + \alpha = (2k - 1)\pi \,{\rm for}\,k{\rm = 1,2}, \end{equation}
4 θ + α = ( 2 k 1 ) π for k = 1 , 2 ,
(15)

where we have chosen the values of θ for which cos(α + 4θ) = −1.

On the other hand, using Eq. (13), this choice of θ yields the relation

\begin{eqnarray} && \| {{\bf h}^{I,J} (\theta)} \|^2\nonumber\\ &&\quad= \frac{1}{2}[ {( {h^{I,J} } )^2 + ( {g^{I,J} } )^2 } ] - C = \frac{1}{2}[ {( {h^{I,J} } )^2 + ( {g^{I,J} } )^2 } ]\nonumber\\ &&\qquad - \frac{1}{2}\sqrt {[ {( {h^{I,J} } )^2 - ( {g^{I,J} } )^2 } ]^2 + 2( {{\bf h}^{I,J} \cdot {\bf g}^{I,J} } )^2 } \nonumber\\ &&\quad \le\, \frac{1}{2}[ {( {h^{I,J} } )^2 + ( {g^{I,J} } )^2 } ] - \frac{1}{2}| {( {h^{I,J} } )^2 - ( {g^{I,J} } )^2 } | \nonumber\\ &&\quad= \min[ {( {h^{I,J} } )^2,( {g^{I,J} } )^2 } ], \end{eqnarray}
h I , J ( θ ) 2 = 1 2 [ ( h I , J ) 2 + ( g I , J ) 2 ] C = 1 2 [ ( h I , J ) 2 + ( g I , J ) 2 ] 1 2 [ ( h I , J ) 2 ( g I , J ) 2 ] 2 + 2 ( h I , J · g I , J ) 2 1 2 [ ( h I , J ) 2 + ( g I , J ) 2 ] 1 2 | ( h I , J ) 2 ( g I , J ) 2 | = min [ ( h I , J ) 2 , ( g I , J ) 2 ] ,
(16)

where min() denotes the smaller of the two arguments. The ≤ in Eq. (16) becomes an equality only when (hI, J · gI, J)2 = 0. Thus (hI, J)2 and hence S are strictly decreasing unless hI, J · gI, J = 0. By choosing the pairs of states I and J with the largest rotation angle every iteration, the rotation will always lower the value of S until hI, J · gI, J = 0 holds for all pairs of states I and J, which is the desired result. Using Eq. (9a) and the definition of hJ,K, dI is updated each iteration, until the solution to Eq. (10b) is obtained. As discussed below the resulting hi,j are used to define the rotated branching space.

Equation (6) may have more than one solution corresponding to local solutions of Eq. (10b). Further changing the phase of a single eigenstate, a reflection not included in Eqs. (7a)–(7c) which changes hI,J but not gI,J, can change the final value of S. This perhaps counterintuitive result was verified by explicit numerical calculations. These two issues are attributable to the range of starting values of the gI,J and hI,J. Methods to deal with these situations are described in the Sec. III.

At R0 a point of N state conical intersection Hd in Eq. (2a) is, to first order in δR, a rigorous diabatic representation. See Eq. (2b). However, it is inappropriate for defining the branching space of this conical intersection since it does not explicitly account for the trace. For points of conical intersection of 3 states, to be considered in Sec. III, the g-h or branching space is defined by the five vectors:8 

\begin{equation} {\bf h}^{[1]} \equiv {\bf h}^{1,2},\quad {\bf h}^{[2]} \equiv {\bf h}^{2,3},\quad {\bf h}^{[3]} \equiv {\bf h}^{1,3}, \end{equation}
h [ 1 ] h 1 , 2 , h [ 2 ] h 2 , 3 , h [ 3 ] h 1 , 3 ,
(17a)
\begin{eqnarray} {\bf h}^{[4]} & \equiv & {\bf g}^{1,2} = - ({\bf h}^{1,1} - {\bf h}^{2,2})/2,\quad \nonumber\\[-7.5pt]\\[-7.5pt] {\bf h}^{[5]} & \equiv & {\bf g}^{1,2,3}= (2{\bf h}^{3,3} - {\bf h}^{1,1} - {\bf h}^{2,2})/3.\nonumber \end{eqnarray}
h [ 4 ] g 1 , 2 = ( h 1 , 1 h 2 , 2 ) / 2 , h [ 5 ] g 1 , 2 , 3 = ( 2 h 3 , 3 h 1 , 1 h 2 , 2 ) / 3 .
(17b)

Also required is the trace vector

\begin{equation} {\bf h}^{[6]} \equiv {\bf s}^{1,2,3} = ({\bf h}^{1,1} + {\bf h}^{2,2} + {\bf h}^{3,3})/3, \end{equation}
h [ 6 ] s 1 , 2 , 3 = ( h 1 , 1 + h 2 , 2 + h 3 , 3 ) / 3 ,
(17c)

which does not contribute to the state splitting and therefore is not a member of the branching space. In this representation, Hd is given by8 

\begin{equation} {\bf H}^{\rm d} {\rm = }\sum\limits_{j = 1}^6 {( {{\bf h}^{[j]} \cdot \delta {\bf R}} ){\bf B}^{[j]} }, \end{equation}
H d = j = 1 6 ( h [ j ] · δ R ) B [ j ] ,
(18)

where

\begin{eqnarray} B^{[1]} &=& \left( {\begin{array}{c@\quad c@\quad c} 0 & 1 & 0 \\[6pt] 1 & 0 & 0 \\[6pt] 0 & 0 & 0 \end{array}} \right),\quad B^{[2]} = \left( {\begin{array}{c@\quad c@\quad c} 0 & 0 & 0 \\[6pt] 0 & 0 & 1 \\[6pt] 0 & 1 & 0 \end{array}} \right),\nonumber\\[-7.5pt]\\[-7.5pt] B^{[3]} &=& \left( {\begin{array}{c@\quad c@\quad c} 0 & 0 & 1 \\[6pt] 0 & 0 & 0 \\[6pt] 1 & 0 & 0 \end{array}} \right),\nonumber \end{eqnarray}
B [ 1 ] = 0 1 0 1 0 0 0 0 0 , B [ 2 ] = 0 0 0 0 0 1 0 1 0 , B [ 3 ] = 0 0 1 0 0 0 1 0 0 ,
(19a)
\begin{eqnarray} B^{[4]} &=& \left( {\begin{array}{c@\quad c@\quad c} { - 1} & 0 & 0 \\[6pt] 0 & 1 & 0 \\[6pt] 0 & 0 & 0 \end{array}} \right),\quad B^{[5]} = \frac{1}{2}\left( {\begin{array}{c@\quad c@\quad c} { - 1} & 0 & 0 \\[6pt] 0 & { - 1} & 0 \\[6pt] 0 & 0 & 2 \end{array}} \right),\nonumber\\[-7.5pt]\\[-7.5pt] B^{[6]} &=& \left( {\begin{array}{c@\quad c@\quad c} 1 & 0 & 0 \\[6pt] 0 & 1 & 0 \\[6pt] 0 & 0 & 1 \end{array}} \right).\nonumber \end{eqnarray}
B [ 4 ] = 1 0 0 0 1 0 0 0 0 , B [ 5 ] = 1 2 1 0 0 0 1 0 0 0 2 , B [ 6 ] = 1 0 0 0 1 0 0 0 1 .
(19b)

Note that Hd does not have the simple form of Eq. (5) since here only h[1] and h[4] are orthogonal. The continuous representation of the branching or g-h space vectors using Eq. (10b) for a practical example is discussed below. For the N state generalization of Eq. (18) see Ref. 8 and also Ref. 29. These works use Lie group techniques to process the B matrices.

The algorithm described in Sec. II B is used here to develop a continuous representation of the branching space for a portion of the three state conical intersection seam in pyrazolyl (CH)3N2, Ref. 5, pictured in Fig. 1. This seam has been implicated in the photoelectron spectrum of the pyrazolide anion, which has been simulated previously.24,25 However to date only isolated points on the seam of three state conical intersections have been discussed.

FIG. 1.

Pyrazolyl at its minimum energy 3-state intersection structure displaying the atom numbering used in this work. Bold numbers indicate the bond lengths (in Å) and angles (in °) of the minimum energy 3-state intersection geometry, while the plain face numbers mark that of the global minimum on the ground state.

FIG. 1.

Pyrazolyl at its minimum energy 3-state intersection structure displaying the atom numbering used in this work. Bold numbers indicate the bond lengths (in Å) and angles (in °) of the minimum energy 3-state intersection geometry, while the plain face numbers mark that of the global minimum on the ground state.

Close modal

The low-lying 1,2,32A states of pyrazolyl exhibit a seam of three-state conical intersections whose description is inherently multi-reference in character and so multi-reference configuration interaction singles and doubles (MR-CISD) wave functions will be employed to describe this seam. The active space for these demonstration calculations is comprised of three orbitals including the two π type orbitals with one additional nodal plane perpendicular to molecular plane, as well as the anti-symmetric linear combination of the two lone pair orbitals one on each of the nitrogen atoms, yielding a 5 electron/3 orbital complete active space (CAS) reference wavefunction, with 3 CSFs that correspond to the 3 electronic states. The molecular orbitals were obtained from a state-averaged multiconfigurational self consistent field (MCSCF) procedure that included all three electronic states with equal weights and employed polarized double-zeta30 atomic basis sets. This basis set, in conjunction with the second-order CI treatment of electron correlation and employing a generalized interacting space restriction approximation, yielded MRCI wavefunctions comprised of 1 127 824 CSFs. All calculations reported in this work were carried out using the COLUMBUS suite of electronic structure codes.31–33 All energies are reported relative to the ground state minimum of pyrazolyl for which E1 = −224.7651475041 a.u.

The minimum energy three state conical intersection structure for pyrazolyl, Q0 with energy E1(Q0) = 4905.19 cm−1, is pictured in Fig. 1. Starting from this point a number of additional points on the seam of three state intersections were determined along each of three paths: (1) the stretching of the C2-H7 bond, (2) the out-of-plane angle of H6 with respect to the molecular plane, and (3) the out-of-plane angle of H7. We desire conical intersections systematically displaced along the path parameter but otherwise randomly placed. Such paths were determined by incrementing the parameter describing the path and then running the conical intersection search without energy optimization. The produced three state conical intersections are degenerate to <0.5 cm‑1. The degenerate energies along these three paths are plotted in Figs. 2(a)–2(c), respectively, to demonstrate their continuity. Nascent g, h, and s vectors and the corresponding vectors after the solution of Eq. (10b), indicated by superscript ∼, are juxtaposed in Figs. 3 and 4, respectively. Figs. 3(a)–3(f) (Figs. 4(a)–4(f)) report the nascent vectors h[j], j = 1–6 (the rotated vectors |${\bf \tilde h}^{[j]}$| h ̃ [ j ] j = 1–6) at the first point along path (3). This point is chosen because while it is nearly identical to the minimum energy 3 state crossing point, Q0, which has C2v symmetry, its geometry carries only C1 point group symmetry, making it a more suitable example of a general point without symmetry. The differences between the corresponding vectors in these figures are seen to be dramatic, except for h[6] which is invariant to any rotation. The orthogonality of |${\bf \tilde h}$| h ̃ 1,2 and |${\bf \tilde g}$| g ̃ 1,2, which is expected, is clearly evident. Also evident is the near C2v symmetry of the rotated vectors, |${\bf \tilde h}^{[j]}$| h ̃ [ j ] j = 1–6, in Fig. 4 which again is expected. No evidence of this symmetry is present in Fig. 3 except of course for h[6].

FIG. 2.

Degenerate energies along displacement paths. Plates (a), (b), (c) report paths 1, 2, 3, respectively.

FIG. 2.

Degenerate energies along displacement paths. Plates (a), (b), (c) report paths 1, 2, 3, respectively.

Close modal
FIG. 4.

(a)–(f) Branching space vectors at the first point along path 3 after orthogonalization.

FIG. 4.

(a)–(f) Branching space vectors at the first point along path 3 after orthogonalization.

Close modal
FIG. 3.

(a)–(f) Branching space vectors at the first point along path 3 before orthogonalization.

FIG. 3.

(a)–(f) Branching space vectors at the first point along path 3 before orthogonalization.

Close modal

Figs. 5–7 consider the continuity of the vectors describing the branching space along the seam paths (1), (2), (3), respectively. At each seam point Qi these figures report (as individual markers) ||h[j](Qi) − h[j](Q0)|| the norm of the difference between the nascent vectors h[j],j = 1–6 and the corresponding reference vectors, h[j](Q0), as well as (solid lines) |$||{\bf \tilde h}^{[j]} ({\bf Q}_i) - {\bf \tilde h}^{[j]} ({\bf Q}_0)||$| h ̃ [ j ] ( Q i ) h ̃ [ j ] ( Q 0 ) the difference norms for the corresponding rotated quantities. The continuity (lack of continuity) is clearly evident in the rotated (nascent) quantities. For each of the solid lines the difference norm is small and a continuous function of the displacement. For the individual markers the difference from the origin vector is a discontinuous function of the displacement, except for h[6]. h[6] = |${\bf \tilde h}^{[6]}$| h ̃ [ 6 ] and are continuous as is expected of the trace.

FIG. 5.

Along path 1: Norm of differences of branching space vectors from the rotated reference vectors in nascent and orthogonalized representations.

FIG. 5.

Along path 1: Norm of differences of branching space vectors from the rotated reference vectors in nascent and orthogonalized representations.

Close modal
FIG. 7.

Along path 3: Norm of differences of branching space vectors from the rotated reference vectors in nascent and orthogonalized representations.

FIG. 7.

Along path 3: Norm of differences of branching space vectors from the rotated reference vectors in nascent and orthogonalized representations.

Close modal
FIG. 6.

Along path 2: Norm of differences of branching space vectors from the rotated reference vectors in nascent and orthogonalized representations.

FIG. 6.

Along path 2: Norm of differences of branching space vectors from the rotated reference vectors in nascent and orthogonalized representations.

Close modal

A comment on the irregular behavior of the nascent states and the corresponding g and h vectors is in order. The “discontinuities” pictured in Figs. 5–7 are a necessary consequence of the search algorithm and the conical topography. Since exact degeneracies are never achieved with a numerical procedure the computed, approximately degenerate, well defined eigenstates, the nascent representation, fall on a continuum of locations near the apex of the three state cone. Because of the conical topography these points give rise to distinct values of gI,J. From Eq. (9c) we see that the distinct values of gI,J correspond to distinct mixings of the N states at the conical intersection. Since at neighboring points on the seam completely different positions on the cones are possible, indeed likely, the nascent states representations behave irregularly along the seam. Our orthogonalization returns all mixtures to a common representation.

The observed continuity is not achieved by simply solving Eq. (10b) as described in Sec. II. This is a consequence of the multiple local solutions to that equation. Instead, a preliminary transformation is first carried out to minimize the quantity

\begin{equation} \Delta = \sum\limits_{I,J} {\| {{\bf h}^{U,I,J} ({\bf Q}_i) - {\bf \tilde h}^{I,J} ({\bf Q}_{ref})} \|^2 }, \end{equation}
Δ = I , J h U , I , J ( Q i ) h ̃ I , J ( Q r e f ) 2 ,
(20)

where |${\bf \tilde h}^{I,J} ({\bf Q}_{ref})$| h ̃ I , J ( Q r e f ) are the rotated vectors at a reference geometry where an orthogonalized representation is already achieved through the aforementioned rotation procedure and hU, I, J is hI,J constructed from U transformed eigenvectors as in Eq. (9b). Every iteration in the determination of U, the phase of the U transformed wave functions may be adjusted to accommodate the need to alter the wave function phase as noted at the end of Sec. II B. Δ has a unique minimum which is determined in a very straightforward manner with numerical gradients and line search techniques. After this preliminary transformation convergence of Eqs. (10a) and (10b) is readily achieved. This preconditioning by ensuring that the same minimum of S is obtained results in a representation that is continuous along the whole seam. In this work only a single reference geometry was used. More sophisticated selection schemes are possible. In future work it will be interesting to study the topography, particularly the local minima, of S.

The nascent and orthogonalized representations of the g-h space of a point of conical intersection of three states contain the same information about the local topography. However the orthogonalized form has several uses for which the nascent form is not suitable. Here we have focused on the geometry dependence of the representation. Only the orthogonalized representation offers the possibility of extrapolating or interpolating that information to infer changes in the local topography associated with motion in the seam space. In a qualitative sense the orthogonalized representation accomplishes this by aligning the coordinate axes at adjacent geometries evincing only the true geometric variations in the branching space vectors. The rotated vectors were also shown to evince the relevant point group symmetry.

The orthogonalized representation finds application in our construction of quasi diabatic representations of coupled adiabatic potential energy surfaces.28 In this case since we fit both energy gradients and derivative couplings we need a way to compare the conical parameters of the fit and ab initio representations of a single point of three state conical intersection. The orthogonalized representation enables that comparison. For more details on the use of conical intersections in coupled surface fitting, see Ref. 28. Further, as in the case of two state conical intersections, the orthogonalized g-h space vectors, by avoiding spurious changes due to the incomplete definition which afflicts the nascent representation, can be used to extrapolate the iterative procedure for locating conical intersections. See Ref. 23 for details in the two state case. This issue will be considered in detail in a future study.

At a point of conical intersection of N states, the electronic states are defined up to the N(N − 1)/2 parameters defining an orthogonal mixing of the degenerate states. This flexibility renders the branching space vectors discontinuous functions of nuclear coordinates along the seam. In this work the orthogonalization approach used remove this flexibility in the two state case is generalized to N state conical intersections. Using the three state seam of conical intersection in pyrazolyl as an example, the orthogonalization condition, Eq. (6), or its equivalent Eq. (10b), was shown to efficiently eliminate the discontinuous nuclear coordinate dependence of the branching space vectors in the seam space, that occurs for the nascent eigenstates. The possible issue of local solutions to the defining equations, which also exists in the two state case, where the g and h vectors can switch, was noted. A method to avoid the issue was demonstrated. Other uses of the rotated representation, including improving the algorithm for locating points of conical intersection, were briefly noted.

This work was supported by National Science Foundation Grant No. CHE-1361121 to D.R.Y.

1.
K. A.
Kistler
and
S.
Matsika
,
J. Chem. Phys.
128
,
215102
(
2008
).
2.
S.
Matsika
, in
Conical Intersections, Theory, Computation and Experiment
, edited by
W.
Domcke
,
D.
Yarkony
, and
H.
Köppel
(
World Scientific
,
Singapore
,
2011
).
3.
S.
Matsika
,
J. Phys. Chem. A.
109
,
7538
(
2005
).
4.
S.
Matsika
and
D. R.
Yarkony
,
J. Chem. Phys.
117
,
6907
(
2002
).
5.
S.
Matsika
and
D. R.
Yarkony
,
J. Am. Chem. Soc.
125
,
12428
(
2003
).
6.
M. S.
Schuurman
and
D. R.
Yarkony
,
J. Chem. Phys.
124
,
124109
(
2006
).
7.
M. S.
Schuurman
and
D. R.
Yarkony
,
J. Chem. Phys.
124
,
244103
(
2006
).
8.
M. S.
Schuurman
and
D. R.
Yarkony
,
J. Phys. Chem. B.
110
,
19031
(
2006
).
9.
S.
Han
and
D. R.
Yarkony
,
J. Chem. Phys.
119
, 1
1561
(
2003
).
10.
J. D.
Coe
and
T. J.
Martinez
,
J. Am. Chem. Soc.
127
,
4560
(
2005
).
11.
J. D.
Coe
and
T. J.
Martínez
,
J. Phys. Chem. A
110
,
618
(
2006
).
12.
J.
Gonzalez-Vazquez
and
L.
Gonzalez
,
ChemPhysChem
11
,
3617
(
2010
).
13.
P.
Krause
and
S.
Matsika
,
J. Chem. Phys.
136
,
034110
(
2012
).
14.
S.
Matsika
and
D. R.
Yarkony
,
J. Am. Chem. Soc.
125
,
10672
(
2003
).
15.
H. C.
Longuet-Higgins
,
Proc. R. Soc. London A
344
,
147
(
1975
).
16.
S. P.
Keating
and
C. A.
Mead
,
J. Chem. Phys.
82
,
5102
(
1985
).
17.
D. E.
Manolopoulos
and
M. S.
Child
,
Phys. Rev. Lett.
82
,
2223
(
1999
).
18.
J.
Samuel
and
A.
Dhar
,
Phys. Rev. Lett.
87
,
260401
(
2001
).
19.
G. J.
Atchity
,
S. S.
Xantheas
, and
K.
Ruedenberg
,
J. Chem. Phys.
95
,
1862
(
1991
).
20.
D. R.
Yarkony
,
J. Chem. Phys.
112
,
2111
(
2000
).
21.
D. R.
Yarkony
,
J. Chem. Phys.
114
,
2601
(
2001
).
22.
X.
Zhu
and
D. R.
Yarkony
,
Mol. Phys.
108
,
2611
(
2010
).
23.
D. R.
Yarkony
,
Faraday Discuss.
127
,
325
(
2004
).
24.
T.
Ichino
,
A. J.
Gianola
,
W. C.
Lineberger
, and
J. F.
Stanton
,
J. Chem. Phys.
125
,
084312
(
2006
).
25.
M. S.
Schuurman
and
D. R.
Yarkony
,
J. Chem. Phys
129
,
064304
(
2008
).
26.
D. R.
Yarkony
,
J. Phys. Chem. A
101
,
4263
(
1997
).
27.
C. A.
Mead
,
J. Chem. Phys.
78
,
807
(
1983
).
28.
X.
Zhu
and
D. R.
Yarkony
,
J. Chem. Phys.
140
,
024112
(
2014
).
29.
B.
Lasorne
,
Adv. Math. Phys.
2014
,
795730
(
2014
).
30.
T. H.
Dunning
 Jr.
,
J. Chem. Phys.
90
,
1007
(
1989
).
31.
H.
Lischka
,
M.
Dallos
,
P.
Szalay
,
D. R.
Yarkony
, and
R.
Shepard
,
J. Chem. Phys.
120
,
7322
(
2004
).
32.
M.
Dallos
,
H.
Lischka
,
P.
Szalay
,
R.
Shepard
, and
D. R.
Yarkony
,
J. Chem. Phys.
120
,
7330
(
2004
).
33.
H.
Lischka
,
R.
Shepard
,
I.
Shavitt
,
R.
Pitzer
,
M.
Dallos
,
T.
Müller
,
P. G.
Szalay
,
F. B.
Brown
,
R.
Alhrichs
,
H. J.
Böhm
,
A.
Chang
,
D. C.
Comeau
,
R.
Gdanitz
,
H.
Dachsel
,
C.
Erhard
,
M.
Ernzerhof
,
P.
Höchtl
,
S.
Irle
,
G.
Kedziora
,
T.
Kovar
,
V.
Parasuk
,
M.
Pepper
,
P.
Scharf
,
H.
Schiffer
,
M.
Schindler
,
M.
Schüler
, and
J.-G.
Zhao
, COLUMBUS, An ab initio Electronic Structure Program,
2003
.