Equation-of-motion coupled-cluster singles and doubles (EOM-CCSD) is a reliable and popular approach to the determination of electronic excitation energies. Recently, we have developed a rank-reduced CCSD (RR-CCSD) method that allows the ground-state coupled-cluster energy to be determined with low-rank cluster amplitudes. Here, we extend this approach to excited-state energies through a RR-EOM-CCSD method. We start from the EOM-CCSD energy functional and insert low-rank approximations to the doubles amplitudes. The result is an approximate EOM-CCSD method with only a quadratic number (in the molecular size) of free parameters in the wavefunction. Importantly, our formulation of RR-EOM-CCSD preserves the size intensivity of the excitation energy and size extensivity of the total energy. Numerical tests of the method suggest that accuracy on the order of 0.05–0.01 eV in the excitation energy is possible with 1% or less of the original number of wavefunction coefficients; accuracy of better than 0.01 eV can be achieved with about 4% or less of the free parameters. The amount of compression at a given accuracy level is expected to increase with the size of the molecule. The RR-EOM-CCSD method is a new path toward the efficient determination of accurate electronic excitation energies.

## I. INTRODUCTION

The equation-of-motion^{1–4} coupled-cluster (EOM-CC) approach is among the most robust electronic structure methods for the determination of molecular excited states. When truncated at the level of single and double excitations—equation-of-motion coupled-cluster singles and doubles (EOM-CCSD)—this method provides accurate, black-box excitation energies at comparable computational cost to the underlying coupled-cluster singles and doubles (CCSD) treatment of the ground state.^{5–7} Unfortunately, this computational cost is quite high and scales with the sixth-power of system size, O(*N*^{6}), which severely limits the size of systems that can be studied with EOM-CCSD. Furthermore, the number of parameters in the EOM-CCSD wavefunction scales with the fourth-power of system size, O(*N*^{4}), leading to a storage bottleneck associated with the parameterization of the double excitations. In the present work, we develop an approach for a rank-reduced EOM-CCSD (RR-EOM-CCSD) method that offers the potential for reductions in the scaling of both storage and computational requirements with respect to system size.

Although there have been attempts to reduce the computational expense of the EOM-CCSD method, progress in this area has primarily been limited to methods for reducing multiplicative prefactors of the computational cost as opposed to methods that attempt to reduce the asymptotic scaling. Typically, approaches of this sort do not significantly increase the maximum size of systems that can be studied, but they do make computations on smaller molecules more routine. Some of the early attempts to improve the performance of EOM-CCSD were due to Koch and co-workers via an integral direct formulation.^{8–10} This can avoid bottlenecks associated with the storage of electron repulsion integrals (ERIs), and, by exploiting sparsity in the atomic orbital basis, the complexity of certain terms may be reduced for sufficiently large molecules. More recently, Krylov and co-workers have applied Cholesky^{11–15} and resolution-of-the-identity decompositions (RIs)^{16–20} of the electron repulsion integrals (ERIs) in the context of EOM-CCSD to reduce the storage bottleneck associated with certain classes of integrals.^{21} They report about a factor of two reduction in computational time for EOM-CCSD using these low-rank factorizations of the ERIs (with additional efficiency gains possible through the truncation of the virtual space).^{21} Recently, an improved algorithm for performing Cholesky decompositions of the ERIs has been developed that may lead to new Cholesky-based coupled-cluster algorithms with improved performance.^{15} For many years, it has been known that pseudospectral integral approximations^{22–26} allow for reduced-scaling factorizations of certain classes of diagrams that are common to numerous many-body methods.^{27–31} We applied these ideas to ground state CCSD wavefunctions through the introduction of tensor hypercontraction (THC)^{32–34} approximations to reduce the expense associated with the rate-limiting contraction of the cluster amplitudes and four-virtual-index ERIs (sometimes known as the particle-particle-ladder diagram).^{35} Subsequently, Neese and co-workers have applied this idea directly to EOM-CCSD.^{36,37} Since this reduces the asymptotic scaling of certain terms (and not just the prefactor), pseudospectral-based approaches are far more profitable than RI-like integral factorizations in the context of EOM-CCSD. Neese and co-workers report speedups of a factor of 10 due to this approximation.^{36}

For ground-state coupled-cluster wavefunctions, local methods (particularly those based on pair natural orbitals, PNOs) have become extremely popular in recent years.^{38–42} These methods allow coupled-cluster energies to be determined for systems with hundreds of atoms on a fairly routine basis. However, extensions of local methods to electronic excited states have proven to be far more challenging.^{43} Since excited-state wavefunctions are formulated in terms of nonlocal changes to the ground state wavefunction, the sparsity patterns that exist on the ground state may not be relevant for excited states. The need to describe multiple excited states simultaneously further complicates matters as the sparsity pattern may differ for each excited-state. Many recent efforts have focused on approaches that leverage a local CCSD wavefunction (or wavefunctions) to construct more compact parameterizations of the excited state wavefunction.^{44–47} In particular, pair natural orbitals have been applied to similarity-transformed equation-of-motion coupled-cluster singles and doubles (STEOM-CCSD).^{48} This method allows a reasonable treatment of excitation energies to be obtained from the singles-singles block of the STEOM Hamiltonian. As a result, STEOM-CCSD is more compatible with PNOs since the local approximations do not need to be extended to the excited state wavefunction parameters. There has been some preliminary work on the formulation of a local EOM-CCSD method based on PNOs, but it is unclear how effective PNO-based approaches can be for excited states; this remains an active area of research.^{43,49}

It has been understood for many years that the *t*-amplitudes defining the CCSD wavefunction exhibit rank-sparsity. In particular, the number of doubles amplitudes is quartic in the molecular size but can be accurately described with a quadratic number of parameters—this can be demonstrated by singular value decomposition (SVD).^{33,50–52} For large enough molecules, the number of parameters needed may grow only linearly with molecular size. There have been numerous attempts to develop approximate CCSD methodologies that exploit this rank-sparsity.^{33,50–55} The difficulty in developing such a method is that the *t*-amplitudes are defined through a nonlinear system of equations that is most naturally expressed in terms of full-rank tensors and it is not immediately obvious how to solve these equations in terms of low-rank quantities without involving full-rank intermediates. While it has been demonstrated previously^{33,50,52,55,56} that it is possible to solve for the CCSD (or CCSD-like) wavefunction in a compressed representation, none of these earlier approaches have been extended to electronic excited states. Recently, we have developed a rank-reduced CCSD method by starting from an energy Lagrangian containing only low-rank quantities.^{56} A particularly important feature of our Lagrangian formulation is that it provides a path forward from ground state energies to molecular excited states and other response properties. In the present work, we extend our rank-reduced CCSD method to excited states via EOM-CCSD. We demonstrate for a series of representative examples that chemical accuracy can be achieved with a compressed EOM-CCSD wavefunction containing only a quadratic number of free parameters, without ever accessing or enumerating the full quartic set of doubles amplitudes.

## II. THEORY

### A. Notation

We will adhere to the following convention for the labeling of different types of indices.

: Indexes occupied molecular orbitals.*i,j,k,l*: Indexes virtual molecular orbitals.*a,b,c,d*: Indexes any molecular orbital.*p,q,r,s*: Indexes density fitting (or Cholesky decomposition) auxiliary functions.*A,B,C,D*: Indexes the rank-reduced coupled cluster factors.*V,W,X,Y*: Indexes the quadrature points relevant to Laplace transformation of the orbital energy denominator.*w*

