Calculations of non-adiabatic couplings within equation-of-motion coupled-cluster framework: Theory, implementation, and validation against multi-reference methods

Transitions between electronic states play a key role in excited-state processes. They are essential in photochemistry, photobiology, and the operation of photovoltaic materials. These transitions are responsible for radiationless relaxation, a process in which electronic energy is converted into nuclear motions. Radiationless relaxation protects living organisms from the damaging effects of UV radiation by converting energy imparted by solar photons into heat and efficiently quenching unwanted excited-state chemistry.¹ Radiationless relaxation is a key step in vision—it allows photoexcited rhodopsin’s chromophore to return to the initial state, upon completing the isomerization process, so the cycle can be repeated over and over again. In other instances, radiationless relaxation can be viewed as a nuisance. For example, in imaging applications that use fluorescent tags, radiationless relaxations compete with fluorescence (radiative relaxation), leading to reduced optical output.² In solar energy harvesting, radiationless relaxation is also often a parasitic process as it competes with charge separation or reactions producing solar fuels.³ However, non-adiabatic transitions can be exploited to increase the efficiency of solar cells via multi-exciton generation. In organic photovoltaics, this process is called singlet fission;^4,5 it involves a non-adiabatic transition between the initially excited bright state and a dark state of multi-exciton character.⁶ 

Within the Born-Oppenheimer separation of nuclear (R) and electronic (r) coordinates, the exact molecular wave function is represented by the following ansatz:

where Φ_I(r; R) and $ξIK(R)$ denote electronic and nuclear wave functions, respectively. The electronic wave functions are defined as the solutions of the electronic Schrödinger equation,

where the Hamiltonian, which includes all terms (electron kinetic energy T_e and Coulomb interactions between all charged particles) except nuclear kinetic energy T_n, depends on nuclear positions via V_en. Equation (2) is solved at each nuclear configuration, giving rise to eigenstates and eigenenergies depending parametrically on R; at each R, the electronic wave functions Φ_I form an orthonormal basis. Ansatz (1) gives rise to the following equations for nuclear wave functions:

where A denotes a particular nucleus and $EK$ is the total energy of wave function (1). Within the Born-Oppenheimer (or adiabatic) approximation, the terms on the right-hand side of Eq. (4) are neglected giving rise to separable wave function, $ΦIξIK$⁠. In this case, the nuclear motions are confined to a single potential energy surface, E_K(R), and nuclear wave functions corresponding to each electronic state are independent (or uncoupled) from nuclear wave functions corresponding to other electronic states. Physically, this means that electrons adjust to moving nuclei’s positions instantly so that the electronic state of the system never changes. Electronic transitions are only possible beyond adiabatic approximation. They arise due to the terms on the right-hand side of Eq. (4), which couple nuclear motions on different potential energy surfaces (PESs). These terms originate in the parametric dependence of the electronic wave functions on the nuclear coordinates and are called derivative couplings. They become large when the changes in electronic wave functions due to nuclear motions are large, which is quantified by the derivative with respect to nuclear coordinates. These terms lead to electronic inertia so that electrons are no longer following nuclei adiabatically; rather, upon some nuclear motions, electrons can remain frozen at their configuration leading to the change of electronic state. Of the two derivative terms, only the first term contributes to the diagonal Born-Oppenheimer energy correction, which is computed when high-accuracy thermochemical data are desired.^7,8 This term may cause discontinuous behavior in calculations of non-adiabatic dynamics in adiabatic representation⁹ and is omitted in the original surface-hopping model.¹⁰ Recently, Meek and Levine have shown that this term can be accurately accounted for within multiple-spawning calculations¹¹ on adiabatic surfaces by using locally diabatized Gaussian basis functions.¹² Interestingly, numeric simulations of non-adiabatic relaxation in ethylene have shown that the results are nearly identical whether this term is fully accounted for or completely omitted.¹² 

In this paper, we focus on the term called derivative coupling or non-adiabatic coupling (NAC),

As clearly seen from Eq. (5), NAC is a 3N vector, where N is the number of nuclei. d_IJ couples nuclear motions in electronic states I and J (for real-valued wave functions, the diagonal NAC elements vanish, d_II = 0, by virtue of the normalization condition). As one can see from Eq. (4), the total magnitude of the coupling is given by the scalar product of d_IJ and the nuclear momentum, ∇_Rξ. Thus, the probability of electronic transitions depends not only on the electronic coupling d_IJ but also on how fast the nuclei are moving as well as the direction of their motion. Obviously, for infinitesimally slow nuclei, the coupling vanishes and the adiabatic approximation is restored. The symmetry of d_IJ is determined by the symmetries of Φ_I and Φ_J and is equal to the product of the two irreps. Consequently, the two states are only coupled by the nuclear motions of that symmetry (otherwise, the scalar product vanishes). For example, in a $C2v$ molecule, electronic states from the A₂ and B₁ irreps are coupled by B₂ vibrations (and the symmetry of d_IJ vector is B₂). An important aspect of non-adiabatic dynamics in polyatomic systems is that, in contrast to the one-dimensional case, one can distinguish between molecular motions inducing the transition (those along d_IJ) and motions that modulate the magnitude of d_IJ (those along which the wave-function composition changes significantly); for concrete examples, see Refs. 6 and 13.

