Rank reduced coupled cluster theory. II. Equation-of-motion coupled-cluster singles and doubles

Hohenstein, Edward G.; Zhao, Yao; Parrish, Robert M.; Martínez, Todd J.

doi:10.1063/1.5121867

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⁶), 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⁴), 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.²¹ 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).²¹ 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.¹⁵ 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).³⁵ 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.³⁶

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.⁴³ 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).⁴⁸ 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.⁵⁶ 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.

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

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,

Ĥ e^{\hat{T}} |Φ_{0}⟩ = E_{CC} e^{\hat{T}} |Φ_{0}⟩ .

(1)

Here, $\hat{T}$ is the cluster excitation operator,

\hat{T} = {\hat{T}}_{1} + {\hat{T}}_{2} + \dots = t_{i}^{a} a_{a}^{†} a_{i} + \frac{1}{4} t_{i j}^{a b} a_{a}^{†} a_{b}^{†} a_{j} a_{i} + \dots .

(2)

After left-multiplication by $e^{- \hat{T}}$ ⁠,

e^{- \hat{T}} Ĥ e^{\hat{T}} |Φ_{0}⟩ = E_{CC} |Φ_{0}⟩,

(3)

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

⟨Φ_{0}| e^{- \hat{T}} Ĥ e^{\hat{T}} |Φ_{0}⟩ = E_{CC},

(4)

⟨) Φ_{i j \dots}^{a b \dots} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ_{0} ⟩) = 0 .

(5)

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.⁵⁹ 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.⁶⁰ 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.,

e^{- \hat{T}} Ĥ e^{\hat{T}} {\hat{R}}^{(n)} |Φ_{0}⟩ = E_{EOM - CC}^{(n)} {\hat{R}}^{(n)} |Φ_{0}⟩ .

(6)

Here, the ${\hat{R}}^{(n)}$ operator is a CI-like linear excitation operator (for the nth electronic state),

{\hat{R}}^{(n)} = r_{0} + r_{i}^{a} a_{a}^{†} a_{i} + \frac{1}{4} r_{i j}^{a b} a_{a}^{†} a_{b}^{†} a_{j} a_{i} + \dots .

(7)

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

⟨Φ_{0}| {\hat{L}}^{(n)} e^{- \hat{T}} Ĥ e^{\hat{T}} = ⟨Φ_{0}| {\hat{L}}^{(n)} E_{EOM - CC}^{(n)},

(8)

where ${\hat{L}}^{(n)}$ is a linear de-excitation operator corresponding to the n-th state,

{\hat{L}}^{(n)} = l_{0} + l_{i}^{a} a_{i}^{†} a_{a} + \frac{1}{4} l_{i j}^{a b} a_{i}^{†} a_{j}^{†} a_{b} a_{a} + \dots .

(9)

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

E_{E O M - C C}^{(n)} = \frac{⟨Φ_{0}| {\hat{L}}^{(n)} e^{- \hat{T}} Ĥ e^{\hat{T}} {\hat{R}}^{(n)} |Φ_{0}⟩}{⟨Φ_{0}| {\hat{L}}^{(n)} {\hat{R}}^{(n)} |Φ_{0}⟩},

(10)

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

δ_{m n} = ⟨Φ_{0}| {\hat{L}}^{(m)} {\hat{R}}^{(n)} |Φ_{0}⟩ .

(11)

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

\frac{\partial E_{E O M - C C}^{(n)}}{\partial r_{i j \dots}^{a b \dots}} = 0,

(12)

\frac{\partial E_{E O M - C C}^{(n)}}{\partial l_{i j \dots}^{a b \dots}} = 0 .

(13)

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.⁷ Furthermore, amplitudes of the $\hat{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,

{\hat{R}}^{(0)} = 1 .

(14)

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,⁶³ by left-multiplying the (right-handed) ground-state EOM-CC eigenvalue problem [Eq. (3)] by the ${\hat{R}}^{(n)}$ operator and subtracting it from the (right-handed) excited-state EOM-CC eigenvalue problem [Eq. (6)],

\begin{matrix} e^{- \hat{T}} Ĥ e^{\hat{T}} {\hat{R}}^{(n)} |Φ_{0}⟩ - {\hat{R}}^{(n)} e^{- \hat{T}} Ĥ e^{\hat{T}} |Φ_{0}⟩ \\ = E_{E O M - C C}^{(n)} {\hat{R}}^{(n)} |Φ_{0}⟩ - E_{E O M - C C}^{(0)} {\hat{R}}^{(n)} |Φ_{0}⟩, \end{matrix}

(15)

[e^{- \hat{T}} Ĥ e^{\hat{T}}, {\hat{R}}^{(n)}] |Φ_{0}⟩ = (E_{E O M - C C}^{(n)} - E_{E O M - C C}^{(0)}) {\hat{R}}^{(n)} |Φ_{0}⟩ .

(16)

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⁶), and the storage requirements (for the doubles amplitudes, $r_{i j}^{a b}$ ⁠) scale with the fourth-power of system size, O(N⁴). In the present work, we will reformulate the EOM-CCSD method to exploit a low-rank factorization of the doubles amplitudes of the form

r_{i j}^{a b} = U_{i a}^{V} R^{V W} U_{j b}^{W} .

(17)

Here, the $U_{i a}^{V}$ 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⁴) to O(N²) by replacing the $r_{i j}^{a b}$ tensor with the R^VW matrix. At present, we will work exclusively with unitary projectors, viz.,

δ_{V W} = U_{i a}^{V} U_{i a}^{W} .

(18)

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 ${\hat{R}}^{(n)}$ and ${\hat{L}}^{(n)}$ operators,

{\hat{R}}^{(n)} = r_{0} + r_{i}^{a} a_{a}^{†} a_{i} + \frac{1}{4} U_{i a}^{V} R^{V W} U_{j b}^{W} a_{a}^{†} a_{b}^{†} a_{j} a_{i},

(19)

{\hat{L}}^{(n)} = l_{0} + l_{i}^{a} a_{i}^{†} a_{a} + \frac{1}{4} U_{i a}^{V} L^{V W} U_{j b}^{W} a_{i}^{†} a_{j}^{†} a_{b} a_{a} .

(20)

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,

\begin{matrix} L (r_{0}, r_{i}^{a}, R^{V W}, l_{0}, l_{i}^{a}, L^{V W}; t_{i}^{a}, t_{i j}^{a b}, U_{i a}^{V}) = \frac{⟨Φ_{0}| {\hat{L}}^{(n)} e^{- \hat{T}} Ĥ e^{\hat{T}} {\hat{R}}^{(n)} |Φ_{0}⟩}{⟨Φ_{0}| {\hat{L}}^{(n)} {\hat{R}}^{(n)} |Φ_{0}⟩} . \end{matrix}

(21)

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,

\frac{\partial L}{\partial l_{0}} = 0; \frac{\partial L}{\partial l_{i}^{a}} = 0; \frac{\partial L}{\partial L^{V W}} = 0 .

(22)

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²v² × o²v²), the eigenvectors are found using a nonsymmetric variant of Davidson diagonalization.⁶⁴ For EOM-CCSD, the equations that must be satisfied can be written in the following form:

⟨) Φ_{0} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ_{0} ⟩) r_{0} + ⟨) Φ_{0} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ_{k}^{c} ⟩) r_{k}^{c} + ⟨) Φ_{0} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ_{k l}^{c d} ⟩) r_{k l}^{c d} = E_{E O M - C C} \cdot r_{0},

(23)

⟨) Φ_{i}^{a} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ_{0} ⟩) r_{0} + ⟨) Φ_{i}^{a} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ_{k}^{c} ⟩) r_{k}^{c} + ⟨) Φ_{i}^{a} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ_{k l}^{c d} ⟩) r_{k l}^{c d} = E_{E O M - C C} \cdot r_{i}^{a},

(24)

