Approximating electronically excited states with equation-of-motion linear coupled-cluster theory

Optical spectroscopy is an important and ubiquitous tool in modern experimental studies. Improvements over the years in the use of lasers in optical spectroscopy has brought the field to an impressive level of accuracy and precision that is difficult to match using modern theoretical methods. Even so, theory has an important role to play in optical spectroscopy. By assisting in excited state assignments and predicting the density of excited states expected in a given energy range, computational results can be a powerful tool for spectroscopists.

Widespread use of time dependent density functional theory¹ (TD-DFT) for the calculation of excited states will continue for a long time to come due to the attractive low computational scaling inherent to the method. However, reliability and accuracy problems² in TD-DFT keep the theory from replacing the more computationally expensive but highly accurate multireference configuration interaction³ (MRCI) or equation of motion coupled-cluster^4,5 (EOM-CC) theories. The most commonly used form of EOM-CC theory is the equation-of-motion coupled-cluster singles and doubles (EOM-CCSD) variant,⁶ which limits ground and excited state excitations to singles and doubles only. This approximation is qualitatively consistent and approaches quantitative accuracy in many cases where single excitations are dominant. For doubly excited states, perturbative⁷ (T) or iterative⁸ (T-n) triple excitations are required^9,10 to obtain quantitative agreement with experiment.

The development of approximate excited state methods based on many-body perturbation and coupled-cluster theory has been an active and fruitful endevour since the emergence of the field.¹¹ Many approaches make use of a mean-field starting reference based on configuration interaction singles¹² and add perturbative corrections to that reference.^13–18 Other approaches use a many-body treatment of propagator theory such as the algebraic-diagrammatic construction^19,20 (ADC) scheme or by using approximate coupled-cluster linear-response^21–23 theory with the CCn methods.²⁴ Approximate methods based on equation-of-motion theory can either add a post hoc EOM-CC perturbative correction^7,25,26 for triples or an a priori perturbative approximation.^27–29 We will use perturbation theory here to develop a new perturbative EOM-CC method that also includes infinite-order effects in the ground state wavefunction. The goal will be to provide a consistent and accurate way to obtain excitation energies from either a linear coupled-cluster doubles (LCCD) or coupled-cluster singles and doubles (LCCSD) ground state wavefunction. Thus providing a route toward obtaining electronically excited states from alternative ansatz ground state wavefunctions such as the related molecular cluster perturbation theory embedding method³⁰ which dissociates to the LCCD limit.

In this work, we will use coupled-cluster perturbation theory³¹ (CCPT) to derive a general EOM-CCPT framework in Sec. II. The various EOM-CC approximations will then be systematically derived from EOM-CCPT by a choice in the Hamiltonian perturbation partitioning. Numerical results of these approximations compared against very accurate theoretical and experimental spectra are presented in Sec. IV.

Central to coupled-cluster theory⁴ is the similarity transformation of the electronic Hamiltonian by the exponential wave operator,

The excitation cluster operators (here limited to single and double excitations) are given by

and the normal ordered electronic Hamiltonian is

Here f is the one-electron Fock matrix, $v ̄$ the antisymmetrized two-electron integrals, {⋯} denotes normal ordering of the enclosed operators, and the one- (⁠$F ˆ$⁠) and two-particle (⁠$W ˆ$⁠) are defined appropriately. Throughout we reserve the indices i, j, k, … and a, b, c, … for occupied and virtual orbitals, respectively, while the indices p, q, … may refer to either occupied or virtual orbitals. Because the Hamiltonian only contains one- and two-particle operators, the Baker-Campbell-Hausdorff expansion of Eq. (1) naturally truncates after four commutators giving

Starting with an appropriate single-reference mean-field reference, |ϕ₀〉, the coupled-cluster Schrödinger equation is expressed with Eq. (5) as

The CC correlation energy given by projecting on the left by the reference function

while the necessary cluster amplitudes are obtained by solving the equations generated by projecting on the left with the auxiliary space

where 〈ϕ_g| is the set of g-fold excited determinants

The conventional approach to computing excited states through the equation-of-motion method,⁴ which directly computes the kth excitation energy (ω_k) relative to the ground state coupled-cluster wavefunction, is also obtained with Eq. (5) as

Here, ω_k is the kth general eigenvalue, and $R ˆ ( k )$ (⁠$L ˆ ( k )$⁠) is the kth right (left) general eigenvector solution to the non-Hermitian matrix, $H ̄$⁠. The 〈…〉_C denotes that only fully connected contractions are included. These eigenvector solutions are defined as the linear excitation and de-excitation operators (limited again to only singles and doubles),

with the biorthogonalization constraint

