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,  (1) where (ψ̃k,ψ k ) 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 Schrodinger 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 1 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 1, 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 1 is a valid option.
In a recent paper, Shamasundar 21 found, for the rovibrational Schrodinger 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 Haẗtig 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 Schrodinger 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 and where ⟨ψ̃m ψ ⟩ n = δ mn . The projection space is defined as where τ μ and τμ with μ > 0 are the excitation operators relative to the Hartree−Fock state ⟩ HF and ⟨μ|ν⟩ = δ μν . The cluster operator is defined as where {t μ } are scalars called cluster amplitudes. The right ground state is assumed to have the exponential form  (12) where E n is the energy of the nth state and μ ν ̅ = ⟨ ̅ ⟩ μν H .
Here and throughout, we use bold font (X) to denote vectors (X μ ) and matrices (X μν ). In summary, eqs 10−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 Schrodinger 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 Schrodinger equation and projecting out the electronic components; the result is a set of nuclear Schrodinger 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 (ψ̃n and ψ n ). Hence, we have a left and a right total wave function given by the Born−Huang expansions n n n (14) with associated left and right nuclear wave functions χñ and χ n , and where we have assumed biorthonormal electronic states in the third equality: ψ ψ ⟨̃⟩ m n = δ mn . In eqs 13 and 14, the electronic and nuclear coordinates are denoted by r and R, 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 Schrodinger 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 Schrodinger equation is done onto the left electronic basis, leading to computationally tractable expressions that scale as expected for the given model (e.g., N ( ) 6 for CCSD, where N is the size of the system). By inserting the Ψ in eq 13 into the time-dependent Schrodinger equation and projecting it onto the left electronic basis, we get a coupled set of equations for the right nuclear wave functions χ n . These nuclear Schrodinger equations can be expressed as  (18) where we have suppressed the R and t dependence for readability. The nonadiabatic coupling vectors in eq 18 are given in the biorthonormal basis In analogous fashion, we derive the nuclear Schrodinger equations for the left nuclear wave functions from the complex conjugated Schrodinger equation Inserting eq 14 into 21 and projecting onto the right electronic basis leads to  (24) The nuclear Schrodinger equations can be expressed in the more compact matrix notation where E is a diagonal matrix with the electronic energies on the diagonal, I is the identity matrix, χ is a vector containing the right nuclear wave functions, and G I and F I 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 Schrodinger 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 R n in general, some have suggested N R n = (N L n ) −1 or N R n = N L 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 Schrodinger 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 R n , while the right nuclear wave functions are multiplied by N R 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 as 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 Schrodinger equation onto the electronic basis, we get electronic overlap matrix elements. In particular 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 31 into the Schrodinger equation and project onto the transformed left electronic wave functions. The result is the right nuclear Schrodinger equation If the total wave function is invariant and ψ′ = ψM, then we must have nuclear wave functions that cancel the transformation of the electronic wave functions, meaning that Indeed, with χ′ as given in eq 36, we have ∑ ∑ ∑ implying that the scalar couplings transform as The gradient and Laplacian of χ′ is derived in the same way as for the electronic states, giving we can write  I  I  I I   I  I  I I  I  I   I  II   I  I  I I   2  In other words, with χ′ = M −1 χ, the right nuclear Schrodinger equation simplifies to ( 2 ) which, upon premultiplication by N − † , is seen to be equivalent to the original right nuclear Schrodinger equation in eq 25.
Since all of the steps we have made are reversible, we have shown that χ is a solution to the untransformed nuclear Schrodinger equation if and only if χ′ is a solution to the transformed Schrodinger 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 as We let i ∈ I mean that x i is one of the three coordinates of the Ith nucleus (x, y, or z). Note that x ( ) mn depends on x 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 x but also on a set of wave function parameters λ (which themselves depend on x). Written out in terms of the parameters, we have The κ operator accounts for orbital rotations, meaning changes in the Hartree−Fock orbitals, where, by definition, we have κ(x 0 ) = 0. Following the standard recipe, we add the equations (denoted by mn ) that determine the parameters as constraints with associated Lagrangian multipliers (denoted by γ) where λ and γ are determined for every x by the stationarity conditions The derivatives of this Lagrangian are identical to the derivatives of the frozen overlap (since = 0 mn ). 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 have Furthermore, if we let Clearly, F mn I and G mn I are similar in complexity to the energy gradient and Hessian. However, G mn I is somewhat simpler than the energy Hessian because the first derivatives of the parameters (λ (i) ) can be considered one at a time.
To proceed, we must define the Lagrangian mn in detail. The conditions mn include all equations that must be solved to evaluate the overlap mn . These are (a) the Hartree−Fock equations, (b) the amplitude equations, and (c) the eigenvalue where we have introduced multipliers associated with the different sets of equations, κ ̅ , ζ, as well as β n and E̅ n . We have also introduced the Brillouin condition Furthermore, the similarity transformed the Hamiltonian in Ω and ̅ is given by ) and the nth electronic energy is defined as With mn defined, we can now consider the equations for the zeroth order multipliers. These are determined from the zeroth order terms of the λ stationarity, eq 60. 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 x 0 . Differentiation with respect to the state parameters gives Next, we consider stationarity with respect to t. This can be expressed as and where we have introduced the notation Finally, we have stationarity with respect to κ, which can be written as and where quantities at x 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 F mn I that originate from the many-body operators (see Appendix A and refs 13 and 31).
The vector coupling given in eq 82 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 82 but rather be averaged and expressed with normalized states. As we have shown, eq 82 is a valid choice due to norm invariance and represents the vector coupling in the right nuclear Schrodinger equations. For the left Schrodinger 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 t, we have where In the case of κ, we similarly have The biorthonormality condition implies while the eigenvalue condition implies Here, we have defined With the derivatives of the parameters determined, let us next consider f α (i) and H αβ , see eq 64. Recall that the α and β indices refer to the parameters λ α and λ β . The gradient f is given by the zeroth order equations for the multipliers, that is, eqs 71, 75, and 79, with λ = λ 0 but allowing for x ≠ x 0 . Partially differentiating these terms with respect to x i gives f (i) . The blocks of the ∑ α f α (i) λ α (i) contributions to G mn I can be written Next, we have terms involving right state and the cluster amplitudes and orbital rotations, i.e.
Written in compact notation, the scalar coupling is where we have let The final term in eq 108 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 Lagrangian which gives the slightly different n stationarity condition Here, we again select E̅ n to make the first term vanish, giving Other than this change, the derivation of the scalar coupling is virtually unchanged. Terms involving differentiation of nn has the left state ⟨ n in the bracket instead of ⟨ m (e.g., in the stationarity conditions for the zeroth order multipliers). In particular, the expression in eq 102 is valid with m = n.
Unlike for the vector coupling, there is no convenient relationship between G mn I and G̃m n I . To derive the latter quantity, we may consider the Lagrangian The m stationarity then gives 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 84 for the vector coupling. For the derivative of the parameters, we have the same equations for t (i) and κ (i) . For the derivative of m , we must solve the equation Finally, G̃n n I is obtained in a manner similar to G nn I , see eq 109 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 Schrodinger equations assuming biorthonormal projection onto the electronic basis. Combined with expressions derived for the vector and scalar couplings, these nuclear Schrodinger 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 connection for nonadiabatic coupling elements was also discussed by Christiansen. 13 To evaluate derivatives at x 0 , we relate the basis at x 0 to some basis at x = x 0 + Δx. Suppose the molecular orbitals (MOs) at x 0 are where C αm are orbital coefficients and χ α are atomic orbitals. The unmodified MOs (UMOs) are defined by freezing the orbital coefficients The UMOs are not orthonormal, however where b(x) is the anti-Hermitian operator with b pq (x) defined such that U(x) = exp(−b(x)). The many-body operators can be expanded as λ