⟨) Φ_{i j}^{a b} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ_{0} ⟩) r_{0} + ⟨) Φ_{i j}^{a b} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ_{k}^{c} ⟩) r_{k}^{c} + ⟨) Φ_{i j}^{a b} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ_{k l}^{c d} ⟩) r_{k l}^{c d} = E_{E O M - C C} \cdot r_{i j}^{a b} .

(25)

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,

|) Φ^{V W} ⟩) = |) Φ_{i j}^{a b} ⟩) U_{i a}^{V} U_{j b}^{W} .

(26)

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,

⟨) Φ_{0} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ_{0} ⟩) r_{0} + ⟨) Φ_{0} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ_{k}^{c} ⟩) r_{k}^{c} + ⟨) Φ_{0} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ^{X Y} ⟩) R^{X Y} = E_{EOM - CC} \cdot r_{0},

(27)

⟨) Φ_{i}^{a} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ_{0} ⟩) r_{0} + ⟨) Φ_{i}^{a} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ_{k}^{c} ⟩) r_{k}^{c} + ⟨) Φ_{i}^{a} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ^{X Y} ⟩) R^{X Y} = E_{EOM - CC} \cdot r_{i}^{a},

(28)

⟨) Φ^{V W} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ_{0} ⟩) r_{0} + ⟨) Φ^{V W} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ_{k}^{c} ⟩) r_{k}^{c} + ⟨) Φ^{V W} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ^{X Y} ⟩) R^{X Y} = E_{EOM - CC} \cdot R^{V W} .

(29)

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 $U_{i a}^{V}$ ⁠. 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 $U_{i a}^{V}$ 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⁵⁶ 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.⁵⁸ 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₀ = 1, $r_{i}^{a}$ = 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, $⟨) Φ_{i}^{a} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ_{0} ⟩)$ and $⟨) Φ^{V W} |) e^{- \hat{T}} Ĥ e^{\hat{T}} |) Φ_{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.⁶⁵ 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, $U_{i a}^{V}$ ⁠. These projectors can be defined as a subset of the eigenvectors of some generating matrix,

G_{i j}^{a b} = U_{i a}^{V} g^{V} U_{j b}^{V} .

(30)

The $G_{i j}^{a b}$ generator is a fourth-order tensor that is treated as a ov × ov matrix for the purpose of eigendecomposition. The low-rank projectors, $U_{i a}^{V}$ ⁠, are the subset of the eigenvectors chosen from the eigenvalues with the largest magnitude. That is, we introduce some threshold, $ε$ ⁠, 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,

G_{i j}^{a b} = r_{i j}^{a b} .

(31)

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 $r_{i j}^{a b}$ amplitudes without requiring knowledge of the exact EOM-CCSD wavefunction.

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

ω \cdot r_{i j}^{a b} = ⟨) Φ_{i j}^{a b} |) [\hat{F}, {\hat{R}}_{2}] + [Ĥ, {\hat{R}}_{1}] |) Φ_{0} ⟩),

(32)

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

G_{i j}^{a b} = \frac{⟨) Φ_{i j}^{a b} |) [Ĥ, {\hat{R}}_{1}] |) Φ_{0} ⟩)}{ϵ_{i} + ϵ_{j} - ϵ_{a} - ϵ_{b} + ω},

(33)

that requires as input only a reasonable approximation for the $r_{i}^{a}$ amplitudes and the excitation energy, ω. In terms of both accuracy and practicality, the best source of $r_{i}^{a}$ 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 $r_{i}^{a}$ 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).⁶⁷ 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 $r_{i}^{a}$ 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 $r_{i}^{a}$ amplitudes. The process of solving the eigenvalue problem and updating the projector is repeated until the $r_{i}^{a}$ 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)⁶⁸ 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⁴) effort, provided that a low-rank decomposition of the electron repulsion integrals of the form

(p q |r s)) \approx {\tilde{L}}_{p q}^{A} {\tilde{L}}_{r s}^{A}

(34)

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

{\tilde{R}}_{i a}^{A} = r_{i}^{b} {\tilde{L}}_{b a}^{A} - r_{j}^{a} {\tilde{L}}_{i j}^{A}

(35)

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

\begin{matrix} \frac{1}{ε_{i} + ε_{j} - ε_{a} - ε_{b} + ω} & = & - \int_{0}^{\infty} e^{- [ε_{a} + ε_{b} - ε_{i} - ε_{j} - ω] t} d t \\ \approx & - \sum_{w} λ_{w} e^{- α_{w} [ε_{a} + ε_{b} - ε_{i} - ε_{j} - ω]} = τ_{i}^{w} τ_{j}^{w} τ_{a}^{w} τ_{b}^{w} \end{matrix}

(36)