For the exact wave functions, NAC vector d_IJ can be expressed as a matrix element of the derivative of the electronic Hamiltonian, H^el,

Thus, one can expect large values of NAC when (i) the electronic Hamiltonian depends strongly on nuclear positions (large h_IJ) and (ii) the energy gap between the two states is small (as in the regions of PESs’ near-degeneracy). The two vectors, h_IJ and d_IJ, are parallel. Vector h_IJ, called the non-adiabatic force matrix element, can be described as an interstate generalization of the nuclear gradient,

The second equality holds only when the Hellman-Feynman condition is satisfied.

This connection to the nuclear gradient can be exploited for practical calculations of NACs based on the following observation.¹⁴ Consider two electronic states, Φ_I and Φ_J, and a fictitious summed state Φ_(I+J) ≡ Φ_I + Φ_J. The energy (defined as an expectation value) of Φ_(I+J) equals E_I + E_J, but Φ_(I+J) is not an eigenstate—thus, the Hellman-Feynman theorem does not hold for it. One can define (and compute) Hellman-Feynman gradients for all three states, which are G_K ≡ ⟨Φ_K|∇_RH|Φ_K⟩, K = I, J, I + J. At the same time,

Thus

and, consequently,

Equation (10) provides an insight into the nature of non-adiabatic force matrix element. It is the non-Hellman-Feynman part of the gradient of the summed state, i.e., the difference between the expectation value of ∇_RH and the sum of the gradients of the two states. When the two states belong to different irreps, this is the part of the summed-state gradient which is not fully symmetric but has symmetry of the product of the two irreps.

Because the magnitude of NACs becomes large when the energy gap between the two electronic states is small, the F − 2 dimensional seams along which the surfaces cross are of special significance in non-adiabatic processes (F denotes the number of internal degrees of freedom; for a nonlinear molecule, F = 3N − 6). The minimum-energy crossing points^15–17 along these seams (called MECPs) play a role akin to transition states in transition-state theory (TST)—their energetic accessibility is one of the factors controlling the rate of non-adiabatic transitions. In the non-adiabatic extension of TST, MECP structures are used to compute the density of states, electronic couplings, and Arrhenius factors that together define the rate of non-adiabatic reaction flow. The degeneracy seams are also called conical intersections, as the topology of the two PESs around such seams has a conical shape. The conical shape arises because the degeneracy between the two states, I and J, is lifted along the following two coordinates,^17,18 h_IJ [defined by Eq. (29)] and g_IJ ≡ G_I − G_J. These two orthogonal vectors give rise to the g − h plane defining the conical intersection (Fig. 1). Thus, h is needed for characterizing the conical intersection and for computing the rate of non-adiabatic transitions (we note that efficient algorithms^19,20 for locating MECPs also use h).

What is needed for computing MECPs and NACs by electronic structure methods? First, the method should be capable of describing multiple electronic states on the same footing, so the relative state energies and, consequently, the degeneracy seams are correctly described. Second, the method should be able to handle multi-configurational wave functions in a sufficiently flexible way, correctly reproducing the changes in the adiabatic wave functions upon changes in nuclear geometries. Third, the expressions for the exact states should be reformulated for approximate wave functions using techniques similar to those used in analytic gradient and response theory. Pioneering studies of conical intersections and non-adiabatic phenomena employed multi-reference configuration interaction (MRCI) method.^21–25 The formalism and computer implementations for the calculation of NACs using complete active space self-consistent field (CASSCF) wave functions^26,27 and within multi-state (CASSCF corrected by second-oder perturbation theory) framework^28,29 have also been reported. Calculations of MECPs and NACs by using more economical methods, configuration interaction singles (CISs) and time-dependent density functional theory (TDDFT), have been presented.^30–34 These single-reference methods use linear parameterization of target states and are, therefore, capable of describing conical intersections between excited states in the vicinity of the Franck-Condon region, where the reference state retains its single-configurational character. However, the intersections between the ground and excited states cannot be treated by these approaches, owing to single-determinantal reference state.³⁵ This issue is addressed by the spin-flip (SF) approach in which a high-spin reference is used, and the target states (both the ground state and excited states) are described as spin-flipping excitations.^36–39 Recently, calculations of NACs and MECPs using SF-CIS and SF-DFT have been reported^34,40,41 and used in ab initio surface-hopping dynamics.⁴² 

Equation-of-motion coupled-cluster (EOM-CC) methods^43–45 share many common features with CI. The target states are described by the CI-like excitation operators (R) from the CC reference state,

where $0$ is the reference determinant and amplitudes of T satisfy CC equations for the reference state. Amplitudes R are found by solving a non-Hermitian eigenproblem,