A generalized Einstein summation convention is used where summation over repeated indices is implied.

### B. Canonical EOM-CCSD

We will begin with the coupled-cluster (CC) form of the Schrödinger equation,

Here, $T^$ is the cluster excitation operator,

After left-multiplication by $e\u2212T^$,

equations defining the ground state energy and t-amplitudes can be obtained via left-projection into a basis of Slater determinants,

This is the conventional definition of the ground-state CC energy and wavefunction. There are several techniques for extending CC to molecular excited states; the two most common are the equation-of-motion coupled-cluster method (EOM-CC)^{1–7} and coupled-cluster linear response theory (CCLR).^{57,58} These two methods produce *identical* excitation energies and only begin to differ in their definition of one-electron properties. Both approaches have advantages and disadvantages. CCLR provides a better treatment of one-electron properties—for example, size intensive transition dipole moments—although additional computational effort is often required.^{59} On the other hand, it is more natural to extend the EOM-CC formalism to approaches that do not conserve the electron number or the spin of the reference state.^{60} In the present work, we focus on the EOM-CC formalism.

In the traditional formulation of EOM-CC, the energies of the ground and excited states are the eigenvalues of the similarity transformed Hamiltonian viz.,

Here, the $R^n$ operator is a CI-like linear excitation operator (for the *n*th electronic state),

Equivalently, EOM-CC can be formulated as a left-hand eigenvalue problem,

where $L^n$ is a linear de-excitation operator corresponding to the *n*-th state,

The EOM-CC energy functional can be written as^{61,62}

and the right- and left-hand eigenvectors form a biorthogonal basis, viz.,

The EOM-CC energy is variational with respect to the *r*- and *l*-amplitudes,

This usual formulation of the EOM-CC wavefunction is the starting point of our present work.

In EOM-CCSD, which is the focus of the present work, only single and double excitations are included.^{7} Furthermore, amplitudes of the $T^$ operator are defined by the left-projection shown in Eq. (5). Other truncations are certainly possible (as well as other definitions of the *t*-amplitudes). However, several important and desirable properties stem from the consistency of this definition. As a direct consequence of Eq. (5), the matrix elements coupling the ground and excited states in the similarity transformed Hamiltonian are zeroed. This decoupling of the ground and excited states results in an analytic form for the right-hand EOM-CC ground-state wavefunction,

Equation (5) also forms the basis of the relationship between EOM-CC and CCLR because the EOM-CC ground state is a stationary reference on which CCLR is formulated.

Another property resulting from a consistent definition of the linear and cluster excitation operators is the size intensivity of the EOM-CC excitation energy. This can be demonstrated, as shown by Bartlett,^{63} by left-multiplying the (right-handed) ground-state EOM-CC eigenvalue problem [Eq. (3)] by the $R^n$ operator and subtracting it from the (right-handed) excited-state EOM-CC eigenvalue problem [Eq. (6)],

The size intensivity of the excitation energies is preserved because they can be written entirely in terms of commutators between the linear excitation operator and the similarity transformed Hamiltonian (i.e., only linked diagrams contribute to the excitation energy). Terms containing the ground-state correlation energy are formally removed from the definition of the *r*-amplitudes. It is essential to realize that the proper scaling of energies in EOM-CC depends entirely on the definition of the *t*-amplitudes; it must be the case that Eq. (5) defines the EOM-CC ground-state wavefunction.

### C. Rank-reduced formulation of EOM-CCSD

One of the principal limitations of the EOM-CCSD method is its computational expense. The computational cost of the method scales with the sixth-power of system size, O(*N*^{6}), and the storage requirements (for the doubles amplitudes, $rijab$) scale with the fourth-power of system size, O(*N*^{4}). In the present work, we will reformulate the EOM-CCSD method to exploit a low-rank factorization of the doubles amplitudes of the form

Here, the $UiaV$ tensors are low-rank projectors from the full *ov* space to a compressed auxiliary space. If the projector can be chosen such that the size of the auxiliary space (the range of *V*) grows linearly with system size (obviously, the full *ov* space grows quadratically), the number of free parameters in the EOM-CCSD wavefunction can be reduced from O(*N*^{4}) to O(*N*^{2}) by replacing the $rijab$ tensor with the *R*^{VW} matrix. At present, we will work exclusively with unitary projectors, viz.,

This restriction is a practical consideration; extending the formalism to nonunitary projectors is straightforward.

Our goal is to develop a rank-reduced EOM-CCSD method (RR-EOM-CCSD) that directly solves for the compressed EOM-CCSD wavefunction. To that end, we introduce rank-reduced ansätze for the $R^n$ and $L^n$ operators,

It should be noted that the choice to use a single projector for both the *r*- and *l*-amplitudes is not required and is only done here for simplicity; the use of separate projectors could potentially yield somewhat greater compression; however, this possibility will not be explored in the present work. Now, we may write the EOM-CC energy functional in terms of the rank-reduced quantities, remembering that this energy functional is a Lagrangian depending explicitly on the *r*- and *l*-amplitudes and depending parametrically on the *t*-amplitudes and low-rank projector,

To generate equations for the rank-reduced *r*-amplitudes, we may simply set the partial derivative of the energy (Lagrangian) with respect to the *l*-amplitudes to zero,

This results in a rank-reduced analog of the usual EOM-CCSD equations.

The working equations defining the canonical EOM-CCSD *r*-amplitudes are typically formulated as a matrix eigenvalue problem. Due to the size of the matrix (approximately *o*^{2}*v*^{2} × *o*^{2}*v*^{2}), the eigenvectors are found using a nonsymmetric variant of Davidson diagonalization.^{64} For EOM-CCSD, the equations that must be satisfied can be written in the following form:

Here, we are solving for the total EOM-CCSD energy, rather than the excitation energy. When we apply rank-reduction to the *r*- and *l*-amplitudes, the partial derivatives of the rank-reduced EOM-CC Lagrangian lead to a new set of rank-reduced EOM-CC equations that are projected into a compressed space. For a compact presentation of these new equations, we introduce the compressed doubles space,

In the present work, we assume the underlying functions to be singly and doubly substituted Slater determinants; extension to a basis of spin-adapted configuration state functions is straightforward. Using the compressed doubles space, we may write the RR-EOM-CC eigenvalue problem in terms of the compressed wavefunction coefficients,

Now, it is obvious that the RR-EOM-CC energies and wavefunctions are the eigenvalues and eigenvectors of the compressed similarity transformed Hamiltonian. The Lagrangian formulation allows the compressed wavefunction to be determined, directly, from a standard matrix eigenvalue problem.

An alternative formulation of RR-EOM-CC could allow the Lagrangian (energy functional) to depend explicitly on the low-rank projector. In that case, one would solve variationally for both the compressed doubles amplitudes, *R*^{VW}, and the corresponding projector $UiaV$. This might lead to a more compact representation of the wavefunction since the full variational flexibility of the rank-reduced ansatz could be utilized. However, simultaneous optimization of both *R*^{VW} and $UiaV$ cannot be cast as a simple eigenvalue problem. We do not pursue this possibility in the present work. Similarly, we do not yet consider combining our previous work on RR-CCSD^{56} to obtain the ground state for RR-EOM-CCSD. In that context, it might be possible to define a Jacobian in terms of the rank-reduced ground-state wavefunction to obtain excitation energies within a response formalism.^{58} The use of RR-CCSD to generate the RR-EOM-CCSD ground state will be considered in a future work.