(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,

G_{i j}^{a b} = - τ_{i}^{w} τ_{j}^{w} τ_{a}^{w} τ_{b}^{w} ({\tilde{L}}_{i a}^{A} {\tilde{R}}_{j b}^{A} + {\tilde{R}}_{i a}^{A} {\tilde{L}}_{j b}^{A}) .

(37)

We may further simplify this expression by writing a ov × 2N_DF supermatrix containing the ${\tilde{R}}_{i a}^{A}$ and ${\tilde{L}}_{i a}^{A}$ quantities (treated as $o v \times N_{DF}$ matrices),

Y_{i a}^{w B} = τ_{i}^{w} τ_{a}^{w} [{\tilde{L}}_{i a}^{A} {\tilde{R}}_{i a}^{A}],

(38)

Z_{i a}^{w B} = τ_{i}^{w} τ_{a}^{w} [{\tilde{R}}_{i a}^{A} {\tilde{L}}_{i a}^{A}] .

(39)

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

G_{i j}^{a b} = - Y_{i a}^{w B} Z_{j b}^{w B},

(40)

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

G_{i k}^{a c} G_{k j}^{c b} = U_{i a}^{V} {(g^{V})}^{2} U_{j b}^{V} .

(41)

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 $G_{i j}^{a b}$ matrix is symmetric, we may expand the square of the generator as

G_{i k}^{a c} G_{k j}^{c b} = Y_{i a}^{w B} Z_{k c}^{w B} Z_{k c}^{w^{'} B^{'}} Y_{j b}^{w^{'} B^{'}} .

(42)

We treat the Z tensor as an ov × 2N_DFN_w matrix and express it in terms of its singular value decomposition (SVD),

Z_{i a}^{w B} = U_{i a}^{V} s^{V} {\tilde{V}}_{w B}^{V} .

(43)

Inserting this decomposition into Eq. (41),

G_{i k}^{a c} G_{k j}^{c b} = Y_{i a}^{w B} {\tilde{V}}_{w B}^{V^{'}} s^{V^{'}} \underset{δ_{V^{'} W^{'}}}{\underset{︸}{U_{k c}^{V^{'}} U_{k c}^{W^{'}}}} s^{W^{'}} {\tilde{V}}_{w^{'} B^{'}}^{W^{'}} Y_{j b}^{w^{'} B^{'}},

(44)

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

{\tilde{Z}}_{i a}^{V} = Y_{i a}^{w B} {\tilde{V}}_{w B}^{V} s^{V},

(45)

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

G_{i k}^{a c} G_{k j}^{c b} = {\tilde{Z}}_{i a}^{V} {\tilde{Z}}_{j b}^{V} .

(46)

Now, the projector (eigenvectors of $G_{i j}^{a b}$ ⁠) may be obtained via an SVD of the $\tilde{Z}$ tensor,

{\tilde{Z}}_{i a}^{V} = U_{i a}^{W} s^{W} {\bar{V}}_{V}^{W} .

(47)

The $U_{i a}^{V}$ singular vectors are the eigenvectors of $G_{i j}^{a b}$ ⁠, and the singular values are the absolute values of the $G_{i j}^{a b}$ 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,

\begin{matrix} [\begin{matrix} r^{(n)} & R^{⊥} \end{matrix}] \bar{H} [\begin{gathered} r^{(n)} \\ R^{⊥} \end{gathered}] = [\begin{matrix} E^{(n)} & * \\ 0 & {\bar{H}}^{⊥} \end{matrix}] . \end{matrix}

(48)

Here, r⁽ⁿ⁾ are the vectors corresponding to the nth eigenvalue of the similarity transformed Hamiltonian, $\bar{H}$ ⁠, and R^⊥ is a subspace that is orthogonal to the eigenvectors that have been determined. The rⁿ⁺¹ vector is obtained as an eigenvector of ${\bar{H}}^{⊥}$ ⁠. This process can be repeated until ${\bar{H}}^{⊥}$ 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

{\bar{H}}_{1} = \bar{H} - σ r_{1} v^{†} .

(49)

Here, ${\bar{H}}_{1}$ is the deflated similarity transformed Hamiltonian, σ is a shift parameter, r₁ is the first right-hand eigenvector of the similarity transformed Hamiltonian, and v is an arbitrary vector chosen such that v^†r₁ = 1. One may expect that v should be chosen to be the left-hand eigenvector of the similarity transformed Hamiltonian corresponding to r₁; 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₁, 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 $r_{i}^{a}$ 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 $r_{i}^{a}$ amplitudes; several RR-EOM-CCSD states can then be determined using a single projector constructed from the state-averaged amplitudes. The $r_{i}^{a}$ amplitudes from these RR-EOM-CCSD wavefunctions also serve as a good initial guess for the self-consistent determination of state-specific projectors.

ALGORITHM 1

An algorithm for solving for RR-EOM-CCSD electronic states.

1:

Solve HF equations

2:

Solve CCSD equations

3:

for n States do

4:

Initialize: k=0

5:

Construct guess

r_{i}^{a, (n)}

amplitudes

6:

Form

G_{i j}^{a b, (n)}

using

r_{i}^{a, (n)}

⁠: Eq.

7:

Form

U_{i a}^{V, (n)}

by diagonalizing

r_{i j}^{a b, (n)}

⁠: Eq. ()

8:

if n > 0 then

9:

R^{⊥} = Î - \sum_{n = 0,1 \dots n - 1} {\hat{R}}^{(n)} |Φ_{0}⟩ ⟨Φ_{0}| {\hat{R}}^{(n) †}

10:

end if

11:

Solve RR-EOM-CCSD: Eqs.

12:

Set:

E_{k}^{(n)} = E_{RR - EOM - CC}^{(n)}

13:

Update

r_{i}^{a, (n)}

amplitudes

14:

loop

15:

Rebuild

G_{i j}^{a b, (n)}

using

r_{i}^{a, (n)}

⁠: Eq.

Rebuild

U_{i a}^{V, (n)}

by diagonalizing

G_{i j}^{a b, (n)}

⁠: Eq. ()

16:

if n > 0 then

17:

R^{⊥} = Î - \sum_{n = 0,1 \dots n - 1} {\hat{R}}^{(n)} |Φ_{0}⟩ ⟨Φ_{0}| {\hat{R}}^{(n) †}

18:

end if

19:

Solve RR-EOM-CCSD: Eqs.

20:

Set:

E_{k}^{(n)} = E_{RR - EOM - CC}^{(n)}

21:

Update

r_{i}^{a, (n)}

amplitudes

22:

if

|E_{k + 1}^{(n)} - E_{k}^{(n)}|

then

23:

return

E_{RR - EOM - CC}^{(n)}

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.⁶⁴ The formation of the residual follows from Eqs. (27)–(29). In practice, we solve for the excitation energy directly in terms of $r_{i}^{a}$ and R^VW and determine r₀ 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

{\tilde{δ}}_{i}^{a} = - \frac{ρ_{i}^{a}}{f_{a a} - f_{i i} - \tilde{ω}},

(50)

{\tilde{δ}}_{i j}^{a b} = - \frac{ρ_{i j}^{a b}}{f_{a a} + f_{b b} - f_{i i} - f_{j j} - \tilde{ω}},

(51)

where $\tilde{ω}$ 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,

{\tilde{δ}}^{V W} = - {[⟨Φ^{V W}| \hat{F} |Φ^{X Y}⟩ - δ_{V X} δ_{W Y} \tilde{ω}]}^{- 1} ρ^{X Y} .

(52)

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

U_{i a}^{X} U_{j b}^{Y} (f_{a a} + f_{b b} - f_{i i} - f_{j j} - \tilde{ω}) U_{i a}^{V} U_{j b}^{W} {\tilde{δ}}^{V W} = - ρ^{X Y} .

(53)

Now, defining a new matrix,

A^{V W} = U_{i a}^{V} (f_{a a} - f_{i i} - \frac{1}{2} \tilde{ω}) U_{i a}^{W},

(54)

and remembering that the $U_{i a}^{V}$ projectors are unitary, we arrive at our final working form for the correction vectors,

A^{V X} {\tilde{δ}}^{X W} + {\tilde{δ}}^{V X} A^{X W} = - ρ^{V W} .

(55)

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.⁷² This is analogous to the equations that are encountered for the update of the doubles t-amplitudes in rank-reduced CCSD.⁵⁶

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,

σ = \bar{H} \hat{b} |Φ_{0}⟩ - \sum_{n = 0,1, \dots n - 1} {\hat{R}}^{(n)} |Φ_{0}⟩ ⟨Φ_{0}| {\hat{R}}^{(n) †} \bar{H} \hat{b} |Φ_{0}⟩ .

(56)

Here, $\hat{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.⁷⁸ 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, $r_{i j}^{a b}$ ⁠. 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,

\begin{matrix} E_{EOM - C C S D} \cdot r_{i j}^{a b} = ⟨a b| |c j⟩ r_{i}^{c} + ⟨a b| |i c⟩ r_{j}^{c} - ⟨k b| |i j⟩ r_{k}^{a} - ⟨a k| |i j⟩ r_{k}^{b} + \dots . \end{matrix}

(57)

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,

E_{EOM - CCSD} = \frac{⟨Φ_{0}| {\hat{R}}^{(n) †} e^{- \hat{T}} Ĥ e^{\hat{T}} {\hat{R}}^{(n)} |Φ_{0}⟩}{⟨Φ_{0}| {\hat{R}}^{(n) †} {\hat{R}}^{(n)} |Φ_{0}⟩} .

(58)

We find that the EOM-CCSD energies are recovered quickly with decreasing eigenvalue cutoffs. Accuracy better than 10⁻⁴E_h (roughly 0.003 eV) can be achieved by retaining all eigenvalues with a magnitude greater than 10⁻³. The resulting accuracy would be sufficient for most practical applications. Remarkably, this level of accuracy can be reached with less than 10% of the $r_{i j}^{a b}$ 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 $r_{i j}^{a b}$ ⁠. 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.

FIG. 1.

$FIG. 1. The absolute value of the eigenvalues of the rijab amplitudes are plotted using a relative index where the negative eigenvalue with the largest magnitude has an index of 0 and the largest positive eigenvalue has a relative index of 1 (upper panel). This analysis is performed on the first singlet excited state for a series of linear alkenes of increasing size. The error in the EOM-CCSD excitation energy is shown when the rijab amplitudes are approximately reconstructed using eigenpairs that are above a cutoff on the magnitude of the eigenvalue (middle panel). The fraction of the ov space that is retained for a particular choice of cutoff is shown in the lower panel. All computations are performed at the EOM-CCSD/cc-pVDZ level of theory.$

View large Download slide

The absolute value of the eigenvalues of the $r_{i j}^{a b}$ amplitudes are plotted using a relative index where the negative eigenvalue with the largest magnitude has an index of 0 and the largest positive eigenvalue has a relative index of 1 (upper panel). This analysis is performed on the first singlet excited state for a series of linear alkenes of increasing size. The error in the EOM-CCSD excitation energy is shown when the $r_{i j}^{a b}$ amplitudes are approximately reconstructed using eigenpairs that are above a cutoff on the magnitude of the eigenvalue (middle panel). The fraction of the ov space that is retained for a particular choice of cutoff is shown in the lower panel. All computations are performed at the EOM-CCSD/cc-pVDZ level of theory.

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, $U_{i a}^{V}$ ⁠, 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⁶⁶ and leaves open the possibility to form the $U_{i a}^{V}$ projectors with O(N⁴) 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 $(r_{i}^{a})$ ⁠. 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⁻² eV with cutoffs on the eigenvalues of 10⁻⁴. 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⁻⁶ 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.

FIG. 2.

$FIG. 2. The accuracy of the rijab=ΦijabĤ,R^‘1Φ0/ϵa+ϵb−ϵi−ϵj−ω amplitudes (shown diagrammatically in the upper panel) for use in the construction of a projector for RR-EOM-CCSD is assessed for the lowest singlet excited states of a series of linear alkenes. The projector of this type constructed using the exact ria amplitudes (solid lines) is compared to the projector that is formed iteratively (dashed lines). The accuracy of the RR-EOM-CCSD/cc-pVDZ excitation energy relative to the exact EOM-CCSD/cc-pVDZ is shown in the upper panel, and the fraction of the ov space included in the RR-EOM-CCSD computation is shown in the lower panel.$

View large Download slide

The accuracy of the $r_{i j}^{a b} = ⟨Φ_{i j}^{a b} |[Ĥ, {\hat{R}}_{‘ 1}]| Φ_{0}⟩ / (ϵ_{a} + ϵ_{b} - ϵ_{i} - ϵ_{j} - ω)$ amplitudes (shown diagrammatically in the upper panel) for use in the construction of a projector for RR-EOM-CCSD is assessed for the lowest singlet excited states of a series of linear alkenes. The projector of this type constructed using the exact $r_{i}^{a}$ amplitudes (solid lines) is compared to the projector that is formed iteratively (dashed lines). The accuracy of the RR-EOM-CCSD/cc-pVDZ excitation energy relative to the exact EOM-CCSD/cc-pVDZ is shown in the upper panel, and the fraction of the ov space included in the RR-EOM-CCSD computation is shown in the lower panel.

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⁻² or 10⁻³, cutoffs on the eigenvalues are applied. In this test, we converge the energy change to 10⁻⁸ 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 $\bar{H}$ (microiterations) than does traditional EOM-CCSD, but we have yet to develop an optimal procedure that exploits this possibility.

FIG. 3.

View large Download slide

Convergence of the fixed-point iterations to simultaneously optimize the RR-EOM-CCSD excited-state wavefunction and projector of the $r_{i j}^{a b}$ amplitudes. The main panel shows the energy change between iterations using a variety of cutoffs, $ϵ$ ⁠, on the magnitude of the eigenvalues. The procedure was assigned as converged when the energy change became smaller than 10⁻⁸ E_h. The inset shows the error in the RR-EOM-CCSD excitation energy relative to the full EOM-CCSD result during the iterative procedure. These computations were performed for the lowest singlet electronic state of the hexatriene molecule at the EOM-CCSD/cc-pVDZ level of theory.

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₃ and S₄ excited states of pyridine are determined significantly more accurately than the excitation energies of transitions to S₁ and S₂. 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 $r_{i j}^{a b}$ ⁠) 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₄ 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.

FIG. 4.

View large Download slide

Accuracy of RR-EOM-CCSD/cc-pVDZ relative to the full EOM-CCSD computation for the lowest four singlet states of pyridine as a function of the cutoff on the eigenvalues of the projector defined in Eq. (33).

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₁ state of butadiene in the presence of methanol molecules. We consider the noninteracting case where the methanol molecules are distant (10⁶ Å 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⁻⁵ E_h and that the number of auxiliary indices in the rank reduced doubles amplitudes remains constant as methanol molecules are placed 10⁶ Å 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.

TABLE I.

Excitation energies of butadiene/methanol clusters corresponding to the first excited singlet state computed at the RR-EOM-CCSD/cc-pVDZ level of theory using a cutoff of 10⁻⁴ in the construction of the doubles amplitude projector. Errors are reported relative to the full RR-EOM-CCSD/cc-pVDZ excitation energy.

Solvent molecules	Excitation energy (E_h)	Error (E_h)	N_aux	N_ov
Butadiene/methanol infinitely separated
0	0.247 892	1.89 $\times$ 10⁻⁴	263	781
1	0.247 893	1.91 $\times$ 10⁻⁴	263	1 980
2	0.247 893	1.90 $\times$ 10⁻⁴	263	3 725
3	0.247 895	1.92 $\times$ 10⁻⁴	263	6 016
4	0.247 893	1.91 $\times$ 10⁻⁴	263	8 853
5	0.247 893	1.91 $\times$ 10⁻⁴	263	12 236
6	0.247 898	1.95 $\times$ 10⁻⁴	263	16 165
Butadiene/methanol cluster
0	0.247 891	1.88 $\times$ 10⁻⁴	263	781
1	0.247 052	1.64 $\times$ 10⁻⁴	341	1 980
2	0.245 645	1.73 $\times$ 10⁻⁴	423	3 725
3	0.244 516	1.64 $\times$ 10⁻⁴	489	6 016
4	0.243 428	1.52 $\times$ 10⁻⁴	569	8 853
5	0.242 378	1.51 $\times$ 10⁻⁴	647	12 236
6	0.241 494	1.60 $\times$ 10⁻⁴	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 $\times$ 10⁻⁴	263	781
1	0.247 893	1.91 $\times$ 10⁻⁴	263	1 980
2	0.247 893	1.90 $\times$ 10⁻⁴	263	3 725
3	0.247 895	1.92 $\times$ 10⁻⁴	263	6 016
4	0.247 893	1.91 $\times$ 10⁻⁴	263	8 853
5	0.247 893	1.91 $\times$ 10⁻⁴	263	12 236
6	0.247 898	1.95 $\times$ 10⁻⁴	263	16 165
Butadiene/methanol cluster
0	0.247 891	1.88 $\times$ 10⁻⁴	263	781
1	0.247 052	1.64 $\times$ 10⁻⁴	341	1 980
2	0.245 645	1.73 $\times$ 10⁻⁴	423	3 725
3	0.244 516	1.64 $\times$ 10⁻⁴	489	6 016
4	0.243 428	1.52 $\times$ 10⁻⁴	569	8 853
5	0.242 378	1.51 $\times$ 10⁻⁴	647	12 236
6	0.241 494	1.60 $\times$ 10⁻⁴	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 $r_{i j}^{a b}$ 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⁻⁴ and 10⁻⁶ in the construction of the projectors. For the 10⁻⁴ 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⁻⁶. 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⁻⁴ cutoff.

TABLE II.

Excitation energies of the lowest few singlet states of several common dye molecules computed at the RR-EOM-CCSD/cc-pVDZ level of theory. Structures of the dye molecules are shown in Fig. 5.

		Excitation	Error (eV)		Number of factors
Molecule	State	energy (eV)	10^–4	10^–6	10^–4	10^–6	Full
Lawsone	S₁	3.439	0.037	0.009	487	1301	5 344
	S₂	3.912	0.029	0.007	510	1330	5 344
	S₃	4.166	0.034	0.008	558	1545	5 344
Indophenol	S₁	3.170	0.042	0.009	592	1681	7 511
	S₂	3.296	0.036	0.008	638	1781	7 511
Alizarin	S₁	3.544	0.032	0.008	685	2146	10 120
	S₂	3.663	0.020	0.003	750	2129	10 120
Harmine	S₁	4.426	0.021	0.005	788	2117	9 120
	S₂	4.892	0.021	0.005	814	2195	9 120
Indigo	S₁	2.698	0.020	0.004	827	2467	12 576
	S₂	3.382	0.035	0.010	753	2066	12 576
	S₃	3.521	0.027	0.006	871	2496	12 576
Acriflavine	S₁	3.256	0.024	0.006	834	2324	10 458
	S₂	3.646	0.027	0.006	812	2234	10 458
Methyl yellow	S₁	3.063	0.019	0.006	555	1579	10 879
	S₂	3.904	0.023	0.005	788	2230	10 879
Neutral red	S₁	3.427	0.027	0.006	910	2548	13 392
	S₂	3.768	0.029	0.009	885	2532	13 392
	S₃	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₁	3.439	0.037	0.009	487	1301	5 344
	S₂	3.912	0.029	0.007	510	1330	5 344
	S₃	4.166	0.034	0.008	558	1545	5 344
Indophenol	S₁	3.170	0.042	0.009	592	1681	7 511
	S₂	3.296	0.036	0.008	638	1781	7 511
Alizarin	S₁	3.544	0.032	0.008	685	2146	10 120
	S₂	3.663	0.020	0.003	750	2129	10 120
Harmine	S₁	4.426	0.021	0.005	788	2117	9 120
	S₂	4.892	0.021	0.005	814	2195	9 120
Indigo	S₁	2.698	0.020	0.004	827	2467	12 576
	S₂	3.382	0.035	0.010	753	2066	12 576
	S₃	3.521	0.027	0.006	871	2496	12 576
Acriflavine	S₁	3.256	0.024	0.006	834	2324	10 458
	S₂	3.646	0.027	0.006	812	2234	10 458
Methyl yellow	S₁	3.063	0.019	0.006	555	1579	10 879
	S₂	3.904	0.023	0.005	788	2230	10 879
Neutral red	S₁	3.427	0.027	0.006	910	2548	13 392
	S₂	3.768	0.029	0.009	885	2532	13 392
	S₃	3.785	0.044	0.011	767	2263	13 392

FIG. 5.

View large Download slide

Structures of the dye molecules included in the tests of RR-EOM-CCSD shown in Table II: (a) lawsone, (b) indophenol, (c) alizarin, (d) harmine, (e) indigo, (f) acriflavine, (g) methyl yellow, and (h) neutral red.

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⁰ 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⁰ discontinuities. Even with the size of the auxiliary space held constant, C⁰ 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⁰ continuity but may have discontinuous derivatives (i.e., lack C¹ 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⁰ 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⁻² and decreasing to 10⁻⁶. Absolute errors on the order of 0.01 eV are not achieved until the cutoff is reduced to 10⁻⁵; however, cutoffs at or below 10⁻³ 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₁ and S₂ 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₁ and S₂ states that causes the adiabatic states to change character. Again, we find that an eigenvalue cutoff of 10⁻³ or less produces visually smooth potential curves (here, for both S₁ and S₂). 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⁻⁵ 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.

FIG. 6.

View large Download slide

Excited-state proton transfer potential of the 10-hydroxybenzo[h]quinoline (HBQ) molecule on the S₁ electronic state computed at the EOM-CCSD/cc-pVDZ level of theory. RR-EOM-CCSD/cc-pVDZ energies are computed using decreasing cutoffs $(ϵ)$ in the construction of the ${\hat{R}}_{2}$ amplitude projector. The top panel shows the energy (in eV) relative to the S₀ enol minimum. The lower panel shows the error (in eV) of RR-EOM-CCSD relative to EOM-CCSD. The proton transfer coordinate, R_PT, is defined as the difference between the O–H and N–H bond distances: R_PT = R_OH − R_NH.

FIG. 7.

View large Download slide

Excited-state proton transfer potential of the 2-(2′-hydroxyphenyl)benzothiazole (HBT) molecule on the S₁ electronic state computed at the EOM-CCSD/cc-pVDZ level of theory. RR-EOM-CCSD/cc-pVDZ energies are computed using decreasing cutoffs $(ϵ)$ in the construction of the ${\hat{R}}_{2}$ amplitude projector. The top panel shows the energy (in eV) relative to the S₀ enol minimum. The lower panel shows the error (in eV) of RR-EOM-CCSD relative to EOM-CCSD. The proton transfer coordinate, R_PT, is defined as the difference between the O–H and N–H bond distances: R_PT = R_OH − R_NH.

FIG. 8.

View large Download slide

Torsional potential of the S₁ and S₂ electronic states of 4-(dimethylamino)benzonitrile (DMABN) computed at the EOM-CCSD/cc-pVDZ level of theory. RR-EOM-CCSD/cc-pVDZ energies are computed using decreasing cutoffs $(ϵ)$ in the construction of the ${\hat{R}}_{2}$ amplitude projector. The top panel shows the energy (in eV) relative to S₀ with a torsional angle of 0°. The lower panel shows the error (in eV) of RR-EOM-CCSD relative to EOM-CCSD.

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⁻² 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₄ 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.

FIG. 9.

View large Download slide

Excitation energies of the lowest four singlet excited states of pyridine computed with EOM-CCSD in several basis sets. The upper panel reports the error in excitation energy of RR-EOM-CCSD relative to EOM-CCSD using a cutoff of $ϵ$ = 10⁻⁶. The lower panel reports the compression of the doubles amplitudes for each state. The compression is calculated as a ratio of the number of factors used to construct the RR-EOM-CCSD doubles amplitudes divided by the number of factors in the full ov space.

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⁴) to O(N²) 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⁵) scaling, provided that an RI factorization of the ERIs is also applied.⁵⁶ 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⁶) scaling, the EOM-CCSD excitation energies could be evaluated with O(N⁵) 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⁶) 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