Depending on the specific choice of the reference and excitation operators, different types of target states can be accessed, as illustrated in Fig. 2. Most commonly used variant is EOM for excitation energies in which the spin- and particle-conserving excitation operators are used, and the reference corresponds to the ground state of the system.^46–48 EOM operators that change the number of electrons (or their spins) open access to various open-shell wave functions.^36,49–53 EOM-CC is capable of describing many types of multi-configurational wave functions within single-reference formalism.^43,54 Because of its multi-state nature, EOM-CC can treat iterations and crossings between multiple interacting states. Thus, EOM-CC is an attractive platform for computing MECPs and NACs.^55–62 

Formulation of NACs within EOM-CC theory have been described for EE (excitation energy) and IP (ionization potential) variants.^14,55 Here we extend the theory by including the EOM-EA (EOM for electron affinity) method. We also extend the list of examples for EOM-EE and EOM-IP methods and provide detailed comparisons with CASSCF and MRCI. Our focus is on NACs rather than MECPs, as the former has not received as much attention as the latter. The outline of this paper is as follows. Section II presents the formalism for computing NAC forces for inexact states and within non-Hermitian framework, as needed for EOM-CC wave functions. Section III summarizes computational details. The benchmark results and analysis are presented in Sec. IV.

For inexact states, calculations of properties, such as nuclear gradients, require extra steps, relative to simply computing an expectation value of the appropriate derivative of the Hamiltonian using approximate wave functions. Additional terms (often referred to as “Pulay terms”), which include derivatives of the wave function parameters, arise because the Hellman-Feynman theorem is not satisfied. The efficient evaluation of these terms can be derived using the Lagrangian approach.⁶³ For example, in CC/EOM-CC theory, the analytic gradient can be formulated as a contraction of derivatives of the Hamiltonian with one- and two-particle density matrices which include an unrelaxed part (the equivalent of the expectation-value calculation) and amplitude and orbital-response parts.^50,64–67 The latter two contributions are computed by solving auxiliary response equations. This formalism has been extended by Stanton and co-workers⁵⁵ for calculating NAC forces, which they refer to as quasidiabatic coupling strength; see Eq. (57) in Ref. 55,

where x denotes a particular Cartesian coordinate with respect to which the derivative is taken. Here the first term corresponds to the full derivative of the Hamiltonian (i.e., including the derivatives of the molecular orbitals), and the second two terms include EOM-CC amplitude response contributions. The full derivation of the EOM-CC NAC force and the comparison with quasidiabatic couplings λ defined in Ref. 55 is given in  Appendix A.

In the CI formalism, it is convenient to cast the gradient and NAC force expressions in the generalized Hellman-Feynman form,^31,34

where $H$ denotes the projected Hamiltonian (e.g., in the case of CIS, it is the Hamiltonian projected into the space of singly excited determinants) and superscript x denotes the full derivative (which includes the derivative of the Hamiltonian as well as the derivative of the projector operator). Given the CI-like form of the EOM-CC equations, the expression for the NAC force in EOM-CC theory can be cast in the same form,

where $H¯$ is the similarity transformed Hamiltonian projected onto the space of the reference, singly, and doubly excited determinants, and $hIJx$ denotes the x component of h_IJ. By taking the derivative, one arrives to formally exact h_IJ, which are identical to a calculation based on the non-Hermitian generalization of Eq. (10).

For approximate wave functions, the connection between the NAC vector (d_IJ) and the NAC force (h_IJ) is not as simple as in Eq. (6). Applying the definition given by Eq. (5) and breaking the derivative into two terms, $∇R=∇RC+∇Rϕ$⁠, where ϕ denotes the derivative of molecular and atomic orbitals and C denotes all the rest, one arrives at

The orbital-derivative part, $dIJϕ$⁠, can be evaluated, giving rise to the “proper” d_IJ, which has been done within MRCI^21,23,24 and EOM-CC frameworks.¹⁴ The expression for $dIJϕ$ is given in  Appendix B. Interestingly, full d_IJ computed for approximate wave functions is not translationally invariant. That is, one can induce non-adiabatic transitions by simply moving the entire molecule around, which is unphysical and problematic for dynamics simulations. Subotnik and co-workers³¹ have developed a correction restoring translational invariance of d_IJ within the CIS framework; the correction was called ETFs (electronic translation factors). Later, Zhang and Herbert have shown that the correction is equal to the orbital contribution part and cancels it out exactly.³⁴ By virtue of Eq. (10), it is clear that h_IJ and, consequently, $dIJC$ are translationally invariant, as they are expressed in terms of expectation values and gradients. Thus, the unphysical violation of the translational invariance can only come from $dIJϕ$⁠, and by omitting this unphysical term, one naturally arrives to the ETF-corrected NACs. Therefore, in the present study, we focus on the translationally invariant NAC,

Another rational⁵⁵ for omitting $dIJϕ$⁠, as has been done in the definition of quasi-diabatic couplings,⁵⁵ is that all important couplings leading to rapid variation of the wave function come from the diagonalization, i.e., the EOM or CI part.

The theory of calculations of NAC within the EOM-CC framework has been described in detail in Refs. 14, 55, and 68. Here we employ the approach of Tajti and Szalay who generalized Eq. (10) to the EOM-CC theory. Since EOM-CC theory is non-Hermitian, the left and right EOM states are not the conjugate of each other but form a biorthogonal set,