It is clear from this formulation that the ground state, right-hand wavefunction has *r*_{0} = 1, $ria$ = 0 and *R*^{VW} = 0. This follows from the fact that we have not modified the equations defining the *t*-amplitudes and they are taken to be solutions of the usual CCSD equations [Eq. (5)]. As a result, the matrix elements coupling the reference determinant to the singly and doubly excited subspaces of the compressed similarity transformed Hamiltonian, $\Phi iae\u2212T^\u0124eT^\Phi 0$ and $\Phi VWe\u2212T^\u0124eT^\Phi 0$, are trivially zero. This property guarantees that the ground state RR-EOM-CCSD energy is the usual CCSD energy. Furthermore, the structure of the compressed similarity transformed Hamiltonian allows for the removal of the ground state energy and the excitation energy to be solved directly [following from the matrix equivalent of Eq. (16)]. Finally, this proves that the size intensivity of EOM-CCSD excitation energies has been exactly preserved in our rank-reduced formulation.

One perspective on RR-EOM-CCSD is that the similarity transformed Hamiltonian itself is compressed. Another equivalent perspective is that the usual EOM-CCSD equations are solved with a compressed number of free parameters. In that case, the variational energy functional, Eq. (10), is minimized using parameters within the auxiliary space rather than parameters in the full *ov* space. In that way, RR-EOM-CCSD can be thought of as a constrained optimization of the EOM-CCSD wavefunction, where the projector serves to constrain the search space. Perhaps the most important result stemming from this alternative perspective is that the RR-EOM-CCSD energies provide an upper bound to the exact EOM-CCSD energy, which is simply to say that the RR-EOM-CCSD energies are variational. From the practical perspective, this property could potentially be exploited to develop extrapolation strategies to obtain more accurate RR-EOM-CCSD energies by increasing the number of free parameters included in the wavefunction. These ideas are similar to recent work applying model-order-reduction (MOR) to the coupled-cluster Green’s function.^{65} There, compression of the similarity transformed Hamiltonian was accomplished via the application of MOR, and extrapolated results obtained from the reduced-order models were found to be useful in this context.

### D. Low-rank projectors for RR-EOM-CCSD

The accuracy and efficiency of the RR-EOM-CCSD method depends entirely on the choice of the low-rank projectors, $UiaV$. These projectors can be defined as a subset of the eigenvectors of some generating matrix,

The $Gijab$ generator is a fourth-order tensor that is treated as a *ov* × *ov* matrix for the purpose of eigendecomposition. The low-rank projectors, $UiaV$, are the subset of the eigenvectors chosen from the eigenvalues with the largest magnitude. That is, we introduce some threshold, $\epsilon $, and include all eigenvectors with |*g*^{V}| > ϵ in the projector. The optimal choice for the generator, in the sense of the 2-norm, is the matrix containing the exact doubles amplitudes,

Although this choice offers the most efficient compression of the EOM-CCSD wavefunction, it requires the full-rank EOM-CCSD wavefunction to be known prior to the construction of the low-rank approximation. Therefore, we must develop a proxy that captures the essential physics of the $rijab$ amplitudes without requiring knowledge of the exact EOM-CCSD wavefunction.

Practically speaking, choosing a projector requires balancing the accuracy of the approximate $rijab$ amplitudes against the efficiency of constructing the projector. Taking inspiration from approximate coupled-cluster methods such as CC2,^{66} we introduce a simple approximation for the doubles amplitudes,

where *ω* is the excitation energy. This allows us to write a closed-form expression for the generator,

that requires as input only a reasonable approximation for the $ria$ amplitudes and the excitation energy, *ω*. In terms of both accuracy and practicality, the best source of $ria$ amplitudes with which to construct the generator defined in Eq. (33) seems to be RR-EOM-CCSD itself. Instead of taking the projector to be a parameter on which the RR-EOM-CCSD Lagrangian depends [Eq. (21)], we define the projector to be a function of the $ria$ amplitudes. Since the projector now depends on part of the eigenvector, the RR-EOM-CCSD equations become nonlinear. This allows all of the information available during the solution of the RR-EOM-CCSD equations to be used in the construction of the projector defined by the generator in Eq. (33). This approach, particularly the choice of approximate doubles amplitudes, is closely related to the correlated natural transition orbital framework for low-scaling excitation energy calculations (CorNFLEx).^{67} In that approach, the approximate doubles amplitudes given by Eq. (33) are used to generate a corresponding one-particle density matrix and, subsequently, natural transition orbitals that form a reduced orbital space. This method has been applied to improve the efficiency of CC2 excitation energy computations.

We solve the nonlinear eigenvalue problem that results from the use of the RR-EOM-CCSD $ria$ amplitudes in the definition of the projector by fixed-point iteration. We linearize the problem by holding the projector constant and solving the RR-EOM-CCSD eigenvalue problem. Then, the projector is updated to reflect the new $ria$ amplitudes. The process of solving the eigenvalue problem and updating the projector is repeated until the $ria$ amplitudes used to construct the projector are self-consistent with the eigenvectors of the compressed similarity transformed Hamiltonian. At present, we fully solve a series of linearized problems until the change in the eigenvalue falls below a specified convergence threshold. We initialize the procedure by taking singles amplitudes from a configuration interaction singles (CIS)^{68} wavefunction or an average over several CIS states. There are possible improvements to the efficiency of this approach that we have yet to consider. One such possibility is to solve the nonlinear eigenvalue problem directly. A simpler possibility would be to devise an algorithm that uses looser convergence criteria for early iterations and only tightly converges the linearized problems when a good estimate of the projector is available.

One of the primary motivations for the generator presented in Eq. (33) is that its eigenvectors can be found with O(*N*^{4}) effort, provided that a low-rank decomposition of the electron repulsion integrals of the form

is available. This may be obtained via density fitting or Cholesky decomposition approximations.^{11,12,16–18} Defining a new intermediate,

and a Laplace transformation of the shifted orbital energy denominator,^{69,70}

(where λ_{w}/α_{w} indexes one of the N_{w} quadrature weights/points), it is possible to express the generator amplitudes of Eq. (33) in a simple form,

We may further simplify this expression by writing a *ov* × 2*N*_{DF} supermatrix containing the $R\u0303iaA$ and $L\u0303iaA$ quantities (treated as $ov\u2009\xd7\u2009NDF$ matrices),

Here, we use *B* to index the columns of the supermatrix. Now, the generator may be written as

To obtain the projector, we will find the eigenvalues/eigenvectors of the square of the generator (again, treated as an *ov* × *ov* matrix),

Since we will eigendecompose the square of the generator, the sign of the eigenvalues, *g*^{V}, is lost. However, we consider the absolute value of the eigenvalue when constructing the projector, so this is not problematic in the present case.

Because the $Gijab$ matrix is symmetric, we may expand the square of the generator as

We treat the *Z* tensor as an *ov* × 2*N*_{DF}*N*_{w} matrix and express it in terms of its singular value decomposition (SVD),

Inserting this decomposition into Eq. (41),

we can see that through the formation of a new intermediate,

the square of the generator may be expressed in a simple form,

Now, the projector (eigenvectors of $Gijab$) may be obtained via an SVD of the $Z\u0303$ tensor,

The $UiaV$ singular vectors are the eigenvectors of $Gijab$, and the singular values are the absolute values of the $Gijab$ eigenvalues.