Note that because $R ˆ ( k )$ is an excitation operator, it will commute with the ground state cluster amplitudes,

As is always possible in Hamiltonian based perturbation theory we can represent the complete Hamiltonian in terms of a zeroth order Hamiltonian $H ˆ 0$ and perturbation $V ˆ$

where α → 1 is an order parameter used for convenience. Equation (17) can be directly inserted into any coupled-cluster similarity transformation or energy functional and expanded to the desired order in α. Examples other than the usual form given by Eq. (5) would include the Hermitian^32,33 (⁠$e T ˆ † H ˆ e T ˆ$⁠) or unitary³⁴ (⁠$e τ ˆ † H ˆ e τ ˆ$⁠, where $τ ˆ = T ˆ − T ˆ †$⁠) expansions, but they do not terminate without truncation.

To facilitate the decomposition into a reference and perturbation operator, we rearrange the Hamiltonian into particle excitation rank form

where the [n] superscript denotes that the operator changes particle rank from right to left by n. From this point the CC perturbation framework used is the same as we have employed in the past.^18,30,31 Here the exponential operator $e T ˆ$ is expanded to the appropriate order in terms of nth order cluster operators $T ˆ ( n )$⁠. The computation of these perturbative cluster operators is done order by order using standard perturbation theory. This is easily illustrated by expanding Eq. (5) using the partitioned Hamiltonian of Eq. (17). The similarity transformed Hamiltonian to arbitrary order in α is given by

Throughout this work we denote the exclusive order in α of an operator with the superscript (n) and we use {n} to denote an expansion inclusive of all orders in α up to n. With this notation, the general nth order similarity transformed Hamiltonian can be expressed as a compact series

The nth order ground state CCPT cluster amplitudes are completely defined from insertion of the exclusive similarity transformed Hamiltonian, Eq. (19), into (8),

Using these nth order ground state amplitudes in Eq. (7) will give the n + 1 order correlation energy

beginning with n = 2. From this the exact coupled-cluster correlation energy, relative to the mean-field energy

is

or equivalently

To obtain excitation energies relative to the nth order correlation energy without perturbatively expanding the $L ˆ$ and $R ˆ$ operators, we use the nth order inclusive similarity transformed Hamiltonian (Eq. (20)) in Eq. (10) rather than the exclusive Eq. (19), giving

Contributions from the nth order wavefunction to the spectrum obtained from Eq. (27) are possible. However, any such nth order amplitude will be neglected if it does not contribute to the corresponding nth order energy.

We have used the particle rank conserving partitioning of the Hamiltonian,

in previous work^18,30,31 to great success. Alternatively we can also choose a partitioning scheme that introduces linear coupling between cluster amplitudes, or the traditional Møller-Plesset (MP) partitioning which we will discuss next.

The standard approach in electronic structure theory is to use generalized many-body perturbation theory (GMBPT),³⁵ based on the MP partitioning of the Hamiltonian, given by

Using Eqs. (30) and (31) in Eq. (19), the first two orders of $H ̄ ( n )$ are, respectively,

and

Inserting Eq. (32) into Eq. (21) gives the quintessential GMBPT(1) wavefunction,

from which

The GMBPT(2) wavefunction is obtained from Eq. (33) similarly. The inclusive second-order similarity transformed Hamiltonian,

is the GMBPT generalization of previous^27–29 MBPT approximations to $H ̄$⁠. The spin-orbital equations for the $H ̄ { 2 }$ matrix elements^27,28 using Einstein notation are

where P(pq) is the anti-permutation operator

The difference between GMBPT and MP theory occurs when $f i a$⁠, $f j i$⁠, and $f b a ≠0$ as they would be for canonical Hartree-Fock (HF) orbitals. The first introduces a first-order singles wavefunction (Eq. (34)), while the second and third terms have to be summed to all order as in coupled-cluster theory or better transformed away by a semi-canonical occupied-occupied and virtual-virtual rotation.

The first-order similarity transformed Hamiltonian using the alternative particle rank conserving partitioning of the Hamiltonian (denoted PR0 for brevity, see Eqs. (28) and (29)) is

Using Eq. (55) in the amplitude projection Equation (21) yields the decoupled, linear-amplitude wavefunction equations

The PR0 wavefunction is infinite-order in $W ˆ [ 0 ]$⁠, which is an advantage for describing potential energy surfaces^31,36 and systematically including intermolecular many-body perturbations.³⁰ The second-order similarity transformed Hamiltonian is

from which the second-order amplitudes are computed. Starting from a canonical Hartree-Fock reference eliminates $T ˆ 1 ( 1 )$ in Eqs. (55) and (58), yielding the LCCD similarity transformed Hamiltonian