1.

H. J.

Monkhorst

, “

Calculation of properties with the coupled-cluster method

,”

Int. J. Quantum Chem.

12

,

421

(

1977

).

https://doi.org/10.1002/qua.560120850

Google Scholar

Crossref

2.

D.

Mukherjee

and

P. K.

Mukherjee

, “

A response-function approach to the direct calculation of the transition-energy in a multiple-cluster expansion formalism

,”

Chem. Phys.

39

,

325

(

1979

).

https://doi.org/10.1016/0301-0104(79)80153-6

Google Scholar

Crossref

3.

K.

Emrich

, “

An extension of the coupled cluster formalism to excited states (I)

,”

Nucl. Phys. A

351

,

379

(

1981

).

https://doi.org/10.1016/0375-9474(81)90179-2

Google Scholar

Crossref

4.

S.

Ghosh

and

D.

Mukherjee

, “

Use of cluster expansion techniques in quantum chemistry. A linear response model for calculating energy differences

,”

Proc. - Indian Acad. Sci., Chem. Sci.

93

,

947

(

1984

).

https://doi.org/10.1007/BF02840339

Google Scholar

Crossref

5.

H.

Sekino

and

R. J.

Bartlett

, “

A linear response, coupled-cluster theory for excitation energy

,”

Int. J. Quantum Chem.