The norms of left or right EOM states are arbitrary, which leads to an ambiguity in calculations of the interstate matrix elements. The first (minor) issue is that IJ matrix elements are not equal to JI ones,

The second issue is that the magnitudes of these elements depend on the arbitrary normalization. One possible solution is to consider geometric average,⁴⁸ 

While the signs of A_IJ and A_JI are arbitrary, they are not independent: biorthonormality condition requires that if the sign of right EOM state I is flipped, then the sign of the left EOM state I must be flipped as well; consequently, both A_IJ and A_JI change signs. Thus, one can use the sign of either A_IJ or A_JI to assign a sign to A_IJ.

An equally justified approach of defining A_IJ would be to consider an arithmetic average,¹⁴ 

where weights W_K depend on the norms of the left and right EOM states,¹⁴ $NIL$ and $NJR$⁠,

The definition and the expressions of the norms are given in  Appendix C. Tajti and Szalay¹⁴ pointed out that in the limit of exact solution, the EOM values of interstate matrix elements defined by Eq. (26) should agree with full configuration interaction (FCI) ones.

We note that the sign of the NAC vector (or any other interstate matrix element) is not uniquely defined because the phases of adiabatic states are not defined by the electronic Schödinger equation, i.e., flipping the sign of eigenstate Φ_I(r; R) does not affect Eq. (2). The sign issue is related to a well-known geometric phase effect,⁶⁹ which requires proper handling of the phase of the electronic wave function when solving the nuclear problem for the overall wave function defined by Eq. (1).

Tajti and Szalay have shown that the EOM NAC force defined by using Eq. (26) can be computed from the gradients of the two states and a fictitious summed state,

where left and right amplitudes for the summed state are

Here we introduced normalization factors $12$ for convenience so that the separable part of the respective density matrices (reference contribution to G_K’s) is handled correctly by the regular EOM-CC code (we note that one can introduce weights W_K when computing the right summed state and use simple summation of the left amplitudes; the results are numerically identical). This approach, which allows for a straightforward evaluation of the EOM-CC NACs by trivial modifications of the analytic gradient code,⁶⁷ is exploited in this work.⁷⁰ Our benchmark calculations showed that the computed W_K are always nearly equal to one and have a minute effect on the computed NAC forces: for all considered cases, the effect on the NAC was less than 10⁻⁶. Thus, for the sake of simplicity in our final implementation, we use W_I = W_J = 1. All reported values use these unit weights.

In the context of conical intersections, the non-Hermiticity of EOM-CC theory might be a concern.⁷¹ As pointed out by Köhn and Tajti, a straightforward extension of Wigner-Neumann derivation⁷² to non-symmetric model Hamiltonians suggests that the dimensionality of a conical intersection between the states of the same symmetry becomes F-3 instead of F-2. That means that the apex of the cone in Fig. 1 becomes a cylinder. They illustrated this troublesome behavior by the EOM-EE calculations of excited states in formaldehyde, where in a small region around conical intersection, EOM-EE-CCSD roots become complex (with real parts of the energies degenerate). Yet, the exact EOM-CC theory, which is equivalent to FCI, describes the intersections correctly. This paradox has been recently explained by Koch and co-workers;⁷³ these authors also suggested a possible solution restoring the correct topology of conical intersections, as in Fig. 1, in inexact EOM-CC theories.⁷⁴ 

MRCI calculations were performed by Columbus.⁷⁵ EOM-CC calculations were performed with Q-Chem.^76,77 The Cfour calculations were performed using the modified version, as described in  Appendix A.¹⁰² An EOM-CC wave function analysis was performed using the libwa module^78,79 of Q-Chem.

CASSCF and MRCI calculations in Columbus report both h_IJ and d_IJ, but within MRCI, h_IJ depends on the choice of orbital resolution, while the total coupling d_IJ does not.²⁴ In order to avoid this dependence, we use d_IJ and extract h_IJ from it. To do so, we project out the part violating translational invariance and then multiply the so corrected d_IJ by the respective energy gap,

where α = x, y, z and index n runs over all atoms. Although the NAC values corrected by ETF and by projecting out the translationally non-invariant part are not identical, they are reasonably close.³¹ 

The Cartesian geometries and the values of NAC vectors are given in the supplementary material. The details of EOM-CC and MR calculations are given below.

Cis-1,3,5-hexatriene, 1,3-cyclohexadiene, and uracil cations. The cc-pVDZ basis was used in all calculations. In EOM-IP-CCSD calculations, core electrons were frozen. The complete details of the MR calculations can be found in Ref. 80. For the uracil and hexatriene cations, the state-averaged (SA) CASSCF method was used. For hexatriene, the active space included six π-orbitals and the corresponding five electrons [denoted (5,6)] and state-averaging was performed over three states. For uracil, the active space was (13,10), comprising all eight π-orbitals and two lone pairs on the oxygen atoms; the averaging was performed over four states. For cyclohexadiene, a restricted active space SCF (RASSCF) was used in order to obtain the same ordering of states as in EOM-IP-CCSD. The details of validation of these expansions and correct state ordering are given in Ref. 80.

