Biorthonormal Formalism for Nonadiabatic Coupled Cluster Dynamics.

In coupled cluster theory, the electronic states are biorthonormal in the sense that the left states are orthonormal to the right states. Here, we present an extension of this formalism to a left and right total molecular wave function. Starting from left and right Born-Huang expansions, we derive projected Schrödinger equations for the left and right nuclear wave functions. Observables may be extracted from the resulting wave function pair using standard expressions. The formalism is shown to be invariant under electronic basis transformations, such as normalization of the electronic states. Consequently, the nonadiabatic coupling elements can be expressed with biorthonormal electronic wave functions. Calculating normalization factors that scale as full configuration interaction is not necessary, contrary to claims in the literature. For nonadiabatic nuclear dynamics, we need expressions for the derivative couplings in the biorthonormal formalism. These are derived in a Lagrangian framework.

Introduction

Nonadiabatic coupling elements account for electron–nucleus interactions that are neglected in the Born–Oppenheimer[1] (BO) approximation. These elements couple different electronic states through the nuclear kinetic energy operator. While mostly negligible in ground-state chemistry, coupling elements are required when considering molecular dynamics in excited electronic states. Excited-state dynamics often involves regions of nuclear space where electronic states are nearly or exactly degenerate, causing a breakdown of the BO separation.[2,3] Accurately describing nonadiabatic coupling elements is therefore important for reliable predictions in photochemistry.

Coupled cluster theory is one of the most accurate electronic structure methods, both for ground- and excited-state properties,[4−7] but it has not found widespread use for predicting nonadiabatic dynamics. This is primarily because standard coupled cluster methods give a nonphysical description of regions close to electronic degeneracies or conical intersections.[8−10] This issue can be traced to the method’s non-Hermiticity, which seems to imply that coupled cluster methods cannot be used for nonadiabatic dynamics. However, this is not the case. As we have shown in recent work, the method can be constrained to give a correct physical description of excited-state conical intersections while retaining the standard non-Hermitian formalism and presumably its accuracy.[11,12] These developments may lead to renewed interest in simulations of nonadiabatic dynamics that use coupled cluster theory to describe the electronic structure. While ground-state intersections are not easily treated, the method is expected to accurately describe relaxation between excited states.

Nonadiabaticity, as described by coupled cluster methods, has been considered by several authors. The first-order derivative coupling (or vector coupling) was first derived by Christiansen,[13] who applied the Z-vector substitution method[14] on a biorthonormal expression for the coupling,where (ψ̃,ψ) refers to the left and right kth electronic states and I refers to a nucleus. However, Christiansen’s paper[13] did not include an implementation of the coupling. The vector coupling was later rederived by Tajti and Szalay[15] by differentiating the corresponding m-to-n transition element of the electronic Hamiltonian. Their derivation is closely related to that given by Ichino et al.(16) for the quasidiabatic interstate coupling. Tajti and Szalay[15] also gave an implementation at the coupled cluster singles and doubles (CCSD) level.[17] These papers on the vector coupling[13,15] did not include a discussion of the nuclear Schrödinger equations in coupled cluster theory, where the coupling elements enter.

The correct formula for the vector coupling has been a subject of some controversy. Tajti and Szalay[15] argued that the biorthonormal formula in eq is incorrect. As they correctly noted, the vector coupling changes with the norm of the left and right states. A similar observation had been made in an earlier paper on the diagonal BO correction.[18] Since the vector coupling varies with the norm of the states, the full-CC vector coupling, as given by eq , is different from the full configuration interaction (full-CI) limit, where the left and right states are identical and usually normalized. They therefore suggested that normalizing the states was necessary. Furthermore, since the derivative can either act on the left or right state, they suggested using an average of the two.[15] If true, these observations are troubling: because of the normalization factors for the right states, they imply that computing the vector coupling has a computational cost that scales as full-CI. In practice, the normalization factors are therefore approximated. However, it is unfortunate if one must resort to approximations other than the truncation level of the coupled cluster method (e.g., singles and doubles in CCSD). The need for normalization factors was also assumed in the recent CCSD implementation by Faraji et al.(19)