26

,

255

(

1984

).

https://doi.org/10.1002/qua.560260826

Google Scholar

Crossref

6.

J.

Geertsen

,

M.

Rittby

, and

R. J.

Bartlett

, “

The equation-of-motion coupled-cluster method: Excitation energies of Be and CO

,”

Chem. Phys. Lett.

164

,

57

(

1989

).

https://doi.org/10.1016/0009-2614(89)85202-9

Google Scholar

Crossref

7.

J. F.

Stanton

and

R. J.

Bartlett

, “

The equation of motion coupled-cluster method. A systematic biorthogonal approach to molecular-excitation energies, transition-probabilities, and excited-state properties

,”

J. Chem. Phys.

98

,

7029

(

1993

).

https://doi.org/10.1063/1.464746

Google Scholar

Crossref

8.

H.

Koch

,

O.

Christiansen

,

R.

Kobayashi

,

P.

Jørgensen

, and

T.

Helgaker

, “

A direct atomic orbital driven implementation of the coupled-cluster singles and doubles (CCSD) model

,”

Chem. Phys. Lett.

228

,

233

(

1994

).

https://doi.org/10.1016/0009-2614(94)00898-1

Google Scholar

Crossref

9.

O.

Christiansen

,

H.

Koch

,

A.

Halkier