Excited states of uracil. The cc-pVDZ basis was used in all calculations. In the EOM-EE-CCSD calculations, core electrons were frozen. An active space of 12 electrons in 9 orbitals (12,9) was used in all cases, comprising all π orbitals and one lone pair on oxygen, leading to 2520 reference CSFs (configuration state functions). The CASSCF calculations were performed using three-states averaging (SA3-CASSCF) in most cases. For additional tests along mode 25, SA4-CASSCF was performed. Single excitations from the active space were used in the MRCI expansion (denoted MR-CIS) using SA-CASSCF orbitals. In an additional MRCI calculation denoted MRCI1 (for mode 25), a larger expansion that included also all single excitations from all σ orbitals was used. The chosen active space and the MRCI expansions have been previously benchmarked and used in many previous studies of uracil.^81–83 

Na$(NH3)n$, n = 1…3. A mixed double-zeta basis set (aug-cc-pVDZ on heavy atoms and cc-pVDZ on hydrogens) was used for sodium-doped ammonia clusters. All electrons were active in the EOM-EA-CCSD calculations. In MRCI calculations, orbitals obtained from CASSCF averaged over four states (SA4-CASSCF) were used. The active space in CASSCF and MRCI consisted of 1-electron-in-4-orbitals (3s and 3p of Na). The MRCI expansion included all single and double excitations from all orbitals except the core, giving rise to 893 703, 5 304 376, and 18 072 665 CSFs for the complexes with 1, 2, and 3 NH₃, respectively.

We validated our implementation by comparing EOM-IP-CCSD and EOM-EE-CCSD NAC forces against a modified Cfour implementation of quasidiabatic couplings⁵⁵ (see  Appendix A) for the selected examples. In all cases, the NAC forces computed by the two programs agreed within 6 decimal points. EOM-EA-CCSD NAC forces were validated against EOM-EE using the diffuse orbital trick.⁸⁴ The Cartesian geometries, details of the benchmark calculations, and the computed values of the NAC forces are given in the supplementary material. We note that since the values of the NAC force depend on the molecular orientation, for proper comparison between the methods, the same Cartesian geometries must be used.

For reliable treatment of non-adiabatic dynamics, an electronic structure method should be capable of describing accurately the shape and relative energies of the PESs to ensure the correct description of the MECP as well as the magnitude and the direction of the NAC force. Many previous studies focused on the MECPs.^{34,35,60–62,85–88} In this paper, we focus on the NAC force vector and compare its magnitude and direction between h computed by EOM-CC and multireference methods, e.g., CASSCF, RASSCF, and MRCI. To put the differences between these two high-level methods into a context, for selected examples, we also compare EOM-CC with CIS.

We consider the norm of the difference and the angle between the NAC forces computed by different methods,

Here A and B denote two different methods. For selected examples, we also report the following quantity (computed in the Franck-Condon region) related to the rate of non-adiabatic transition:⁸⁰ 

This is the dot product between the NAC force and the nuclear gradient on the upper (i.e., initially excited) PES. The correlation between r_IJ and the rate on non-adiabatic relaxation, suggested by the last term in Eq. (4), has been illustrated by full-dimensional simulations of non-adiabatic dynamics in several molecules using the surface-hopping approach.⁸⁰ 

A recent study has investigated ultrafast non-adiabatic dynamics in three radical cations: cis-1,3,5-hexatriene, 1,3-cyclohexadiene, and uracil cations.⁸⁰ We use these three cations as a benchmark for comparing EOM-IP NAC force against MR values. In all three cases, we compute electronic properties at the equilibrium geometry of the neutral ground state. The dynamical simulations⁸⁰ have shown that non-adiabatic transitions in these systems occur on a femtosecond time scale and mostly in the Franck-Condon region and that the computed rate of non-adiabatic relaxation correlates well with r_IJ of Eq. (36).

Tables I–III summarize the results for uracil, hexatriene, and cyclohexadiene cations. The state energies, characters, and relevant molecular orbitals are given in the supplementary material (Figs. S1–S3 and Tables S1–S3). As discussed in Ref. 80, for these systems, EOM-IP-CCSD energies and state characters agree well with the multi-reference values. We begin by comparing NAC forces computed by EOM-IP-CCSD and CASSCF/RASSCF. We note that since d_IJ depends on the energy gap between the states, this quantity is more sensitive to the variations of the state energies computed by different methods than the NAC force. As one can see from Tables I–III, the magnitudes of h, as quantified by its norm, are in good agreement between EOM-IP and CASSCF. For most states, the norms of h are within 10% from each other. The only significant difference is observed for the NAC force between states D₀-D₂ and D₁-D₃ in uracil, where the two values differ by ∼30%. Considering Δ and cos θ, Eqs. (34) and (35) afford a more detailed comparison. While the norms of Δ vary from 0.01 to 0.03, the values of cos θ are consistently close to 1 (the smallest value is 0.91) confirming that the NAC forces computed by the two methods are nearly parallel.