### E. Multiple electronic states in RR-EOM-CCSD

In the preceding discussion, we present a formulation of RR-EOM-CCSD that uses a single projector to compress the similarity transformed Hamiltonian. Multiple electronic states may be obtained from this compressed Hamiltonian, and the usual biorthonormality of those states [Eq. (11)] is preserved. However, this presupposes the existence of a single projector that can simultaneously describe multiple electronic states. Although it is certainly possible to construct such a projector, we have found, in practice, that much of the efficiency is lost when multiple states are considered. Essentially, state-specific projectors are so efficient in their description of one particular state that they are unable to describe any other electronic state with similar accuracy. We have found that the size of the auxiliary space would grow nearly linearly with the number of electronic states that are included in the construction of the projector. For some particular applications (such as a chromophore in solution or a protein environment), the larger auxiliary space might not be problematic. In general and for the smaller molecules currently under consideration, however, we must apply state-specific projectors in order to obtain the amount of compression that is desired.

In order to obtain multiple excited states within the RR-EOM-CCSD framework while also using state-specific projectors, we develop an algorithm based on deflation that allows each state to be determined individually. The idea of this approach is to determine one eigenpair at a time and then remove it from the problem,

Here, *r*^{(n)} are the vectors corresponding to the *n*th eigenvalue of the similarity transformed Hamiltonian, $H\xaf$, and *R*^{⊥} is a subspace that is orthogonal to the eigenvectors that have been determined. The *r*^{n+1} vector is obtained as an eigenvector of $H\xaf\u22a5$. This process can be repeated until $H\xaf\u22a5$ vanishes and all of the eigenvalues have been determined; at this point, a unitary matrix has been constructed which transforms the original matrix into an upper triangular form [the dense upper triangle is represented by the asterisk on the right-hand-side of Eq. (48)] with the eigenvalues of the matrix on the diagonal.

This approach is known in general as Wielandt deflation, which, for one vector, may be expressed simply as

Here, $H\xaf1$ is the deflated similarity transformed Hamiltonian, *σ* is a shift parameter, *r*_{1} is the first right-hand eigenvector of the similarity transformed Hamiltonian, and *v* is an arbitrary vector chosen such that *v*^{†}*r*_{1} = 1. One may expect that *v* should be chosen to be the left-hand eigenvector of the similarity transformed Hamiltonian corresponding to *r*_{1}; this choice is known as Hotelling deflation and preserves both the left- and right-hand eigenvectors during the deflation procedure. If only the eigenvalues are of interest, as in the present case, it is useful to choose *v* = *r*_{1}, which preserves the Schur vectors of the similarity transformed Hamiltonian but not the eigenvectors. One can see that the deflation of the similarity transformed Hamiltonian outlined in Eq. (48) is equivalent to the construction of a Schur decomposition one vector at a time. A more detailed description of deflation-based approaches to eigenvalue problems can be found in Ref. 71.

If we obtain multiple RR-EOM-CCSD states using state-specific projectors, each state is an eigenvector of a different compressed similarity transformed Hamiltonian. As a result, the deflation-based approach described above cannot be directly applied. We must make an approximation, which is to project the previously determined RR-EOM-CCSD states into the subspace defined by the new projector and to solve for the next state in the orthogonal subspace. This is an approximation because the old RR-EOM-CCSD states do not correspond to eigenvalues of the new compressed similarity transformed Hamiltonian. As a result, the decoupling of the previously determined eigenvalues from the rest of the matrix shown in Eq. (48) does not occur. Therefore, the RR-EOM-CCSD energies are not eigenvalues of their compressed similarity transformed Hamiltonian. Note that this approximation is removed when the auxiliary space defined by the projector spans the complete *ov* space. Numerical tests will determine whether or not this approach to the solution of the RR-EOM-CCSD energies is reasonable.

### F. Implementation details of RR-EOM-CCSD

In Algorithm 1, we sketch the algorithm used for solving for an arbitrary number of RR-EOM-CCSD excitation energies. The guess $ria$ amplitudes can be obtained from a CIS calculation; if the CIS state resembles the final EOM-CCSD state of interest, this is often sufficient to seed the procedure. Alternatively, it is possible to use an average of several CIS wavefunctions to construct guess $ria$ amplitudes; several RR-EOM-CCSD states can then be determined using a single projector constructed from the state-averaged amplitudes. The $ria$ amplitudes from these RR-EOM-CCSD wavefunctions also serve as a good initial guess for the self-consistent determination of state-specific projectors.

1: | Solve HF equations |

2: | Solve CCSD equations |

3: | for n States do |

4: | Initialize: k=0 |

5: | Construct guess $ria,n$ amplitudes |

6: | Form $Gijab,n$ using $ria,n$: Eq. (33) |

7: | Form $UiaV,(n)$ by diagonalizing $rijab,n$: Eq. (30) |

8: | if n > 0 then |

9: | $R\u22a5=\xce\u2212\u2211n=0,1\u2026n\u22121R^n\Phi 0\Phi 0R^n\u2020$ |

10: | end if |

11: | Solve RR-EOM-CCSD: Eqs. (27)–(29) |

12: | Set: $Ekn=ERR-EOM-CCn$ |

13: | Update $ria,n$ amplitudes |

14: | loop |

15: | Rebuild $Gijab,n$ using $ria,n$: Eq. (33) |

Rebuild $UiaV,(n)$ by diagonalizing $Gijab,n$: Eq. (30) | |

16: | if n > 0 then |

17: | $R\u22a5=\xce\u2212\u2211n=0,1\u2026n\u22121R^n\Phi 0\Phi 0R^n\u2020$ |

18: | end if |

19: | Solve RR-EOM-CCSD: Eqs. (27)–(29) |

20: | Set: $Ekn=ERR-EOM-CCn$ |

21: | Update $ria,n$ amplitudes |

22: | if $Ek+1n\u2212Ekn$ then |

23: | return $ERR-EOM-CCn$ |

24: | end if |

25: | Set: k=k+1 |

26: | end loop |

27: | end for |

1: | Solve HF equations |

2: | Solve CCSD equations |

3: | for n States do |

4: | Initialize: k=0 |

5: | Construct guess $ria,n$ amplitudes |

6: | Form $Gijab,n$ using $ria,n$: Eq. (33) |

7: | Form $UiaV,(n)$ by diagonalizing $rijab,n$: Eq. (30) |

8: | if n > 0 then |

9: | $R\u22a5=\xce\u2212\u2211n=0,1\u2026n\u22121R^n\Phi 0\Phi 0R^n\u2020$ |

10: | end if |

11: | Solve RR-EOM-CCSD: Eqs. (27)–(29) |

12: | Set: $Ekn=ERR-EOM-CCn$ |

13: | Update $ria,n$ amplitudes |

14: | loop |

15: | Rebuild $Gijab,n$ using $ria,n$: Eq. (33) |

Rebuild $UiaV,(n)$ by diagonalizing $Gijab,n$: Eq. (30) | |

16: | if n > 0 then |

17: | $R\u22a5=\xce\u2212\u2211n=0,1\u2026n\u22121R^n\Phi 0\Phi 0R^n\u2020$ |

18: | end if |

19: | Solve RR-EOM-CCSD: Eqs. (27)–(29) |

20: | Set: $Ekn=ERR-EOM-CCn$ |

21: | Update $ria,n$ amplitudes |

22: | if $Ek+1n\u2212Ekn$ then |