,

P.

Jørgensen

,

T.

Helgaker

, and

A.

Sánchez de Merás

, “

Large-scale calculations of excitation energies in coupled cluster theory: The singlet excited states of benzene

,”

J. Chem. Phys.

105

,

6921

(

1996

).

https://doi.org/10.1063/1.471985

Google Scholar

Crossref

10.

H.

Koch

,

A. S.

de Merás

,

T.

Helgaker

, and

O.

Christiansen

, “

The integral-direct coupled cluster singles and doubles model

,”

J. Chem. Phys.

104

,

4157

(

1996

).

https://doi.org/10.1063/1.471227

Google Scholar

Crossref

11.

N. H. F.

Beebe

and

J.

Linderberg

, “

Simplifications in the generation and transformation of two-electron integrals in molecular calculations

,”

Int. J. Quantum Chem.

12

,

683

(

1977

).

https://doi.org/10.1002/qua.560120408

Google Scholar

Crossref

12.

I.

Røeggen

and

E.

Wisløff-Nilssen

, “

On the Beebe-Linderberg two-electron integral approximation

,”

Chem. Phys. Lett.

132

,

154

(

1986

).

https://doi.org/10.1016/0009-2614(86)80099-9

Google Scholar

Crossref

13.

H.

Koch

,

A.

Sánchez de Merás

, and

T. B.

Pedersen

, “

Reduced scaling in electronic structure calculations using Cholesky decompositions

,”

J. Chem. Phys.

118

,

9481

(

2003

).

https://doi.org/10.1063/1.1578621

Google Scholar

Crossref

14.

F.

Aquilante

,

T. B.

Pedersen

, and

R.

Lindh

, “

Low-cost evaluation of the exchange Fock matrix from Cholesky and density fitting representations of the electron repulsion integrals

,”

J. Chem. Phys.

126

,

194106

(

2007

).

https://doi.org/10.1063/1.2736701

Google Scholar

Crossref

PubMed

15.

S. D.

Folkestad

,

E. F.

Kjønstad

, and

H.

Koch

, “

An efficient algorithm for Cholesky decomposition of electron repulsion integrals

,”

J. Chem. Phys.

150

,

194112

(

2019

).

https://doi.org/10.1063/1.5083802

Google Scholar

Crossref

PubMed

16.

J. L.

Whitten

, “

Coulombic potential energy integrals and approximations

,”

J. Chem. Phys.

58

,

4496

(

1973

).

https://doi.org/10.1063/1.1679012

Google Scholar

Crossref

17.

B. I.

Dunlap

,

J. W. D.

Connolly

, and

J. R.

Sabin

, “

On the applicability of LCAO-Xα methods to molecules containing transition metal atoms: The nickel atom and nickel hydride

,”

Int. J. Quantum Chem.

12

,

81

(

1977

).

https://doi.org/10.1002/qua.560120813

Google Scholar

Crossref

18.

B. I.

Dunlap

,

J. W. D.

Connolly

, and

J. R.

Sabin

, “

On some approximations in applications of Xα theory

,”

J. Chem. Phys.

71

,

3396

(

1979

).

https://doi.org/10.1063/1.438728

Google Scholar

Crossref

19.

M.

Feyereisen

,

G.

Fitzgerald

, and

A.

Komornicki

, “

Use of approximate integrals in ab initio theory. An application in MP2 energy calculations

,”

Chem. Phys. Lett.

208

,

359

(

1993

).

https://doi.org/10.1016/0009-2614(93)87156-w

Google Scholar

Crossref

20.

O.

Vahtras

,

J.

Almlöf

, and

M. W.

Feyereisen

, “

Integral approximations for LCAO-SCF calculations

,”

Chem. Phys. Lett.

213

,

514

(

1993

).

https://doi.org/10.1016/0009-2614(93)89151-7

Google Scholar

Crossref

21.

E.

Epifanovsky

,

D.

Zuev

,

X.

Feng

,

K.

Khistyaev

,

Y.

Shao

, and

A. I.

Krylov

, “

General implementation of the resolution-of-the-identity and Cholesky representations of electron repulsion integrals within coupled-cluster and equation-of-motion methods: Theory and benchmarks

,”

J. Chem. Phys.

139

,

134105

(

2013

).

https://doi.org/10.1063/1.4820484

Google Scholar

Crossref

PubMed

22.

R. A.

Friesner

, “

Solution of the Hartree–Fock equations by a pseudospectral method: Application to diatomic molecules

,”

J. Chem. Phys.

85

,

1462

(

1986

).

https://doi.org/10.1063/1.451237

Google Scholar

Crossref

23.

R. A.

Friesner

, “

Solution of the Hartree–Fock equations for polyatomic molecules by a pseudospectral method

,”

J. Chem. Phys.

86

,

3522

(

1987

).

https://doi.org/10.1063/1.451955

Google Scholar

Crossref

24.

J. M.

Langlois

,

R. P.

Muller

,

T. R.

Coley

,

W. A.

Goddard

,

M. N.

Ringnalda

,

Y.

Won

, and

R. A.

Friesner

, “

Pseudospectral generalized valence-bond calculations: Application to methylene, ethylene, and silylene

,”

J. Chem. Phys.

92

,

7488

(

1990

).

https://doi.org/10.1063/1.458184

Google Scholar

Crossref

25.

M. N.

Ringnalda

,

M.

Belhadj

, and

R. A.

Friesner

, “

Pseudospectral Hartree–Fock theory: Applications and algorithmic improvements

,”

J. Chem. Phys.

93

,

3397

(

1990

).

https://doi.org/10.1063/1.458819

Google Scholar

Crossref

26.

R. A.

Friesner

, “

New methods for electronic structure calculations on large molecules

,”

Annu. Rev. Phys. Chem.

42

,

341

(

1991

).

https://doi.org/10.1146/annurev.physchem.42.1.341

Google Scholar

Crossref

PubMed

27.

T. J.

Martinez

,

A.

Mehta

, and

E. A.

Carter

, “

Pseudospectral full configuration interaction

,”

J. Chem. Phys.

97

,

1876

(

1992

).

https://doi.org/10.1063/1.463176

Google Scholar

Crossref

28.

T. J.

Martinez

and

E. A.

Carter

, “

Pseudospectral double excitation configuration interaction

,”

J. Chem. Phys.

98

,

7081

(

1993

).

https://doi.org/10.1063/1.464751

Google Scholar

Crossref

29.

T. J.

Martinez

and

E. A.

Carter

, “

Pseudospectral Moller–Plesset perturbation theory through third order

,”

J. Chem. Phys.

100

,

3631

(

1994

).

https://doi.org/10.1063/1.466350

Google Scholar

Crossref

30.

T. J.

Martinez

and

E. A.

Carter

, “

Pseudospectral multireference single and double excitation configuration interaction

,”

J. Chem. Phys.

102

,

7564

(

1995

).

https://doi.org/10.1063/1.469088

Google Scholar

Crossref

31.

T. J.

Martinez

and

E. A.

Carter

, in

Modern Electronic Structure Theory

, edited by

D. R.

Yarkony

(

World Scientific

,

Singapore

,

1995

), p.

1132

.

Google Scholar

Crossref

32.

E. G.

Hohenstein

,

R. M.

Parrish

, and

T. J.

Martinez

, “

Tensor hypercontraction density fitting. I. Quartic scaling second- and third-order Moller-Plesset perturbation theory

,”

J. Chem. Phys.

137

,

044103

(

2012

).

https://doi.org/10.1063/1.4732310

Google Scholar

Crossref

PubMed

33.

E. G.

Hohenstein

,

R. M.

Parrish

,

C. D.

Sherrill

, and

T. J.

Martinez

, “

Communication: Tensor hypercontraction. III. Least-squares tensor hypercontraction for the determination of correlated wavefunctions

,”