The first main objective of the present paper is to establish that normalization is not necessary. The reason is that normalization is a special case of an invertible transformation of the electronic basis. Such transformations do not change the expansion space in the Born–Huang expansion[20] and therefore do not change the molecular wave function. In particular, the coefficients in the Born-Huang expansion—that is, the nuclear wave functions—absorb the transformation of the electronic states. The vector coupling does depend on normalization, but this quantity is not an observable. Since normalization is not necessary, the biorthonormal formula in eq is a valid option.

In a recent paper, Shamasundar[21] found, for the rovibrational Schrödinger equation, that normalization does not affect the Born–Oppenheimer product wave function. He noted that this finding should also generalize to nonadiabatic dynamics, that is, the predicted dynamics should not depend on the normalization of the underlying electronic wave functions.[21]

Derivative couplings are readily derived with the Lagrangian approach.[6,22−24] Here, we use the Lagrangian introduced for CASCI by Hohenstein[25] to derive ground-to-excited and excited-to-excited state couplings. This Lagrangian is based on an overlap whose geometrical derivatives are nonadiabatic coupling elements;[25] the approach bears some similarities to that used by Hättig et al.,[26,27] who derived transition moments as derivatives of transition moment Lagrangians. We use the Lagrangian approach[25] to derive expressions for the vector coupling as well as the second-order derivative coupling (or scalar coupling)The scalar coupling is often omitted in dynamics simulations, but its potential influence on nonadiabatic dynamics has been considered in recent years.[3]

The second main objective of the paper is to give a framework for nonadiabatic dynamics using coupled cluster methods. In particular, we argue that the biorthonormal formalism for electronic wave functions implies a biorthonormal formalism for the molecular wave function. Hence, we must determine the left and right nuclear wave functions, and the nuclear motion is described by two sets of nuclear Schrödinger equations. The result is a molecular wave function pair (Ψ̃,Ψ). Observable quantities are calculated by the usual biorthonormal formulas.

Theory

Electronic Wave Functions in Coupled Cluster Theory

In the equation of motion coupled cluster formalism, the nth (n = 0, 1, 2, ...) left and right electronic states are expressed as[28]andwhere ⟨ = δ. The projection space is defined aswhere τμ and μ with μ > 0 are the excitation operators relative to the Hartree–Fock state and ⟨μ|ν⟩ = δμν. The cluster operator is defined as where {tμ} are scalars called cluster amplitudes.

The right ground state is assumed to have the exponential formand projecting the Schrödinger equation onto {⟨μ|}μ≥0 gives the ground state equationswhere H̅ = exp(−T)H exp(T). The state amplitudes, and , are determined by making the nth energy stationary under the condition . The result is a set of eigenvalue equationswhere E is the energy of the nth state and . Here and throughout, we use bold font () to denote vectors (Xμ) and matrices (Xμν). In summary, eqs –12 are solved to determine the left and right electronic states. The reader is referred to the literature for more details.[28,29]

Born–Huang Expansion of the Total Wave Function and the Nuclear Schrödinger Equations

The Born–Huang expansion expresses the total wave function in terms of the electronic wave functions. The coefficients of the expansion define the nuclear wave functions. These nuclear wave functions are determined by inserting the expansion in the Schrödinger equation and projecting out the electronic components; the result is a set of nuclear Schrödinger equations. In coupled cluster theory, this implies a biorthonormal description of the total wave function since we can expand in both the left and right electronic wave functions ( and ψ). Hence, we have a left and a right total wave function given by the Born–Huang expansionswith associated left and right nuclear wave functions χ̃ and χ, andwhere we have assumed biorthonormal electronic states in the third equality: = δ. In eqs and 14, the electronic and nuclear coordinates are denoted by and , respectively, and time by t. Expectation values are defined through the standard expression[28,29]

To derive the equations for the nuclear wave functions, one normally projects the total Schrödinger equation on the electronic basis. In this respect, a biorthonormal description is advantageous; for practical coupled cluster models, where the excitation space is truncated to some excitation order, projection of the right Schrödinger equation is done onto the left electronic basis, leading to computationally tractable expressions that scale as expected for the given model (e.g., for CCSD, where N is the size of the system).