23: | return $ERR-EOM-CCn$ |

24: | end if |

25: | Set: k=k+1 |

26: | end loop |

27: | end for |

We solve the RR-EOM-CCSD equations using a nonsymmetric Davidson diagonalization.^{64} The formation of the residual follows from Eqs. (27)–(29). In practice, we solve for the excitation energy directly in terms of $ria$ and *R*^{VW} and determine *r*_{0} ex post facto. Regarding the solution of these equations, the only open question is how the residual should be preconditioned. In EOM-CCSD, the preconditioner is usually taken to be the diagonal Fock matrix contributions to the similarity transformed Hamiltonian. In that case, the Davidson correction vectors are formed as

where $\omega \u0303$ is the current estimate of the eigenvalue (excitation energy) and *ρ* is the residual of the EOM-CCSD eigenvector. In the rank-reduced formulation, however, the preconditioner for the doubles amplitudes is not diagonal,

This can be expanded into a form that makes the system of equations that must be solved more explicit,

Now, defining a new matrix,

and remembering that the $UiaV$ projectors are unitary, we arrive at our final working form for the correction vectors,

The doubles correction vector takes the form of a Lyapunov equation (a special case of the Sylvester equation), which, in the present case, may be solved in closed form.^{72} This is analogous to the equations that are encountered for the update of the doubles *t*-amplitudes in rank-reduced CCSD.^{56}

When multiple electronic states are considered, as outlined previously, we solve for the states individually and solve for each new state in a subspace that is orthogonal to any previously determined states. In practice, this is accomplished by removing any component of the sigma vectors (matrix-vector products of the compressed similarity transformed Hamiltonian with trial vectors in the Davidson diagonalization) that overlaps the previously determined RR-EOM-CCSD wavefunctions,

Here, $b^$ refers to the trial vectors in the Davidson diagonalization. Furthermore, all correction vectors added to the Davidson subspace are orthogonalized against the previously determined RR-EOM-CCSD wavefunctions. These two steps have been determined to be sufficient, in practice, to avoid the variational collapse of the RR-EOM-CCSD equations to lower-lying solutions that have already been determined.

## III. COMPUTATIONAL DETAILS

Our preliminary implementation of the RR-EOM-CCSD method was developed in the LIGHTSPEED-TERACHEM rapid prototyping environment. LIGHTSPEED is a C++/Python library that leverages the highly optimized graphical processing unit (GPU) accelerated capabilities of the Terachem^{73–77} electronic structure package. In all EOM-CCSD and RR-EOM-CCSD computations, the frozen-core approximation is applied. Furthermore, density-fitting (or RI) approximations are applied to the ERIs used in the coupled-cluster computations.^{16–18} The MP2FIT auxiliary basis sets are used to construct the approximate integrals.^{78} Cartesian coordinates corresponding to every EOM-CCSD computation reported in this work can be found in the supplementary material.

## IV. RESULTS AND DISCUSSION