The differences between the values of r, Eq. (36), computed by EOM-IP-CCSD and MR methods are somewhat larger than the differences in ||h|| because these quantities also depend on the energy gaps and gradients. Tables I and III show reasonable correlation, that is, MRCI and EOM-CCSD agree in which cases r are small and in which cases r is large. In the latter set of cases, the discrepancies vary from perfect agreement to a factor of 10. In the case of hexatriene (Table II), we observe good agreement for r involving states D₀-D₂ and huge discrepancies for all cases involving D₃. We traced this problem to different state character: apparently, state 3 in CASSCF calculations has different character compared to EOM-CCSD. CASSCF calculations with 5 states revealed that the state which has the same character as D₃ in EOM-CCSD appears as state number 4. The values for CASSCF state D₄ agree very well with EOM-CC values for D₃ (see the supplementary material, Table S4). This case highlights the utility of using several metrics: the changes in the state’s character were revealed by r more prominently than by slightly increased difference in Δ.

As a benchmark system for EOM-EE-CCSD NACs, we consider ¹ππ* and ¹nπ* states of uracil. Figure 3 shows natural transition orbitals (NTOs) and the respective singular values for the S₀-S₁, S₀-S₂, and S₁-S₂ transitions. The S₁-S₂ transition is well described by one NTO pair (participation ratio is 1.02) and has substantial one-electron character (||γ|| = 0.6). The radiationless relaxation in this system was investigated by several groups with a variety of approaches,^81,89–94 most recently by Zhang and Herbert⁹⁵ using SF-TDDFT.³⁸ This study identified two major pathways for excited-state deactivation: direct ¹ππ* → S₀ and indirect ¹ππ* → nπ* → S₀ relaxation. The authors characterized two relevant MECPs: one between ¹ππ* and ¹nπ* and one between ¹ππ* and S₀. Here we denote the two MECPs as CI(S₂S₁) and CI(S₂S₀). In addition to MECPs, Zhang and Herbert reported the optimized excited-state structure of the ¹ππ* state: MIN(S₁). We use these structures to investigate the variations in NAC forces outside the Franck-Condon region, along geometric distortions most relevant to the radiationless relaxation process. Specifically, we consider the displacements along normal modes connecting the Franck-Condon structure (S₀ equilibrium geometry) with the CI(S₂S₁), CI(S₂S₀), and MIN(S₁) structures. We use ground-state normal modes for this calculation, as often done in parallel-mode double-harmonic calculations of the Franck-Condon factors. The ground-state structure and normal modes are computed by ωB97X-D/cc-pVDZ. We identified normal modes that have the largest displacements between the MIN(S₀) and CI(S₂S₁), CI(S₂S₀), and MIN(S₁) structures from Ref. 95. Table IV lists4t5b the normal modes selected for these calculations; the displacements corresponding to these normal modes are visualized in Fig. S4 in the supplementary material. Table S5 in the supplementary material lists the displacements along all normal modes.

Figure 4 shows potential energy profiles along the selected normal modes. The PES computed by different methods looks qualitatively similar. Along modes 3 and 23, the energies of all three states are nearly flat. Along modes 10 and 25, all energies go up steeply. The energy gap between S₁ and S₂ does not change much; however, the gap between S₁ and S₀ shrinks. Along mode 10, there is a kink in the CASSCF and MRCI curves at points 9 and 10 due to an abrupt change in the state character. Along mode 25, both the variations in PESs and the differences in the PES shapes computed by different methods appear most prominent. This is not surprising because from the 4 modes considered here, mode 25 is the stiffest and the displacement along this mode is large.

Figure 5 shows the absolute values of the NAC force between S₁ and S₂ states along the selected modes for CASSCF, MRCI, CIS, and EOM-CCSD wave functions. The raw data are collected in Tables S6–S9 in the supplementary material. Tables S10–S13 in the supplementary material show differences (Δ and cos θ) between the NAC forces computed by different methods against MRCI.

In the Franck-Condon region (point 0 in Tables S6–S13 of the supplementary material and in Fig. 5 and Fig. S5 of the supplementary material), the magnitude of the difference between EOM and MRCI vectors is Δ = 0.012, which is consistent with the values of the respective ||h||, 0.0217 (EOM) and 0.0078 (MRCI). CASSCF shows the smallest difference (0.004). By looking at cos θ, we see that it is not only the magnitude of h but also the direction of the vector which is different: the value of cos θ^EOM is 0.75, whereas for CASSCF and (surprisingly) CIS, these values are larger (0.87 and 0.83, respectively).

As one can see from Fig. 5, the magnitude of ||h|| does not vary much along modes 3 and 23 (except for an increasing difference between CASSCF and MRCI at points 8-10 along mode 3). To better illustrate the trends, Fig. S5 in the supplementary material shows the values of ||h|| normalized to their respective Franck-Condon values. Along mode 23, all methods show a small decline (less than 5%). Along mode 3, EOM shows a slight increase (1%), whereas MRCI and CIS values drop by about 10%. Table S12 of the supplementary material confirms that along mode 23, the differences between the EOM-CC NAC force vectors and reference MRCI values do not change much, for example, Δ^EOM is constant and cos θ^EOM varies between 0.745 and 0.743. Mode 3 shows more complex behavior: at large displacements, the alignment between EOM and MRCI vectors decreases, dropping below 0.5 at points 8-10 (Table S10 of the supplementary material). Interestingly, Δ^EOM shows only a relatively moderate increase at these displacements. We also note that cos θ^CAS decreases consistently with the trends in the respective ||h||. These increasing differences between MRCI and other methods at large displacements are suggestive of changing the wave function composition. For example, at large displacements, there is a noticeable difference even between CASSCF and MRCI. The reason is that at the CASSCF level, states S₂ and S₃ are very close and interact strongly, changing the wave function and, consequently, the coupling.