By inserting the Ψ in eq into the time-dependent Schrödinger equationand projecting it onto the left electronic basis, we get a coupled set of equations for the right nuclear wave functions χ. These nuclear Schrödinger equations can be expressed aswhere we have suppressed the and t dependence for readability. The nonadiabatic coupling vectors in eq are given in the biorthonormal basis

In analogous fashion, we derive the nuclear Schrödinger equations for the left nuclear wave functions from the complex conjugated Schrödinger equation

Inserting eq into 21 and projecting onto the right electronic basis leads to

The nuclear Schrödinger equations can be expressed in the more compact matrix notationwhere is a diagonal matrix with the electronic energies on the diagonal, is the identity matrix, χ is a vector containing the right nuclear wave functions, and and are matrices consisting of the scalar and vector couplings of the Ith nucleus, respectively. The quantities with a tilde are similarly defined.

This matrix notation has been used to illuminate some relations to gauge theories in the nuclear Schrödinger equations; Pacher et al.(30) found that the vector coupling can be seen to serve a role analogous to the vector potential in electromagnetism. In the present work, it serves as a useful notation for dealing with basis transformations and the vector algebra needed to demonstrate invariance under such transformations.

Basis Invariance and the Special Case of Norm Invariance

It has been suggested that normalizing the electronic states is necessary when calculating nonadiabatic coupling elements.[15] The reason is that the left and right states are biorthonormal in the coupled cluster theory. Compared to the couplings in the full-CI theory, where the states are normalized, the full coupled cluster limit is “incorrect” because the couplings depend on the geometry-dependent normalization constants. While this suggests that one should normalize the states, doing so is not feasible. The computational cost of the normalization factor scales as full-CI for the right electronic states[15]Since one cannot evaluate N in general, some have suggested N = (N)−1 or N = N = 1 as alternatives. The former gives the full-CI limit, while the latter is equivalent to assuming the standard biorthonormality.[15,19]

Biorthonormality is not an issue from the point of view of nonadiabatic dynamics. In fact, normalizing the electronic states is a special case of basis transformation of the electronic basis. As such, the Born–Huang expansion and the nuclear Schrödinger equations are equivalent in the transformed and untransformed bases. Changes in the electronic basis are absorbed in the expansion coefficients, that is, the nuclear wave functions. In the case of normalization, the right electronic wave functions are divided by N, while the right nuclear wave functions are multiplied by N. The total wave function does not change.

More precisely, consider invertible transformations of the left and right electronic bases. In vector notation, these transformations can be expressed aswhere the matrices and are assumed to be smooth invertible matrix functions of the nuclear coordinates. For notational simplicity, we have let the left and right wave function vectors be row vectors. Transformed quantities are denoted by a prime. In the transformed basis, the total left and right wave functions have the Born–Huang expansions

We wish to show that the wave function in the transformed basis is identical to that obtained in the untransformed basis, that is, Ψ′ = Ψ and Ψ̃′ = Ψ̃. The conclusion that follows is that the choice of electronic basis does not change the predictions of the theory. In other words, it is perfectly appropriate to use the standard[28,29] biorthonormal description.

Before proceeding, we define some notation. In the transformed basis, we have to account for the nonunit overlap of the electronic wave functions. Hence, when projecting the time-dependent Schrödinger equation onto the electronic basis, we get electronic overlap matrix elements. In particularSimilarly, the electronic Hamiltonian matrix is not necessarily diagonal

We show the equivalence for the right wave functions. The proof for the left wave function is identical. Following the standard procedure, we now insert the transformed wave function in eq into the Schrödinger equation and project onto the transformed left electronic wave functions. The result is the right nuclear Schrödinger equationIf the total wave function is invariant and ψ′ = ψ, then we must have nuclear wave functions that cancel the transformation of the electronic wave functions, meaning that

Indeed, with χ′ as given in eq , we have

Let us confirm that eq is in fact a solution to the transformed nuclear Schrödinger equation in eq . We begin by relating the old and new nonadiabatic coupling terms. The gradients of the electronic wave functions transform asHence, the vector couplings can be written as

In more compact matrix notation, we have

Similarly, the Laplacian of the electronic wave functions transforms asimplying that the scalar couplings transform as