The effectiveness of the RR-EOM-CCSD approach outlined in the present work relies entirely on the low-rank character of the doubles amplitudes, $rijab$. More specifically, as formulated, the rank-sparsity must exist in the manner shown in Eq. (17), i.e., an *ov*/*ov* type decomposition. To determine to what extent this quantity is low rank, we calculate the lowest singlet state with EOM-CCSD and eigendecompose the resulting doubles amplitudes. The results of this analysis are presented in Fig. 1 (upper panel) for a series of linear alkenes ranging in size from ethylene to hexadecaoctaene. The doubles amplitudes are indeed observed to exhibit rank sparsity and, importantly, rank sparsity that increases with system size. The doubles amplitudes are indefinite with positive and negative wings of the eigenvalue spectrum that are approximately equal in both number and magnitude. The indefinite nature of these amplitudes follows from their leading contributions,

The next question that must be addressed is how accurately the EOM-CCSD excitation energy can be reproduced from an approximate reconstruction of the doubles amplitudes (Fig. 1, middle panel). In this test, we approximate the doubles amplitudes using all eigenvalues with an absolute value above a specified threshold and evaluate the energy variationally,

We find that the EOM-CCSD energies are recovered quickly with decreasing eigenvalue cutoffs. Accuracy better than 10^{−4}*E*_{h} (roughly 0.003 eV) can be achieved by retaining all eigenvalues with a magnitude greater than 10^{−3}. The resulting accuracy would be sufficient for most practical applications. Remarkably, this level of accuracy can be reached with less than 10% of the $rijab$ eigenvalues for a molecule as small as octatetraene; for the largest molecule we considered, less than 5% of the eigenvalues are necessary. It is important to remember that in the context of our low-rank reformulation of EOM-CCSD [using the decomposition shown in Eq. (17)], the doubles amplitudes are compressed by the square of the number of ranks retained to approximate $rijab$. Therefore, the doubles amplitudes of hexadecaoctaene may be compressed to 0.25% of their original size without significant loss of accuracy. This result clearly demonstrates the amount of compression available in the EOM-CCSD wavefunction and is encouraging for the potential success of RR-EOM-CCSD.

Although the compression demonstrated in Fig. 1 (lower panel) is impressive, it requires the full-rank solution of the EOM-CCSD equations to be available and also a rather costly eigendecomposition of the doubles amplitudes. This result motivates the development of simplified projectors, $UiaV$, that may be employed in Eq. (17). In Fig. 2, we test the projectors defined by the generator tensor defined in Eq. (33). As we have noted, this generator is inspired by the success of methods such as CC2^{66} and leaves open the possibility to form the $UiaV$ projectors with O(*N*^{4}) effort. At a first pass, we test what should be the best possible projector defined by Eq. (33), i.e., one that is constructed from the exact EOM-CCSD singles amplitudes $ria$. We find that these projectors perform well enough for practical applications. For the same set of linear alkenes considered earlier, we obtain accuracy in the excitation energy on the order of 10^{−2} eV with cutoffs on the eigenvalues of 10^{−4}. This leads to worse compression than was possible using the eigenvectors of the doubles amplitudes but still can compress the number of wavefunction parameters to 1% for the largest system considered. If subhundredth electronvolts accuracy is required, an eigenvalue cutoff of 10^{−6} appears to be sufficient; this still leads to significant compression of the doubles amplitudes to less than 4% of their original size for hexadecaoctaene. The compromise between the accuracy of the result, amount of compression achieved, and cost of formation leads us to conclude that the projector defined by Eq. (33) is suitable for practical applications of RR-EOM-CCSD.

Another remaining concern is whether or not the projectors defined by Eq. (33) can be determined on-the-fly. It is possible that feedback between errors introduced by the approximated doubles amplitudes and concomitant changes in the singles amplitudes could degrade the performance of the projector. However, as shown in Fig. 2, no such loss of accuracy is observed (compare the dashed and solid lines in the upper and lower panels). The projectors that are determined on-the-fly perform similarly to those that are constructed from the exact singles amplitudes. Their performance is similar both in terms of the accuracy afforded by the projector (Fig. 2, upper panel) as well as the number of eigenvalues above a particular cutoff (Fig. 2, lower panel). Yet another concern is whether or not the on-the-fly determination of the projectors (outlined in Algorithm 1) is stable. If this procedure requires many macroiterations (solutions of the RR-EOM-CCSD equations using a fixed projector) or if the excitation energies diverge, then this projector would not be practically useful (or, at least, the approach used to determine the projector would need to be reconsidered). In Fig. 3, we present a representative example of the convergence (represented by the change in excitation energy) as a function of macroiteration. This test is initialized with singles amplitudes taken from a CIS wavefunction. We observe stable convergence behavior, even when large, 10^{−2} or 10^{−3}, cutoffs on the eigenvalues are applied. In this test, we converge the energy change to 10^{−8} *E*_{h}; this is unnecessarily tight convergence, given the accuracy that can be achieved. The convergence behavior of RR-EOM-CCSD can be quite slow when large cutoffs are applied and becomes much more efficient as the eigenvalue cutoffs are tightened. This is not at all surprising since, in the limit where it is exact, it is obvious that RR-EOM-CCSD requires only one macroiteration (the projector is known *a priori* in this case and is simply the identity). As the cutoffs are tightened, the numerical value of the singles amplitudes becomes less important since a large number of factors are included. In practice, we find that only 3–4 macroiterations are required to reach the maximum achievable level of accuracy. Subsequent iterations are required for a converged result, but the excitation energy does not improve, that is, the changes in the energy are smaller than the deviation between RR-EOM-CCSD and EOM-CCSD. Finally, we note that it should be possible to converge the early macroiterations less tightly (i.e., the Davidson diagonalization) and converge more tightly as the energy change between macroiterations decreases. This may lead to a procedure that does not require significantly more multiplications of trial vectors with $H\xaf$ (microiterations) than does traditional EOM-CCSD, but we have yet to develop an optimal procedure that exploits this possibility.

While we have demonstrated that the RR-EOM-CCSD algorithm outlined in this work performs well for a single electronic state, we must also verify that it performs well when multiple electronic states are considered. To that end, we determine the energy of the lowest four singlet states of pyridine (Fig. 4) using RR-EOM-CCSD, as outlined in Algorithm 1. One potential issue with this approach is that errors introduced in the lower-lying excited states may magnify and lead to larger errors in the higher-lying excited states due to the orthogonalization that is applied. However, this does not seem to be a significant problem in practice. The excitation energies corresponding to the S_{3} and S_{4} excited states of pyridine are determined significantly more accurately than the excitation energies of transitions to S_{1} and S_{2}. This indicates that the accuracy of the RR-EOM-CCSD method depends more on the character of the excited state in question (not all excited states of the same molecule have identical rank in $rijab$) than it does on the order in which the states are determined. A final observation from the RR-EOM-CCSD computation on pyridine is that although the RR-EOM-CCSD energies are upper bounds of the exact EOM-CCSD energy, there is no guarantee that increasing the size of the auxiliary space will lead to lower energies. We do observe a general trend in nearly all numerical tests that tighter thresholds lead to lower, more accurate energies. However, as seen for S_{4} of pyridine, convergence to the exact EOM-CCSD energy is not necessarily monotonic with respect to a decreasing cutoff. Finally, we mention a final caveat on the variational nature of the RR-EOM-CCSD energies: the rigorous upper bound on the excitation energy only applies to the lowest excited state. Errors in the description of lower-lying states may lead to nonvariational behavior in the higher-lying states. It should be noted that if a single (i.e., state-universal) projector is applied, then RR-EOM-CCSD provides an upper bound to all excitation energies. In practice, however, we have yet to encounter a case where the RR-EOM-CCSD energy is lower than the EOM-CCSD energy for any excited state even though state-specific projectors are applied.

Critical to the success of the RR-EOM-CCSD method is the size intensivity of the excitation energy and the error in the excitation energy. There are two limits where size intensivity should be considered. One is the case of noninteracting molecules (i.e., infinitely separated molecules) where the size intensivity of the excitation energy and the error associated with the rank-reduction should be preserved exactly. The other case is what we will call “weakly interacting” systems, for example, a chromophore in solvent. In this case, the excitation will tend to localize on the chromophore, but the effect of the other molecules cannot be neglected. In Table I, we compute RR-EOM-CCSD excitation energies corresponding to the S_{1} state of butadiene in the presence of methanol molecules. We consider the noninteracting case where the methanol molecules are distant (10^{6} Å away) and the weakly interacting case where the geometries are representative of butadiene including the first solvation shell of liquid methanol. We have argued that, on formal grounds, size intensivity should be preserved exactly in RR-EOM-CCSD for the noninteracting case. Empirically, we find this to be true as well. We find that size intensivity in the excitation energy is preserved to better than 10^{−5} *E*_{h} and that the number of auxiliary indices in the rank reduced doubles amplitudes remains constant as methanol molecules are placed 10^{6} Å away from the butadiene. This clearly demonstrates that size intensivity is preserved exactly in the RR-EOM-CCSD method in the case of noninteracting molecules.

Solvent molecules . | Excitation energy (E_{h})
. | Error (E_{h})
. | N_{aux}
. | N_{ov}
. |
---|---|---|---|---|

Butadiene/methanol infinitely separated | ||||

0 | 0.247 892 | 1.89 $\xd7$ 10^{−4} | 263 | 781 |

1 | 0.247 893 | 1.91 $\xd7$ 10^{−4} | 263 | 1 980 |

2 | 0.247 893 | 1.90 $\xd7$ 10^{−4} | 263 | 3 725 |

3 | 0.247 895 | 1.92 $\xd7$ 10^{−4} | 263 | 6 016 |

4 | 0.247 893 | 1.91 $\xd7$ 10^{−4} | 263 | 8 853 |

5 | 0.247 893 | 1.91 $\xd7$ 10^{−4} | 263 | 12 236 |

6 | 0.247 898 | 1.95 $\xd7$ 10^{−4} | 263 | 16 165 |

Butadiene/methanol cluster | ||||

0 | 0.247 891 | 1.88 $\xd7$ 10^{−4} | 263 | 781 |

1 | 0.247 052 | 1.64 $\xd7$ 10^{−4} | 341 | 1 980 |

2 | 0.245 645 | 1.73 $\xd7$ 10^{−4} | 423 | 3 725 |

3 | 0.244 516 | 1.64 $\xd7$ 10^{−4} | 489 | 6 016 |

4 | 0.243 428 | 1.52 $\xd7$ 10^{−4} | 569 | 8 853 |

5 | 0.242 378 | 1.51 $\xd7$ 10^{−4} | 647 | 12 236 |

6 | 0.241 494 | 1.60 $\xd7$ 10^{−4} | 731 | 16 165 |

Solvent molecules . | Excitation energy (E_{h})
. | Error (E_{h})
. | N_{aux}
. | N_{ov}
. |
---|---|---|---|---|

Butadiene/methanol infinitely separated | ||||

0 | 0.247 892 | 1.89 $\xd7$ 10^{−4} | 263 | 781 |

1 | 0.247 893 | 1.91 $\xd7$ 10^{−4} | 263 | 1 980 |

2 | 0.247 893 | 1.90 $\xd7$ 10^{−4} | 263 | 3 725 |

3 | 0.247 895 | 1.92 $\xd7$ 10^{−4} | 263 | 6 016 |

4 | 0.247 893 | 1.91 $\xd7$ 10^{−4} | 263 | 8 853 |

5 | 0.247 893 | 1.91 $\xd7$ 10^{−4} | 263 | 12 236 |

6 | 0.247 898 | 1.95 $\xd7$ 10^{−4} | 263 | 16 165 |

Butadiene/methanol cluster | ||||

0 | 0.247 891 | 1.88 $\xd7$ 10^{−4} | 263 | 781 |

1 | 0.247 052 | 1.64 $\xd7$ 10^{−4} | 341 | 1 980 |

2 | 0.245 645 | 1.73 $\xd7$ 10^{−4} | 423 | 3 725 |

3 | 0.244 516 | 1.64 $\xd7$ 10^{−4} | 489 | 6 016 |

4 | 0.243 428 | 1.52 $\xd7$ 10^{−4} | 569 | 8 853 |

5 | 0.242 378 | 1.51 $\xd7$ 10^{−4} | 647 | 12 236 |

6 | 0.241 494 | 1.60 $\xd7$ 10^{−4} | 731 | 16 165 |

In practice, it is more important that the error in the excitation energy be size intensive in the weakly interacting case. For the weakly interacting butadiene/methanol clusters (shown in Table I), we find that the error in the RR-EOM-CCSD excitation energy remains relatively constant and tends to decrease slightly as the system size is increased. This suggests that RR-EOM-CCSD can be applied to larger systems (molecules in solution, photoactive proteins, etc.) without degradation in accuracy. These results also demonstrate that RR-EOM-CCSD is able to exploit the spatial locality of the excited state without any *a priori* definition of this locality. The iterative refinement of the projectors is a method for discovering spatial locality on-the-fly. This can be seen most clearly for the infinitely separated systems, where the size of the auxiliary space does not change. It is also evident for the largest weakly interacting butadiene methanol cluster where the excitation energy can be reproduced with accuracy better than 0.005 eV while using less than 5% of the space *ov* space and compressing the $rijab$ amplitudes to 0.2% of their original size. This should be compared to the isolated butadiene molecule, where 33% of the *ov* space was included.

After addressing the more formal concerns regarding the accuracy of RR-EOM-CCSD, we turn our focus to more practical questions of accuracy. In Table II, we compute the RR-EOM-CCSD excitation energy for the lowest two or three singlet excited states of eight common dye molecules (depicted in Fig. 5). In each case, we use eigenvalue cutoffs of 10^{−4} and 10^{−6} in the construction of the projectors. For the 10^{−4} cutoff, the excitation energies are typically accurate to within 0.02–0.04 eV. With that larger cutoff, the *ov* space is compressed to 10% or less of the full space. This level of accuracy is probably sufficient for most practical applications since the error due to the rank-reduction is an order of magnitude less than the intrinsic error in the EOM-CCSD excitation energy relative to the best theoretical estimates. This is achieved with 1% or less of the variational parameters present in the full EOM-CCSD wavefunction. Higher accuracy can be reached by decreasing the eigenvalue cutoff to 10^{−6}. In that case, the excitation energies are correct to better than 0.01 eV in all cases. The cost of this higher accuracy is an increase in the size of the auxiliary space by a factor of 2–3 relative to the 10^{−4} cutoff.

. | . | Excitation . | Error (eV) . | Number of factors . | |||
---|---|---|---|---|---|---|---|

Molecule . | State . | energy (eV) . | 10^{–4}
. | 10^{–6}
. | 10^{–4}
. | 10^{–6}
. | Full . |

Lawsone | S_{1} | 3.439 | 0.037 | 0.009 | 487 | 1301 | 5 344 |

S_{2} | 3.912 | 0.029 | 0.007 | 510 | 1330 | 5 344 | |

S_{3} | 4.166 | 0.034 | 0.008 | 558 | 1545 | 5 344 | |

Indophenol | S_{1} | 3.170 | 0.042 | 0.009 | 592 | 1681 | 7 511 |

S_{2} | 3.296 | 0.036 | 0.008 | 638 | 1781 | 7 511 | |

Alizarin | S_{1} | 3.544 | 0.032 | 0.008 | 685 | 2146 | 10 120 |

S_{2} | 3.663 | 0.020 | 0.003 | 750 | 2129 | 10 120 | |

Harmine | S_{1} | 4.426 | 0.021 | 0.005 | 788 | 2117 | 9 120 |

S_{2} | 4.892 | 0.021 | 0.005 | 814 | 2195 | 9 120 | |

Indigo | S_{1} | 2.698 | 0.020 | 0.004 | 827 | 2467 | 12 576 |

S_{2} | 3.382 | 0.035 | 0.010 | 753 | 2066 | 12 576 | |

S_{3} | 3.521 | 0.027 | 0.006 | 871 | 2496 | 12 576 | |

Acriflavine | S_{1} | 3.256 | 0.024 | 0.006 | 834 | 2324 | 10 458 |

S_{2} | 3.646 | 0.027 | 0.006 | 812 | 2234 | 10 458 | |

Methyl yellow | S_{1} | 3.063 | 0.019 | 0.006 | 555 | 1579 | 10 879 |

S_{2} | 3.904 | 0.023 | 0.005 | 788 | 2230 | 10 879 | |

Neutral red | S_{1} | 3.427 | 0.027 | 0.006 | 910 | 2548 | 13 392 |

S_{2} | 3.768 | 0.029 | 0.009 | 885 | 2532 | 13 392 | |

S_{3} | 3.785 | 0.044 | 0.011 | 767 | 2263 | 13 392 |

. | . | Excitation . | Error (eV) . | Number of factors . | |||
---|---|---|---|---|---|---|---|

Molecule . | State . | energy (eV) . | 10^{–4}
. | 10^{–6}
. | 10^{–4}
. | 10^{–6}
. | Full . |

Lawsone | S_{1} | 3.439 | 0.037 | 0.009 | 487 | 1301 | 5 344 |

S_{2} | 3.912 | 0.029 | 0.007 | 510 | 1330 | 5 344 | |

S_{3} | 4.166 | 0.034 | 0.008 | 558 | 1545 | 5 344 | |

Indophenol | S_{1} | 3.170 | 0.042 | 0.009 | 592 | 1681 | 7 511 |

S_{2} | 3.296 | 0.036 | 0.008 | 638 | 1781 | 7 511 | |

Alizarin | S_{1} | 3.544 | 0.032 | 0.008 | 685 | 2146 | 10 120 |

S_{2} | 3.663 | 0.020 | 0.003 | 750 | 2129 | 10 120 | |

Harmine | S_{1} | 4.426 | 0.021 | 0.005 | 788 | 2117 | 9 120 |

S_{2} | 4.892 | 0.021 | 0.005 | 814 | 2195 | 9 120 | |

Indigo | S_{1} | 2.698 | 0.020 | 0.004 | 827 | 2467 | 12 576 |

S_{2} | 3.382 | 0.035 | 0.010 | 753 | 2066 | 12 576 | |

S_{3} | 3.521 | 0.027 | 0.006 | 871 | 2496 | 12 576 | |

Acriflavine | S_{1} | 3.256 | 0.024 | 0.006 | 834 | 2324 | 10 458 |

S_{2} | 3.646 | 0.027 | 0.006 | 812 | 2234 | 10 458 | |

Methyl yellow | S_{1} | 3.063 | 0.019 | 0.006 | 555 | 1579 | 10 879 |

S_{2} | 3.904 | 0.023 | 0.005 | 788 | 2230 | 10 879 | |

Neutral red | S_{1} | 3.427 | 0.027 | 0.006 | 910 | 2548 | 13 392 |

S_{2} | 3.768 | 0.029 | 0.009 | 885 | 2532 | 13 392 | |

S_{3} | 3.785 | 0.044 | 0.011 | 767 | 2263 | 13 392 |

Absolute accuracy is the definitive criterion with which to judge the performance of the RR-EOM-CCSD method (relative to the full EOM-CCSD result). However, it is also of practical importance that the error is relatively constant and, more importantly, relatively smooth when the molecular geometry is changed. That is to say that one would hope to obtain smooth potential energy surfaces from the RR-EOM-CCSD method. Errors on the order of 0.02–0.04 eV might be tolerated if they were nearly constant over some regions of the potential energy surface that is of interest. A random error of the size that was uncorrelated with the molecular geometry, on the other hand, would not be tolerated and much tighter cutoffs (and subsequently larger projectors) would be required. As presently formulated, to achieve C^{0} continuity (or better) over some finite region of the potential surface, the size of the auxiliary space would need to be held constant (here, for simplicity, we will discuss only regions of the EOM-CCSD potential surface that does not contain crossing between electronic states). Discrete jumps in the size of the auxiliary space will unavoidably lead to C^{0} discontinuities. Even with the size of the auxiliary space held constant, C^{0} continuity is not guaranteed, because there are conical intersections between the eigenvalues of the projector’s generator matrix [Eq. (30)]. As a result, the surfaces defined by these eigenvalues will have C^{0} continuity but may have discontinuous derivatives (i.e., lack C^{1} continuity). If a region of the RR-EOM-CCSD potential energy surface is considered that contains conical intersections between the smallest eigenpair of the generator matrix that is included in the projector and the largest excluded eigenpair, the C^{0} continuity of the RR-EOM-CCSD potential energy surface will be lost. In subsequent discussions of smoothness, we will not refer to smoothness in a mathematically rigorous sense, but rather we will consider whether or not any discontinuities are visible to the eye.

In Figs. 6 and 7, we apply RR-EOM-CCSD to the excited-state proton transfer reactions in 10-hydroxybenzo[h]quinoline (HBQ) and 2-(2′-hydroxyphenyl)benzothiazole (HBT), respectively. These potential curves were computed with eigenvalue cutoffs starting at 10^{−2} and decreasing to 10^{−6}. Absolute errors on the order of 0.01 eV are not achieved until the cutoff is reduced to 10^{−5}; however, cutoffs at or below 10^{−3} produce visually smooth potential curves. Perhaps a more demanding test of RR-EOM-CCSD is found in Fig. 8. Here, we apply RR-EOM-CCSD to the S_{1} and S_{2} states of 4-(dimethylamino)benzonitrile (DMABN) as the dimethylamino group is twisted out of plane with the benzene ring. This torsional motion leads to a crossing between the S_{1} and S_{2} states that causes the adiabatic states to change character. Again, we find that an eigenvalue cutoff of 10^{−3} or less produces visually smooth potential curves (here, for both S_{1} and S_{2}). We also observe that the absolute error in each of the surfaces changes with the change in character. This suggests that the absolute error in the RR-EOM-CCSD method is most closely related to the character of the electronic state under consideration, i.e., the error in the diabatic representation is slowly varying. Once again, errors of 0.01 eV or better can be achieved with an eigenvalue cutoff of 10^{−5} or less. The results obtained for these scans of the potential energy surfaces are quite encouraging and suggest that RR-EOM-CCSD could be used in place of EOM-CCSD for many of its common applications.

A final question is how RR-EOM-CCSD performs in larger basis sets and those that include diffuse atomic orbitals. In Fig. 9, we apply RR-EOM-CCSD to compute the lowest four singlet excitation energies of pyridine in cc-pVDZ, cc-pVTZ, cc-pVQZ, aug-cc-pVDZ, and aug-cc-pVTZ basis sets. The performance of RR-EOM-CCSD in all the basis sets is in line with what is observed for general application in smaller basis sets (i.e., Table II). The largest errors that we observe are less than 1.5 × 10^{−2} eV. Due to the relatively small size of the molecule, the compression achieved is not as significant as what is observed for larger molecules. However, the doubles amplitudes are significantly more compressible in larger (and diffuse) basis sets. While the *ov* space can be compressed to only 50%–60% of its conventional size in a cc-pVDZ basis, this improves to 35%–45% in cc-pVTZ and aug-cc-pVDZ basis sets. Even more compression is possible (30%–35%) in cc-pVQZ and aug-cc-pVTZ basis sets. The only somewhat concerning trend is that the error in the excitation energies tends to increase with the size of the basis. At present, the origin of this effect remains unclear; it may be an artifact resulting from the small size of the molecule—more than half of the *ov* space is spanned by the rank-reduced doubles amplitudes in the cc-pVDZ basis. It is also possible that slightly tighter cutoffs are required when larger basis sets are applied. An encouraging result is that the performance in diffuse basis sets is similar to the performance in basis sets of similar size that do not include diffuse orbitals (i.e., aug-cc-pVTZ and cc-pVQZ). This suggests that the performance will not degrade in these basis sets and for states with Rydberg character (here, S_{4} of pyridine has significant Rydberg character). It should be emphasized that the performance in these larger basis sets is satisfactory and allows the target accuracy of 0.01 eV to be achieved.

## V. CONCLUSIONS

Starting from the EOM-CCSD energy functional, we have developed a rank-reduced formulation of the EOM-CCSD method. Introduction of a compressed representation of the EOM-CCSD doubles amplitudes leads to equations for the low-rank wavefunction parameters that do not require the formation of full-rank intermediates. This approach allows excitation energies to be obtained with accuracy on the order of 0.01 eV with roughly 4% of the free parameters available in the full EOM-CCSD method. The compressibility of the wavefunction coefficients is expected to increase for larger molecular systems—especially when the excitation is localized (or delocalized over a constant number of sites). To obtain a projection of the full-rank EOM-CCSD wavefunction into a reduced space, we use the eigenvectors of the approximate doubles amplitudes obtained from a procedure similar to that of CC2. This allows the leading contributions to the EOM-CCSD doubles amplitudes to be captured, while also allowing the eigenvectors to be obtained efficiently. Our RR-EOM-CCSD method provides accurate treatments of multiple low-lying excitations simultaneously. We have demonstrated that the number of free parameters in the EOM-CCSD wavefunction can be reduced from O(*N*^{4}) to O(*N*^{2}) while retaining suitable accuracy for the most common applications of the EOM-CCSD method.

The performance of our RR-EOM-CCSD method—in terms of both accuracy and compression—motivates the development of a fully optimized software implementation that exploits all of the available flexibility in the factorization. As we demonstrated in our previous work on RR-CCSD, a rank-reduced variant of the linearized coupled-cluster doubles method (LCCD) can be implemented with O(*N*^{5}) scaling, provided that an RI factorization of the ERIs is also applied.^{56} This result should transfer directly into RR-EOM-CCSD since both methods are linear in their excitation operators. Following a CCSD computation with O(*N*^{6}) scaling, the EOM-CCSD excitation energies could be evaluated with O(*N*^{5}) effort, making the determination of numerous excited states practical for any system where ground state CCSD computations are possible. Efforts in this direction are currently underway. Owing to the reduction in the size of the wavefunction and the potential for a sub-O(*N*^{6}) implementation, we believe that RR-EOM-CCSD is a promising direction for the accurate determination of electronic excitation energies.

## SUPPLEMENTARY MATERIAL

See the supplementary material for Cartesian coordinates of all molecular geometries used in this work.

## ACKNOWLEDGMENTS

This material is based on the work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Scientific Discovery through Advanced Computing (SciDAC) program.

## REFERENCES

^{+}, CO, and H

_{2}O