For mode 10, all ||h|| slightly increase and the trends among all methods are consistent (except points 9 and 10, where CASSCF and MRCI potential energy curves show kinks). Table S11 of the supplementary material confirms this—the changes in Δ and cos θ are small (we note that cos θ^EOM increases along this mode, reaching 0.94 at point 7). This is a case where the states change character leading to sudden changes in the coupling. The left panel of Fig. S6 in the supplementary material shows the four electronic states of uracil at the SA4-CASSCF level. At point 9, states S₂ and S₃ approach each other and possibly cross. NAC highlights this crossing. As can be seen in the right panel at point 9, the values of ||h|| suddenly switch for both S₂ and S₃ states. This leads to a sudden increase in the S₁-S₂ coupling that we see in Fig. 5.

In contrast to modes 3, 10, and 23, mode 25 shows very large variations between the NAC forces computed by different methods. Regardless of the metric used to quantify the coupling vector (norms, normalized norms, Δ, and cos θ), we clearly see that the discrepancies between the methods along this coordinate are large. Apparently, the changes in the wave functions’ characters are very sensitive to the level of correlation treatment. To further investigate this issue, we performed additional MRCI calculation (denoted MRCI1, see Sec. III) using a larger expansion. The results of the two MRCI calculations are shown in Fig. S7 in the supplementary material. As one can see, ||h|| from the two calculations are quite different: in CASSCF and MRCI (which uses a rather compact expansion), the magnitude of ||h|| first increases by 40% (reaching maximum at point 2) and then slightly decreases, whereas in MRCI1, the trend is reversed: ||h|| first drops by 25% (reaching minimum at point 1), and then increases by 40%. This illustrates the sensitivity of NAC along this mode to dynamical correlation. To understand the changes in states’ characters, we analyzed EOM-CCSD wave functions along mode 25 using NTOs of the S₂-S₁ transition. We found that the changes in the NTOs are subtle, and the transition retains its nπ* character at all displacements. Important quantities, such as participation ratio and ||γ||, do not change much. A close inspection of NTOs reveals that the particle NTO for the S₂-S₁ transition does not change much, retaining its out-of-plane n(O₂) character (see Fig. 3). However, the character of the hole orbital evolves (this is shown in Fig. S8 in the supplementary material) from delocalized π orbital to localized out-of-plane n(O₂). Figure S8 of the supplementary material also shows a continuous decrease in electron density in n(O₁) and an increase in n(O₂). The sharp increase in NAC at point 3 occurs when the amplitude of the hole NTO on O₁ disappears. To further understand these changes, we consider natural orbitals of the S₁ and S₂ states. Figure S9 in the supplementary material compares frontier natural orbitals (NOs) computed for the EOM-CC and MRCI wave functions. As one can see, most of the changes occur in the S₂ state. We note that these changes in S₂ (the flow of density from n(O₁) to n(O₂)) have been observed before in a computational study of radiationless relaxation of uracil and related molecules.⁹⁴ The NOs computed by MRCI and EOM-CC appear to be very similar. The only noticeable difference is that EOM-CC wave function at point 4 shows more open-shell character than MRCI (natural occupations of 1.3 and 0.7 versus 1.5 and 0.6). We conclude that the magnitude of NAC is very sensitive to these small variations in the wave function character, which, in turn, are very sensitive to the correlation treatment. Thus, even when qualitatively wave functions are similar, one can observe quantitative discrepancies among different methods.

Figure 6 shows the MRCI NAC force at the Franck-Condon geometry (point 0) and at displaced geometries along modes 23 and 25. As expected from the symmetries and orbital characters (see Fig. 3) of S₁ (A″, nπ*) and S₂ (A′, ππ*), the NAC vector has a″ symmetry (out-of-plane). Thus, only out-of-plane motions are expected to be effective in promoting non-adiabatic transition between these states. Among the four modes considered here, only mode 10 corresponds to out-of-plane distortion. In terms of Eq. (36), only the out-of-plane component of the gradient will contribute to the rate. The NAC force vector along mode 23 looks qualitatively similar to that at the Franck-Condon structure, whereas displacement along mode 25 results in noticeable changes.

We use sodium-doped clusters as a model system for benchmarking EOM-EA. The low-lying electronic states in these clusters, recently studied theoretically⁹⁶ and experimentally,^97–99 can be described as s- and p-like states of surface-bound electrons, stabilized by a positively charged $Naδ+(NH3)n$ core. Dyson orbitals of the four lowest states are shown in Fig. S10 in the supplementary material (and in Fig. 1 in Ref. 96). The target doublet states are best described by EOM-EA using the closed-shell cation as a reference. Non-adiabatic relaxation of p-like states has been studied in related systems, Ba(Ar)_n clusters.¹⁰⁰ 