The gradient and Laplacian of χ′ is derived in the same way as for the electronic states, givingThus, we have the following contributions on the right-hand side of the nuclear Schrödinger equation

Though somewhat involved, most of the terms cancel when added together. In fact, since

we can writeIn other words, with χ′ = –1χ, the right nuclear Schrödinger equation simplifies to

which, upon premultiplication by –†, is seen to be equivalent to the original right nuclear Schrödinger equation in eq .

Since all of the steps we have made are reversible, we have shown that χ is a solution to the untransformed nuclear Schrödinger equation if and only if χ′ is a solution to the transformed Schrödinger equation. The total right wave function is therefore invariant with respect to transformations of the electronic basis, Ψ′ = Ψ. One consequence of basis invariance is that the nonadiabatic couplings can be expressed in the standard biorthonormal formalism.

In the next section, we derive these in a Lagrangian framework. The couplings have contributions that arise from the geometry-dependence of the many-body operators, where a choice of orbital connection is necessary.[13,31] These terms are described in Appendix A.

Nonadiabatic Coupled Cluster Couplings in a Lagrangian Formalism

To obtain a Lagrangian for the vector coupling in CASCI, Hohenstein[25] defined a partially frozen overlap whose first derivatives are identical to components of the vector coupling. In coupled cluster methods, this overlap can be expressed asin terms of which we haveandWe let i ∈ I mean that x is one of the three coordinates of the Ith nucleus (x, y, or z). Note that depends on 0. We suppress this dependency to keep the notation simple.

The overlap is expressed in terms of coupled cluster wave functions, which depend not only on but also on a set of wave function parameters λ (which themselves depend on ). Written out in terms of the parameters, we havewhereand

The κ operator accounts for orbital rotations, meaning changes in the Hartree–Fock orbitals, where, by definition, we have κ(0) = 0. Following the standard recipe, we add the equations (denoted by ) that determine the parameters as constraints with associated Lagrangian multipliers (denoted by γ)where λ and γ are determined for every by the stationarity conditions

The derivatives of this Lagrangian are identical to the derivatives of the frozen overlap (since ). One advantage of the Lagrangian formalism is that it automatically incorporates the 2n + 1 and 2n + 2 rules for λ and γ, respectively.[6] In particular, we havewhere the final equality follows from stationarity, see eqs and 60. Denoting partial derivatives with respect to geometrical coordinates aswe can write