In this form, the matrix elements of Eq. (59) are identical to Eq. (37).

In order to add linear coupling between cluster operators of different rank (e.g. coupling singles and doubles) the alternative choice in Hamiltonian partitioning (denoted PR1),

can be used. From Eqs. (60) and (61) the first- and second-order similarity transformed Hamiltonians are, in operator form, the same as Eqs. (55) and (58). The motivation for this choice is evident when examining the first-order amplitude equations

As can be seen, the PR1 partitioning choice immediately provides the linearized CCSDT wavefunction ansatz. It is here that we choose to only include linear singles and doubles in the ground state, resulting in

which is the LCCSD similarity transformed Hamiltonian. As with LCCD, the matrix elements of Eq. (65) are identical to Eq. (37) differing only by the $T ˆ$ amplitudes.

To assess both the quality and general behavior of these EOM approximations when computing vertical singlet excitation energies, we employ two established collections of gas phase molecules (see Figure 1): the set of 24 organic molecules from Schreiber et al.^37–39 (referred to hereafter as the Mülheim set) and the set of 11 organic molecules from Caricato et al.^2,40 (referred to here as the Yale set). This Mülheim set contains 121 reference singlet excitation energies³⁷ which use molecular geometries obtained with an MP2/6-31G* optimization in which the first row 1s atomic orbitals are dropped. The extent and balanced nature of the Mülheim set readily lends itself to obtaining statistics between different theoretical methods. To compare the approximate EOM-CC methods proposed here, we take our best theoretical reference to be EOM-CCSDT-3,^41,56 which is known to give both accurate results and is a systematic EOM-CC type theory which will provide a balanced comparison. The Yale set choice in molecules is a subset of the Mülheim set with the addition of acetaldehyde but considers a different selection of excitation energies (69 in all) whose reference values are taken from gas phase experiments (the experiment literature references can be found in Ref. 2). This comparison to experiment provides a different and valuable perspective with which to augment theory vs. theory examinations. Geometries for the Yale set are obtained from an MP2/6-311+G** optimization.

Calculations using the Mülheim set typically³⁷ use the TZVP basis set (3s1p on hydrogen 5s3p1d on first row atoms) from Schäfer et al.⁴² This is a light weight basis set that allows the use of many different kinds of computational methods with only a modest cost in overall accuracy.³⁹ While perfectly acceptable when comparing different methods to each other, this basis is not accurate enough to compare with experimental results. This is especially true when calculating Rydberg states, where the inclusion of many diffuse functions is key.^43,44 When computing the excitation energies for the Yale set, we use the doubly diffuse function augmented Dunning basis,^45,46 d-aug-cc-pVDZ (4s3p on hydrogen and 5s4p3d on first row atoms). This basis is comparable to the 6-311(3+,3+)G** basis used by Caricato et al. in their calculations^2,40,47 with almost no gain in including a third set of diffuse functions.⁴⁸ 

All of the reported electronic structure results were obtained using the serial ACESII⁴⁹ and parallel Aces4⁵⁰ ab initio quantum chemistry packages. Calculations were performed on the University of Florida HiPerGator high performance cluster. Throughout this paper we use RMS to mean “root mean square,” MAD to mean “mean absolute deviation,” SD to mean “standard deviation,” and MD to mean “absolute maximum deviation.”

Developed two decades ago for canonical Hartree-Fock²⁷ EOM-MBPT(2) is a well established perturbation method which has been largely ignored by the community in favor of other approximations such as CC2²⁴ or ADC(2).¹⁹ The EOM-MBPT(2) method remains relevant because of the non-iterative O(n⁵) computational cost of the ground-state wavefunction and the general accuracy of MBPT(2). It is also reasonably accurate for equilibrium vertical excitation energies. This accuracy, within the Mülheim test set, is illustrated in the correlation Figures 2(a), 3(a), and Table V. Compared to the standard EOM-CCSD and EOM-CCSDT-3 methods, EOM-MBPT(2) has an MAD of 0.13 eV and 0.31 eV, respectively, consistent with the known ∼0.2 eV triples correction to EOM-CCSD.

The EOM-MBPT(2) approximation and its Löwdin partitioned variant^29,51 was recently benchmarked using the Yale test set (using the 6-311(3+,3+)G** basis set) against EOM-CCSD and experiment by Goings et al.⁴⁷ Our computed EOM-MBPT(2)/d-aug-cc-pVDZ values, presented in Tables I–IV, yield nearly identical results. The EOM-MBPT(2) method is completely satisfactory when compared to experimental Rydberg states (such as with the alkenes and carbonyl’s) but is much less reliable for the n → π^∗ and π → π^∗ heterocycle valence states. The comparative difficulty in describing valence states is a limitation of EOM theory in general which is only satisfactorily overcome by using methods that include triples or alternative Hamiltonians such as similarity transformed equation-of-motion coupled-cluster theory^52–55 (STEOM-CC).