J. Chem. Phys.

137

,

221101

(

2012

).

https://doi.org/10.1063/1.4768241

Google Scholar

Crossref

PubMed

34.

R. M.

Parrish

,

E. G.

Hohenstein

,

T. J.

Martinez

, and

C. D.

Sherrill

, “

Tensor hypercontraction. II. Least-squares renormalization

,”

J. Chem. Phys.

137

,

224106

(

2012

).

https://doi.org/10.1063/1.4768233

Google Scholar

Crossref

PubMed

35.

R. M.

Parrish

,

C. D.

Sherrill

,

E. G.

Hohenstein

,

S. I. L.

Kokkila

, and

T. J.

Martinez

, “

Communication: Acceleration of coupled cluster singles and doubles via orbital-weighted least-squares tensor hypercontraction

,”

J. Chem. Phys.

140

,

181102

(

2014

).

https://doi.org/10.1063/1.4876016

Google Scholar

Crossref

PubMed

36.

A. K.

Dutta

,

F.

Neese

, and

R.

Izsak

, “

Speeding up equation of motion coupled cluster theory with the chain of spheres approximation

,”

J. Chem. Phys.

144

,

034102

(

2016

).

https://doi.org/10.1063/1.4939844

Google Scholar

Crossref

PubMed

37.

A. K.

Dutta

,

F.

Neese

, and

R.

Izsak

, “

Accelerating the coupled-cluster singles and doubles method using the chain-of-sphere approximation

,”

Mol. Phys.

116

,

1428

(

2018

).

https://doi.org/10.1080/00268976.2017.1416201

Google Scholar

Crossref

38.

F.

Neese

,

A.

Hansen

, and

D. G.

Liakos

, “

Efficient and accurate approximations to the local coupled cluster singles doubles method using a truncated pair natural orbital basis

,”

J. Chem. Phys.

131

,

064103

(

2009

).

https://doi.org/10.1063/1.3173827

Google Scholar

Crossref

PubMed

39.

C.

Riplinger

and

F.

Neese

, “

An efficient and near linear scaling pair natural orbital based local coupled cluster method

,”

J. Chem. Phys.

138

,

034106

(

2013

).

https://doi.org/10.1063/1.4773581

Google Scholar

Crossref

PubMed

40.

M.

Sparta

and

F.

Neese

, “

Chemical applications carried out by local pair natural orbital based coupled-cluster methods

,”

Chem. Soc. Rev.

43

,

5032

(

2014

).

https://doi.org/10.1039/c4cs00050a

Google Scholar

Crossref

PubMed

41.

M.

Schwilk

,

Q. L.

Ma

,

C.

Koppl

, and

H. J.

Werner

, “

Scalable electron correlation methods. 3. Efficient and accurate parallel local coupled cluster with pair natural orbitals (PNO-LCCSD)

,”

J. Chem. Theory Comput.

13

,

3650

(

2017

).

https://doi.org/10.1021/acs.jctc.7b00554

Google Scholar

Crossref

PubMed

42.

Q. L.

Ma

and

H. J.

Werner

, “

Explicitly correlated local coupled-cluster methods using pair natural orbitals

,”

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

8

,

e1371

(

2018

).

https://doi.org/10.1002/wcms.1371

Google Scholar

Crossref

43.

A. K.

Dutta

,

F.

Neese

, and

R.

Izsak

, “

Towards a pair natural orbital coupled cluster method for excited states

,”

J. Chem. Phys.

145

,

034102

(

2016

).

https://doi.org/10.1063/1.4958734

Google Scholar

Crossref

PubMed

44.

A. K.

Dutta

,

M.

Nooijen

,

F.

Neese

, and

R.

Izsak

, “

Automatic active space selection for the similarity transformed equations of motion coupled cluster method

,”

J. Chem. Phys.

146

,

074103

(

2017

).

https://doi.org/10.1063/1.4976130

Google Scholar

Crossref

PubMed

45.

L. M. J.

Huntington

,

M.

Krupicka

,

F.

Neese

, and

R.

Izsak

, “

Similarity transformed equation of motion coupled-cluster theory based on an unrestricted Hartree-Fock reference for applications to high-spin open-shell systems

,”

J. Chem. Phys.

147

,

174104

(

2017

).

https://doi.org/10.1063/1.5001320

Google Scholar

Crossref

PubMed

46.

A. K.

Dutta

,

M.

Nooijen

,

F.

Neese

, and

R.

Izsak

, “

Exploring the accuracy of a low scaling similarity transformed equation of motion method for vertical excitation energies

,”

J. Chem. Theory Comput.

14

,

72

(

2018

).

https://doi.org/10.1021/acs.jctc.7b00802

Google Scholar

Crossref

PubMed

47.

A. K.

Dutta

,

M.

Saitow

,

C.

Riplinger

,

F.

Neese

, and

R.

Izsak

, “

A near-linear scaling equation of motion coupled cluster method for ionized states

,”

J. Chem. Phys.

148

,

244101

(

2018

).

https://doi.org/10.1063/1.5029470

Google Scholar

Crossref

PubMed

48.

M.

Nooijen

and

R. J.

Bartlett

, “

Similarity transformed equation-of-motion coupled-cluster theory: Details, examples, and comparisons

,”

J. Chem. Phys.

107

,

6812

(

1997

).

https://doi.org/10.1063/1.474922

Google Scholar

Crossref

49.

C.

Peng

,

M. C.

Clement

, and

E. F.

Valeev

, “

State-averaged pair natural orbitals for excited states: A route toward efficient equation of motion coupled-cluster

,”

J. Chem. Theory Comput.

14

,

5597

(

2018

).

https://doi.org/10.1021/acs.jctc.8b00171

Google Scholar

Crossref

PubMed

50.

T.

Kinoshita

,

O.

Hino

, and

R. J.

Bartlett

, “

Singular value decomposition approach for the approximate coupled-cluster method

,”

J. Chem. Phys.

119

,

7756

(

2003

).

https://doi.org/10.1063/1.1609442

Google Scholar

Crossref

51.

O.

Hino

,

T.

Kinoshita

, and

R. J.

Bartlett

, “

Singular value decomposition applied to the compression of T-3 amplitude for the coupled cluster method

,”

J. Chem. Phys.

121

,

1206

(

2004

).

https://doi.org/10.1063/1.1763575

Google Scholar

Crossref

PubMed

52.

U.

Benedikt

,

K. H.

Bohm

, and

A. A.

Auer

, “

Tensor decomposition in post-Hartree-Fock methods. II. CCD implementation

,”

J. Chem. Phys.

139

,

224101

(

2013

).

https://doi.org/10.1063/1.4833565

Google Scholar

Crossref

PubMed

53.

U.

Benedikt

,

H.

Auer

,

M.

Espig

,

W.

Hackbusch

, and

A. A.

Auer

, “

Tensor representation techniques in post-Hartree-Fock methods: Matrix product state tensor format

,”

Mol. Phys.

111

,

2398

(

2013

).

https://doi.org/10.1080/00268976.2013.798433

Google Scholar

Crossref

54.

E.

Neuscamman

, “

Communication: A Jastrow factor coupled cluster theory for weak and strong electron correlation

,”

J. Chem. Phys.

139

,

181101

(

2013

).

https://doi.org/10.1063/1.4829536

Google Scholar

Crossref

PubMed

55.

R.

Schutski

,

J. M.

Zhao

,

T. M.

Henderson

, and

G. E.

Scuseria

, “

Tensor-structured coupled cluster theory

,”

J. Chem. Phys.

147

,

184113

(

2017

).

https://doi.org/10.1063/1.4996988

Google Scholar

Crossref

PubMed

56.

R. M.

Parrish

,

Y.

Zhao

,

E. G.

Hohenstein

, and

T. J.

Martinez

, “

Rank-reduced coupled-cluster theory I. Ground state energies and wavefunctions

,”

J. Chem. Phys.

150

,

164118