Furthermore, if we letthen the scalar coupling can be expressed as (see Appendix B)Clearly, and G are similar in complexity to the energy gradient and Hessian. However, G is somewhat simpler than the energy Hessian because the first derivatives of the parameters (λ() can be considered one at a time.

To proceed, we must define the Lagrangian in detail. The conditions include all equations that must be solved to evaluate the overlap . These are (a) the Hartree–Fock equations, (b) the amplitude equations, and (c) the eigenvalue equations for the right state amplitudes. Written out in full, we havewhere we have introduced multipliers associated with the different sets of equations, κ̅, ζ, as well as β and E̅. We have also introduced the Brillouin conditionwhereFurthermore, the similarity transformed the Hamiltonian in Ω and is given by

and the nth electronic energy is defined as

With defined, we can now consider the equations for the zeroth order multipliers. These are determined from the zeroth order terms of the λ stationarity, eq . To keep our notation simple, we will denote the zeroth order terms as γ(0) ≡ γ and λ(0) ≡ λ, where it should be understood from context when these are γ and λ evaluated at 0. Differentiation with respect to the state parameters givesTo solve this equation, we note that if we let

the equation for β becomesThus, we have

Next, we consider stationarity with respect to . This can be expressed aswhereand where we have introduced the notation

Finally, we have stationarity with respect to κ, which can be written aswhereand

With the zeroth order multipliers determined, we can derive the expression for the vector coupling. By partially differentiating , we find thatwhereand where quantities at 0 are denoted as y(0) ≡ y (e.g., we denote T(0) as T). Here, we have assumed the natural connection, for which there are no contributions to that originate from the many-body operators (see Appendix A and refs (13) and (31)).

The vector coupling given in eq has also been identified by other authors. It was derived by Christiansen,[13] who assumed biorthonormality and used Z-vector substitution[14] on the expression for the vector coupling. Tajti and Szalay[15] identified the same expression indirectly using Z-vector substitution on derivatives of Hamiltonian transition elements. However, they also argued[15] that the coupling should not be given by eq but rather be averaged and expressed with normalized states. As we have shown, eq is a valid choice due to norm invariance and represents the vector coupling in the right nuclear Schrödinger equations. For the left Schrödinger equations, we can make use of the identity

Before moving on to the scalar coupling, we note that although the Z-vector substitution method is equivalent to the Lagrangian technique, the latter method gives, in our opinion, an especially elegant way of deriving the coupling elements.

For the scalar coupling, we must determine the first derivatives of the parameters. Equations for these are obtained as the first-order terms of the multiplier stationarity conditions. In the case of , we havewhereIn the case of κ, we similarly have

whereThe biorthonormality condition implies

while the eigenvalue condition impliesHere, we have definedwhere

With the derivatives of the parameters determined, let us next consider fα( and Hαβ, see eq . Recall that the α and β indices refer to the parameters λα and λβ. The gradient is given by the zeroth order equations for the multipliers, that is, eqs , 75, and 79, with λ = λ0 but allowing for ≠ 0. Partially differentiating these terms with respect to x gives (. The blocks of the ∑αfα(λα( contributions to G can be writtenandwhere repeated indices implies summation. For contributions to G involving Hαβ we have, for terms involving and κ,as well asand

Next, we have terms involving right state and the cluster amplitudes and orbital rotations, i.e.Finally, we have the partial derivative of the Lagrangian, which can be written as

Written in compact notation, the scalar coupling iswhere we have letas well asThe final term in eq arises from the orbital connection (see Appendix A).

Throughout the derivations above, we have considered the off-diagonal coupling elements (m ≠ n). The diagonal terms can be derived from the Lagrangianwhich gives the slightly different stationarity conditionHere, we again select E̅ to make the first term vanish, giving

Other than this change, the derivation of the scalar coupling is virtually unchanged. Terms involving differentiation of has the left state in the bracket instead of (e.g., in the stationarity conditions for the zeroth order multipliers). In particular, the expression in eq is valid with m = n.

Unlike for the vector coupling, there is no convenient relationship between G and G̃. To derive the latter quantity, we may consider the LagrangianwhereThe stationarity then gives

from which we again have and thus

The equations for the zeroth order multipliers are derived as before, with the result that the multipliers change their sign, thus giving the result in eq for the vector coupling. For the derivative of the parameters, we have the same equations for ( and κ(. For the derivative of , we must solve the equationwhich is analogous to eq . In contributions involving in G, we obtain similar expressions involving in the case of G̃. The end result iswith

Finally, G̃ is obtained in a manner similar to G, see eq and the surrounding text.

This concludes our derivation of the coupled cluster scalar coupling. To the best of our knowledge, equations for this coupling (with m ≠ n) have not been presented in the literature before. Diagonal terms were also considered by Gauss et al.(18) from a different starting point.

Concluding Remarks

The norm of the electronic states changes the value of nonadiabatic coupling elements but does not change the molecular wave function. The biorthonormal formula assumed by Christiansen[13] is therefore a valid choice for nonadiabatic dynamics using coupled cluster methods. More generally, we have shown that the total wave function is invariant with respect to smooth and invertible transformations of the electronic basis. Of course, the biorthonormal couplings are not directly comparable to the coupling elements of an Hermitian method with normalized states, such as CI or full-CI. However, this reflects the basis-dependence of the couplings and not the validity of the biorthonormal formalism.

We therefore derive a set of nuclear Schrödinger equations assuming biorthonormal projection onto the electronic basis. Combined with expressions derived for the vector and scalar couplings, these nuclear Schrödinger equations serve as a starting point for the application of coupled cluster methods in simulations of nonadiabatic dynamics.

Our derivations have been restricted to the standard coupled cluster theory. However, the Lagrangian formalism is easily extended to similarity constrained coupled cluster methods,[11,12] which are suited to describe relaxation through a conical intersection between excited states. The application to ground-state intersections is less straightforward but may be accessible with approaches that use a different reference than the closed-shell Hartree–Fock state.[32]