By employing the second-order similarity transformed Hamiltonians, Eqs. (59) and (65), we now have an equation-of-motion theory completely consistent with the LCCD (LCCSD) ground state wavefunction. The linear EOM methods have the same iterative O(n⁶) computational scaling as canonical EOM-CCSD. However, with the removal of the quadratic amplitude intermediates, the computational overhead and I/O demands are significantly reduced. Benchmark data from the Mülheim set show a systematic overestimate of excitation energies for both EOM-LCCD and EOM-LCCSD, with a MAD relative to EOM-CCSD of 0.25 and 0.35 eV, respectively (see the correlation Figures 2, 3, and Table V). The comparison of the linear EOM methods against EOM-CCSDT-3 is less satisfactory with a MAD of 0.45 and 0.56 eV for EOM-LCCD and EOM-LCCSD, respectively.

With an acceptable SD of 0.11 and 0.10 eV compared to EOM-CCSD (see Table IV), the excitation energies and relative state orderings obtained with the linear EOM methods are reliably consistent compared to canonical EOM-CCSD. This is further illustrated in Tables I–III where the systematic and consistent overestimation of the excitation energies is balanced by a nearly exact agreement in state ordering and relative state energy differences with a MAD of 0.04, 0.04, and 0.05 eV for EOM-MBPT(2), LCCD, and LCCSD, respectively (see the bottom of Table IV). For approximate methods, getting these relative energies correct is just as useful as quantitatively accurate excitation energies.

We have expanded the coupled-cluster similarity transformed Hamiltonian using general coupled-cluster perturbation theory to obtain the arbitrary order coupled-cluster perturbation theory effective Hamiltonian, Eq. (19). The inclusive (the sum of all orders 0 to n) form of the Hamiltonian is inserted into the standard electronic-excitation equation-of-motion theory to give the completely general EOM-CCPT (Eq. (27)). The result is a way of using the standard equation-of-motion theory to directly compute excitation energies that are consistent with a given CCPT ground state wavefunction.

By choosing the generalized many-body perturbation theory partitioning of the Hamiltonian, we re-derive the well known EOM-MBPT(2)²⁹ (EOM-CCSD(2)²⁷) equations. Two alternative choices in Hamiltonian partitioning, the particle rank conserving and mixing schemes, are also explored. These two schemes, referred to as PR0 and PR1, generate infinite-order amplitudes that correspond to the linearized coupled-cluster doubles (CCD) and CCSD wavefunction. Using the inclusive similarity transformed Hamiltonian derived from the PR0 and PR1 partitioning schemes provides an equation-of-motion method perturbatively consistent with a LCCD and LCCSD reference wavefunction.

These approximate equation-of-motion methods are benchmarked by employing the Mülheim^37–39 and Yale^2,40 small organic molecule data sets. We use the diverse Mülheim test set to benchmark our new methods against the canonical EOM-CCSD and EOM-CCSDT-3 methods, while the Yale test set is employed in benchmarking against experimental spectra. Our methods are found to consistently over estimate excitation energies, relative to EOM-CCSD and EOM-CCSDT-3 by ∼0.25 and ∼0.5 eV, respectively. The precise statistics and correlation plots can be found in Table V and Figures 2 and 3. The very good ∼0.1 eV standard deviation compared with EOM-CCSD suggests that the predicted relative spectra will be much more accurate than the absolute excitation energies. This is supported by studying the benchmark values from the Yale set given in Tables I–IV where the presented approximate EOM methods obtain relative state orderings and energies to within ∼0.04 eV of canonical EOM-CCSD. The systematic overestimation compared to the complete EOM-CCSD suggests that the similarity transformed Hamiltonian (and thus the excited state spectrum) based on the linear CC wavefunction is over-approximated compared to the ground-state description leading to an increased separation. Having explored three different CC ground state wavefunctions based on a linear amplitude approximation (and associated linear amplitude similarity transformed Hamiltonian), it is evident that nonlinear amplitude contributions play a significantly larger role in the excited state spectrum than for the ground state energy. Comparison with the CC2 method, which has all product $T ˆ 1$ and $T ˆ 1 T ˆ 2$ amplitudes but no quadratic $T ˆ 2$ terms, is further instructive as the average CC2 excitation energy is between³⁸ EOM-CCSD and CC3. The ideal second order singles and doubles limited EOM approximation would then apparently lay somewhere between CC2 and linear EOM-CC.