(

2019

).

https://doi.org/10.1063/1.5092505

Google Scholar

Crossref

PubMed

57.

H.

Koch

,

H. J. A.

Jensen

,

P.

Jørgensen

, and

T.

Helgaker

, “

Excitation-energies from the coupled cluster singles and doubles linear response function (CCSDLR). Applications to Be, CH⁺, CO, and H₂O

,”

J. Chem. Phys.

93

,

3345

(

1990

).

https://doi.org/10.1063/1.458815

Google Scholar

Crossref

58.

H.

Koch

and

P.

Jørgensen

, “

Coupled cluster response functions

,”

J. Chem. Phys.

93

,

3333

(

1990

).

https://doi.org/10.1063/1.458814

Google Scholar

Crossref

59.

H.

Koch

,

R.

Kobayashi

,

A.

Sanchez de Merás

, and

P.

Jørgensen

, “

Calculation of size-intensive transition moments from the coupled cluster singles and doubles linear response function

,”

J. Chem. Phys.

100

,

4393

(

1994

).

https://doi.org/10.1063/1.466321

Google Scholar

Crossref

60.

A. I.

Krylov

, “

Equation-of-motion coupled-cluster methods for open-shell and electronically excited species: The hitchhiker’s guide to Fock space

,”

Annu. Rev. Phys. Chem.

59

,

433

(

2008

).

https://doi.org/10.1146/annurev.physchem.59.032607.093602

Google Scholar

Crossref

PubMed

61.

J. F.

Stanton

, “

Many-body methods for excited state potential energy surfaces. I. General theory of energy gradients for the equation-of-motion coupled-cluster method

,”

J. Chem. Phys.

99

,

8840

(

1993

).

https://doi.org/10.1063/1.465552

Google Scholar

Crossref

62.

J. F.

Stanton

and

J.

Gauss

, “

Analytic energy gradients for the equation-of-motion coupled-cluster method: Implementation and application to the HCN/HNC system

,”

J. Chem. Phys.

100

,

4695

(

1994

).

https://doi.org/10.1063/1.466253

Google Scholar

Crossref

63.

I.

Shavitt

and

R. J.

Bartlett

,

Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory

, Cambridge Molecular Science (

Cambridge University Press

,

Cambridge, New York

,

2009

).

Google Scholar

Crossref

64.

E. R.

Davidson

, “

The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real-symmetric matrices

,”

J. Comput. Phys.

17

,

87

(

1975

).

https://doi.org/10.1016/0021-9991(75)90065-0

Google Scholar

Crossref

65.

B.

Peng

,

R.

Van Beeumen

,

D. B.

Williams-Young

,

K.

Kowalski

, and

C.

Yang

, “

Approximate Green’s function coupled cluster method employing effective dimension reduction

,”

J. Chem. Theory Comput.

15

,

3185

(

2019

).

https://doi.org/10.1021/acs.jctc.9b00172

Google Scholar

Crossref

PubMed

66.

O.

Christiansen

,

H.

Koch

, and

P.

Jorgensen

, “

The 2nd-order approximate coupled-cluster singles and doubles model CC2

,”

Chem. Phys. Lett.

243

,

409

(

1995

).

https://doi.org/10.1016/0009-2614(95)00841-q

Google Scholar

Crossref

67.

P.

Baudin

and

K.

Kristensen

, “

Correlated natural transition orbital framework for low-scaling excitation energy calculations (CorNFLEx)

,”

J. Chem. Phys.

146

,

214114

(

2017

).

https://doi.org/10.1063/1.4984820

Google Scholar

Crossref

PubMed

68.

J. B.

Foresman

,

M.

Head-Gordon

,

J. A.

Pople

, and

M. J.

Frisch

, “

Toward a systematic molecular orbital theory for excited states

,”

J. Phys. Chem.

96

,

135

(

1992

).

https://doi.org/10.1021/j100180a030

Google Scholar

Crossref

69.

J.

Almlof

, “

Elimination of energy denominators in Moller-Plesset perturbation-theory by a Laplace transform approach

,”

Chem. Phys. Lett.

181

,

319

(

1991

).

https://doi.org/10.1016/0009-2614(91)80078-c

Google Scholar

Crossref

70.

M.

Haser

and

J.

Almlof

, “

Laplace transform techniques in Moller-Plesset perturbation-theory

,”

J. Chem. Phys.

96

,

489

(

1992

).

https://doi.org/10.1063/1.462485

Google Scholar

Crossref

71.

Y.

Saad

,

Numerical Methods for Large Eigenvalue Problems

(

Society for Industrial and Applied Mathematics

,

2011

).

Google Scholar

Crossref

72.

A. M.

Lyapunov

, “

The general problem of the stability of motion

,”

Int. J. Control

55

,

531

(

1992

).

https://doi.org/10.1080/00207179208934253

Google Scholar

Crossref

73.

I. S.

Ufimtsev

and

T. J.

Martínez

, “

Quantum chemistry on graphical processing units. 1. Strategies for two-electron integral evaluation

,”

J. Chem. Theory Comput.

4

,

222

(

2008

).

https://doi.org/10.1021/ct700268q

Google Scholar

Crossref

PubMed

74.

I. S.

Ufimtsev

and

T. J.

Martinez

, “

Quantum chemistry on graphical processing units. 2. Direct self-consistent-field implementation

,”

J. Chem. Theory Comput.

5

,

1004

(

2009

).

https://doi.org/10.1021/ct800526s

Google Scholar

Crossref

PubMed

75.

I. S.

Ufimtsev

and

T. J.

Martinez

, “

Quantum chemistry on graphical processing units. 3. Analytical energy gradients, geometry optimization, and first principles molecular dynamics

,”

J. Chem. Theory Comput.

5

,

2619

(

2009

).

https://doi.org/10.1021/ct9003004

Google Scholar

Crossref

PubMed

76.

N.

Luehr

,

I. S.

Ufimtsev

, and

T. J.

Martínez

, “

Dynamic precision for electron repulsion integral evaluation on graphical processing units (GPUs)

,”

J. Chem. Theory Comput.

7

,

949

(

2011

).

https://doi.org/10.1021/ct100701w

Google Scholar

Crossref

PubMed

77.

A. V.

Titov

,

I. S.

Ufimtsev

,

N.

Luehr

, and

T. J.

Martinez

, “

Generating efficient quantum chemistry codes for novel architectures

,”

J. Chem. Theory Comput.

9

,

213

(

2013

).

https://doi.org/10.1021/ct300321a

Google Scholar

Crossref

PubMed

78.

F.

Weigend

,

A.

Köhn

, and

C.

Hättig

, “

Efficient use of the correlation consistent basis sets in resolution of the identity MP2 calculations

,”

J. Chem. Phys.

116

,

3175

(

2002

).

https://doi.org/10.1063/1.1445115

Google Scholar

Crossref

2019

Author(s)

Rank reduced coupled cluster theory. II. Equation-of-motion coupled-cluster singles and doubles

I. INTRODUCTION

II. THEORY

A. Notation

B. Canonical EOM-CCSD

C. Rank-reduced formulation of EOM-CCSD

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

E. Multiple electronic states in RR-EOM-CCSD

F. Implementation details of RR-EOM-CCSD

III. COMPUTATIONAL DETAILS

IV. RESULTS AND DISCUSSION

V. CONCLUSIONS

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

REFERENCES

Supplementary Material

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Rank reduced coupled cluster theory. II. Equation-of-motion coupled-cluster singles and doubles

I. INTRODUCTION

II. THEORY

A. Notation

B. Canonical EOM-CCSD

C. Rank-reduced formulation of EOM-CCSD

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

E. Multiple electronic states in RR-EOM-CCSD

F. Implementation details of RR-EOM-CCSD

III. COMPUTATIONAL DETAILS

IV. RESULTS AND DISCUSSION

V. CONCLUSIONS

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

REFERENCES

Supplementary Material

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only