Table V shows energy gaps and the NACs between the three lowest electronic states of the sodium-doped clusters, Na$(NH3)n$⁠, with n = 1-3. The NAC forces for Na(NH₃) and Na$(NH3)2$ are shown in Fig. 7. As one can see, the NAC vectors are dominated by the relative Na-ammonia coordinate. Interestingly, while differences between the methods in the magnitude of ||h|| are consistent with Δ, the value of cos θ appears rather low. The reason is that in these three clusters, the three p-like states are nearly degenerate, which means that small differences in the methods can mix them. Since the mixing of the p-orbitals only affects their direction, such mixing only affects the direction of the NAC forces but not their magnitude. In support of this explanation, we note that the largest cos θ values correspond to the NAC forces between the s-state with the highest p-state (states 1-4) in NaNH₃ and Na$(NH3)3$ in which the gap between state 4 and states 2-3 is the largest.

We presented the theory and implementation of the NAC forces within the EOM-CCSD framework. The main focus of this paper is on the comparison between the EOM-CCSD and MRCI NAC forces. Using selected examples, we analyzed NACs computed using EOM-IP-CCSD, EOM-EE-CCSD, and EOM-EA-CCSD wave functions. In contrast to previous studies, which focused primarily on the structures of conical intersections, here we analyzed the direction of the NAC vector. We note that the direction of NAC vector, not only its magnitude, is very important for the calculations of the rate of non-adiabatic transitions. Our examples illustrate the utility of using different metrics of comparison. The numeric results confirm that the accuracy of NAC forces depends critically on the quality of the underlying wave functions, which are often very sensitive to correlation treatment. Within their domain of applicability, EOM-CC methods provide a viable alternative to MR approaches.

See supplementary material for computational details, additional validation and benchmark calculations, electronically exited states of uracil and sodium-doped ammonia clusters, and relevant Cartesian coordinates.

We are grateful to Professor John Stanton for his help in validating the implementation and for the enlightening analysis of the formalism, partially summarized in  Appendix A. We also thank Dr. Evgeny Epifanovsky and Dr. Xintian Feng from Q-Chem for their help with the implementation. Finally, we thank Dr. Thomas Jagau for his feedback about the manuscript and stimulating discussions.

This work is supported by the Department of Energy through the DE-FG02-05ER15685 grant and by the U.S. Air Force of Scientific Research (AFOSR) under Contract No. FA9550-16-1-0051 (A.I.K.). S.M. acknowledges support from NSF Grant No. CHE-1465138. S.F. was supported by the DAAD fellowship from DFG during her stay at USC.

To derive expressions for the EOM-CC NACs, we begin by differentiating an interstate matrix element of $H¯$⁠, which, obviously, is zero,

By defining $H¯x≡e−TH′eT$ and using the usual notation $f′≡∂f∂x$⁠, we obtain

where p denotes the EOM-CCSD manifold of Slater determinants: p = 0 + S + D (and the complimentary space of higher excitations is denoted by q = T + Q + ⋯).

By using the amplitude-response operator Ξ,

we obtain

Because ⟨0|L_IR_J|0⟩ is a constant,

Thus, we obtain

We now define EOM-CC derivative coupling as

In this definition, we neglect orbital response terms except for those that arise from the differentiation of $H¯$⁠.

By combining Eqs. (A6) and (A7), we obtain the following expression for $dIJC$⁠:

The numerator in this equation is the EOM NAC force from Eq. (15). It is equal to quasidiabatic couplings from Ref. 55 minus the following term:

As our benchmarks illustrate, the magnitude of this term is rather small. Thus, the EOM-CC $dIJC$ is nearly identical to the quasidiabatic couplings λ_IJ from Ref. 55.

The orbital-response part of the full derivative coupling is given by the following expression:¹⁴ 

where $D̃pqIJ$ is the antisymmetric one-particle transition density matrix,

and $⟨ϕp|ϕqx⟩$ is¹⁴ 

where the first term is the derivative of the overlap matrix and the second part is the vector of coupled-perturbed Hartree-Fock coefficients describing the response of the MO coefficients to the moving of the ξ atomic basis functions. The calculation of this term has been described by Gauss et al.⁸ The $dIJϕ$ term violates the translational invariance of the NAC and is omitted in our implementation. It is also omitted in the implementation of quasidiabatic couplings of Stanton and co-workers.⁵⁵ 

The normalization factors for left EOM states are defined as^14,101

This equation can be tackled by inserting the resolution of the identity, $∑xΦxΦx$ (x denotes the excitation level). For EOM-EE-CCSD, only the following three terms survive, giving rise to the following programmable expression:

The $NxR$ is approximated as

This is exact in the case of FCI.¹⁴ 

We note that in EOM-SF version, the first term in the above expression is zero. For EOM-IP and EOM-EA, the expressions for the left state norms are

and