The cluster perturbation series, CPS(D), for coupled cluster singles and doubles excitation energies is considered. It is demonstrated that the second-order model CPS(D-2) is identical to the configuration interaction singles with perturbative doubles, CIS(D) model. The third-order model, CPS(D-3), provides excitation energies of coupled cluster singles and doubles (CCSD) quality in the sense that the difference between CPS(D-3) and CCSD excitation energies is of the same size or smaller than the effect of adding triples corrections to CCSD excitation energies. We further show that the third-order corrections can be efficiently implemented, in particular, when the resolution of the identity approximation is used for integrals. We also show that the CPS(D-3) excitation energies can be determined for system sizes that are far beyond what can be considered in conventional CCSD excitation energy calculations.

## I. INTRODUCTION

During the last decades, coupled cluster (CC) theory has evolved to become the wave-function method of choice for describing single-configuration dominated molecular systems.^{1,2} The coupled cluster singles and doubles (CCSD) model is a robust and reliable wave-function model that gives quantitatively correct results.^{3} It is widely used to calculate both energies and molecular properties. In particular, the calculation of CCSD excitation energies has turned out successfully. CCSD excitation energies are determined as eigenvalues of the CCSD response eigenvalue equation.^{4} Solving the CCSD response eigenvalue equation requires iterative algorithms, which in each iteration have an *N*^{6} computational scaling, where *N* denotes the system size. During the last decades, it has been an active research field to reduce the iterative *N*^{6} computational scaling of CCSD excitation energy calculations without compromising the quality of the excitation energies.

In this paper, we consider a series of excitation energy corrections in orders of the fluctuation potential that added to coupled cluster singles (CCS) excitation energies converges to CCSD excitation energies. The series is a member of the cluster perturbation (CP) series that were developed in Paper II^{5} for excitation energies, where excitation energy corrections are determined that added to the excitation energies of a CC parent state formally converge to the excitation energies of a CC target state. In the context of Paper II,^{5} the CP series that targets the CCSD excitation energies is labeled CPS(D) since excitation energy corrections are added to the excitation energies of a coupled cluster singles (CCS) parent state. In Paper II,^{5} the local convergence was investigated for CP excitation energy series for general CC parent and CC target states, whereas the global convergence of the CP excitation energy series will be examined in a forthcoming paper. In this paper, the convergence of the CPS(D) series for excitation energies is investigated from a cost-benefit perspective.

The CP perturbation series for excitation energies were developed in Paper II.^{5} The derivation of the CPS(D) excitation energy series can be simplified compared to the general derivation presented in Paper II^{5} and performed without reference to the general perturbation framework of CP theory. In Sec. II B, we present this simplified derivation of the CPS(D) series for excitation energies. In Sec. II D, we compare the computational cost for calculating the low-order corrections with the cost of a traditional calculation of CCSD excitation energies, where the CCSD cluster amplitude equation and the CCSD response eigenvalue equation have to be solved explicitly.

CP excitation energy corrections are determined without the need for solving CC target state amplitude and response equations explicitly as in standard CC theory. This makes the evaluation of low-order corrections well suited for an efficient implementation. The third-order correction is in particular well suited for an efficient and massively parallel implementation when the resolution of the identity (RI) approximation for integrals^{6,7} is used, since the third-order correction only contains products of first-order amplitudes and two-electron integrals. In Sec. II E, we describe how the CPS(D-2) and CPS(D-3) excitation energy corrections can be efficiently implemented using the RI approximation for integrals, and we show in Sec. III C that our implementation can be applied to systems of sizes that are far beyond what can be considered in a standard CCSD excitation energy calculation.

In the CPS(D) series, we calculate approximations to CCSD excitation energies by adding perturbation corrections to CCS excitation energies. For the series to be well behaved, it is important that the targeted CCSD excited state is well described at the CCS level. We denote such excitations as single-replacement dominated excitations and describe in Sec. II C how they may be identified in practice.

In Secs. III A and III B, we perform excitation energy calculations on a large variety of molecular systems, which show that the CPS(D-3) model gives excitation energies of CCSD quality. We also compare the accuracy of CPS(D-3) and CCSD excitation energies relative to calculations where the effect of triples is considered and show that very similar performance is obtained for the CPS(D-3) and CCSD model.

The development of response function theory for a CC wave-function was initiated by Monkhorst.^{8} Koch and Jørgensen^{9} later derived explicit expressions for the CC response eigenvalue equations and Koch *et al.*^{4} reported the first CCSD excitation energy calculations in 1990. We note that Symmetry Adapted Cluster Configuration Interaction (SAC-CI) theory is closely related to the CC response eigenvalue problem.^{10} A large variety of approximate models has been developed, where approximations have been introduced in the CCSD Jacobian. Stanton and Bartlett have termed the excitation energies that are determined from the CCSD response eigenvalue equation, equation-of-motion coupled cluster singles and doubles (EOM-CCSD) excitation energies,^{11} and Stanton and Gauss have introduced a second-order approximation to the CCSD response eigenvalue equation and termed the model the EOM-CCSD(2) model.^{12} In a singles-and-doubles partitioned form, the EOM-CCSD(2) eigenvalue equation has been termed the EOM-MBPT2 model (where MBPT stands for many-body perturbation theory).^{13} The EOM-CCSD(2) and EOM-MBPT2 models have an iterative *N*^{6} computational scaling, but with a smaller pre-factor than the CCSD response eigenvalue equation. Introducing a diagonal approximation, containing orbital energy differences in the doubles-doubles block of the Jacobian, leads to the P-EOM-MBPT2 model,^{13} which has an iterative *N*^{5} computational scaling.

Excitation energies in a singles-and-doubles excitation framework may also be determined using the CC2 model.^{14} By introducing the RI approximation for the integrals, Hättig and Weigend have presented a very efficient CC2 implementation for application on large molecular systems.^{15} Excitation energies in a singles-and-doubles excitation framework have also been determined using the second-order polarization propagator approximation (SOPPA) where the excitation energies are also determined using an iterative *N*^{5} scaling algorithm.^{16–18}

In the algebraic-diagrammatic construction (ADC) method, excitation energy corrections are determined in orders of the fluctuation potential and the excitation energies formally converge to the full configuration interaction (FCI) excitation energies.^{19} This series has been implemented through third-order. In zeroth-order, orbital energy differences between the excited and the HF configurations are obtained and in first-order, configuration interaction singles (CIS) excitation energies are determined. In both the second- and the third-order ADC methods—ADC(2)^{19–21} and ADC(3)^{19,22–24}—excitation energies are determined inside a singles-and-doubles excitation framework. In ADC(2), all second-order contributions to the excitation energies are taken into account and the method has an iterative *N*^{5} scaling. In the ADC(3) method, all third-order excitation energy contributions are also added and the method has an iterative *N*^{6} scaling. Since the ADC(3) method is derived inside a singles-and-doubles excitation framework, it is derived inside the same excitation framework as the CCSD Jacobian eigenvalue equation. The ADC(3) scheme may therefore be viewed as a third-order approximation to the CCSD response eigenvalue equation, similarly as the ADC(2) scheme may be viewed as a second-order approximation to the CCSD response eigenvalue equation, with the only difference that the ADC method is developed within a linearly parametrized framework for the time evolution of the ground state, whereas the CC method is developed within an exponentially parametrized framework.

The methods described above have all been developed for calculating the total spectrum of excitation energies by solving eigenvalue equations. Alternatively, excitation energy corrections may be determined one at a time thereby avoiding the need to solve eigenvalue equations explicitly. This is the strategy employed in the CIS(D) model of Head-Gordon^{25} as well as other related models.^{26,27} It is also the strategy employed in this paper for the CPS(D) excitation energy series. The CPS(D-2) model is equal to the CIS(D) model and may be viewed as the equivalent of the iterative *N*^{5} scaling methods, described above, where excitation energies are calculated one by one. In the same way the, CPS(D-3) model may be viewed as the equivalent of the ADC(3) model, where the excitation energies are calculated one by one. The reason for the complex structure of the ADC(3) excitation energy eigenvalue equation compared to the CPS(D-3) model is that the perturbation framework that is used in ADC theory has orbital energy differences as zeroth-order excitation energies, whereas in CP theory the zeroth-order excitation energies are obtained from the CCS model.

A detailed investigation of the performance of the excitation energies of the EOM-MBPT2, P-EOM-MBPT2, CIS(D), and CC2 models compared to CCSD excitation energies has been performed by Goings *et al.*^{28} Extensive benchmarking of the CIS(D), CC2 and EOM-CCSD(2) models has been performed and revealed a serious deficiency of the CIS(D) and the CC2 model for describing Rydberg states.^{29,30} We report a similar benchmark investigation including the CPS(D-3) model in Sec. III B.

## II. THEORY

### A. Excitation energies from the CCSD response eigenvalue equation

CCSD excitation energies may be determined as eigenvalues of the non-Hermitian CCSD Jacobian,^{9}

where **R**_{x} and **L**_{x} are right and left eigenvectors, respectively, for the excited state *x*, and the eigenvalue *ω*_{x} is the excitation energy. The left and right eigenvectors have been chosen to be biorthonormal. The CCSD Jacobian matrix elements are given by^{31}

where *H*_{0} is the electronic Hamiltonian and *T* is the CCSD cluster operator,

The cluster operator contains the cluster amplitudes $t\mu i$, and the many-body excitation operators $\theta \mu i$ perform excitations from occupied to unoccupied orbitals in the Hartree–Fock (HF) reference configuration,

with *i* denoting the excitation level and *μ*_{i} a given excitation at this level. The cluster amplitudes satisfy the cluster amplitude equations,

In a standard CCSD excitation energy calculation, Eq. (5) is first solved to obtain the CCSD amplitudes and the excitation energies are then determined by solving the CCSD response eigenvalue equation in Eq. (1a). Both the cluster amplitude equation and the response eigenvalue equation may be solved using iterative algorithms that have an *N*^{6} computational scaling with the size of the system, *N*.

### B. Simple derivation of the CPS(D) excitation energy series

In this section, we describe a simple derivation of the CPS(D) series for excitation energies. In Sec. II B 1, we derive a perturbation series in orders of the fluctuation potential for the CCSD Jacobian, where the zeroth-order Jacobian is chosen to satisfy the CCS Jacobian eigenvalue equation, and in Sec. II B 2 the Jacobian series is used to determine a series of corrections to a CCS excitation energy that target a CCSD excitation energy.

#### 1. Perturbation series for the CCSD Jacobian

We will now determine a perturbation series for the CCSD Jacobian in orders of the fluctuation potential,

The zeroth-order Jacobian **J**^{(0)} satisfies

where the zeroth-order excitation energy is required to be a CCS excitation energy,

where $\omega xCCS$ satisfies the CCS Jacobian eigenvalue equation,

The matrix elements of the CCS Jacobian are given by^{31}

and since the CCS Jacobian is Hermitian, we also have

In order to identify the Jacobian expansion in Eq. (6), we carry out a Baker–Campbell–Hausdorff (BCH) expansion of the CCSD Jacobian in Eq. (2) and introduce a Møller–Plesset partitioning of the electronic Hamiltonian,^{32,33}

where *f* is the Fock operator and Φ the fluctuation potential. We assume canonical HF orbitals, such that the Fock operator matrix is diagonal. We then obtain

where we have used that contributions from the Fock operator vanish for terms containing more than one commutator and that the expansion terminates after the fourth term due to rank reduction.^{32} The second, third, and fourth terms in Eq. (13) are of at least second-order in Φ and the zeroth- and first-order terms in the Jacobian expansion in Eq. (6) therefore have to be determined from the first term in Eq. (13). The first term is denoted the extended CCS Jacobian and may be partitioned as

where

The singles subblock of **J**^{(0)} is equal to the CCS Jacobian **J**^{CCS} in Eq. (10) and contains both a Fock operator and a fluctuation potential operator. We have thus defined a generalized zeroth-order Jacobian which includes contributions from the fluctuation potential. This is not a conventional choice, however it makes the zeroth-order solutions more accurate and the perturbation corrections smaller. Using **J**^{(0)} in Eq. (15), the Jacobian eigenvalue problem in Eq. (7a) can be written as

where *ε*_{D} is a diagonal matrix containing orbital energy differences in the double excitations space,

and where **J**^{CCS} satisfies the CCS Jacobian eigenvalue equation in Eq. (9). We therefore have that Eq. (8) is satisfied and that the zeroth-order response vectors in Eq. (7) become

Note that we do not consider zeroth-order solutions of the type $\omega x(0)=\epsilon D$, since these solutions describe double-replacement dominated excitations which can only be described properly using models that contain triple excitations.

To identify **J**^{(0)} and **J**^{(1)} in Eq. (21), we have imposed that the CCS Jacobian in the singles subblock of **J**^{(0)} becomes a zeroth-order term and that the singles subblock of **J**^{(1)} vanishes [see Eqs. (15) and (16)]. We have thereby removed direct relaxation effects in the singles subspace through first-order in the Jacobian of Eq. (6). To identify the second- and higher-order terms of the CCSD Jacobian expansion in Eq. (6) from Eq. (21), we will also require that the cluster amplitudes in the *T* operator in Eq. (21) are determined from a perturbation series in which direct relaxation effects in the singles subspace are removed, such that direct relaxation effects are fully removed from the Jacobian expansion. To accomplish this, we perform a BCH expansion of Eq. (5),

where we have used the Brillouin theorem.^{34} We now insert Eq. (14) in the second term of Eq. (22), giving the cluster amplitude equation,

From Eq. (23), the *k*th-order amplitude equations may be written in a two-component form as

where $(.)(k)$ denotes that the terms of order *k* in Φ are picked up and gathered from the expression in parenthesis. Equation (24) confirms that direct relaxation in the single excitation subspace is removed in the cluster amplitude series, since the first term on the right-hand side of Eq. (24a) only introduces a coupling between the single excitations space and the double excitations space, and since sets of linear equations containing the zeroth-order CCS Jacobian are solved in the singles excitation space.

Using the cluster amplitudes from Eq. (24), the Jacobian series in Eq. (6) may be determined from Eq. (21). The cluster amplitudes and Jacobian matrix elements required to calculate excitation energies up to third-order are given in Table I. For higher order expressions, see Tables X and XI.

$t\mu 1(1)$ | = | 0 |

$\u2211\nu 1J\mu 1\nu 1CCSt\nu 1(2)$ | = | $\u2212\u2009\mu 1|\Phi ,T2(1)|HF\u2009$ |

$\epsilon \mu 2t\mu 2(1)$ | = | −⟨ μ_{2}|Φ|HF⟩ |

$\epsilon \mu 2t\mu 2(2)$ | = | $\u2212\u2009\mu 2|\Phi ,T2(1)|HF\u2009$ |

$J\mu i\nu j(0)$ | = | $J\mu i\nu jCCS\delta i1\delta j1+\epsilon \nu j\delta \mu i\nu j\delta i2\delta j2$ |

$J\mu i\nu j(1)$ | = | $\u2009\mu i|\Phi ,\theta \nu j|HF\u2009\delta i1\delta j2+\u2009\mu i|\Phi ,\theta \nu j|HF\u2009\delta i2\delta j1$ |

$\u2009+\u2009\mu i|\Phi ,\theta \nu j|HF\u2009\delta i2\delta j2$ | ||

$J\mu i\nu j(2)$ | = | $\u2009\mu i|\Phi ,T2(1),\theta \nu j|HF\u2009$ |

$J\mu i\nu j(3)$ | = | $\u2009\mu i|\Phi ,T(2),\theta \nu j|HF\u2009$ |

$t\mu 1(1)$ | = | 0 |

$\u2211\nu 1J\mu 1\nu 1CCSt\nu 1(2)$ | = | $\u2212\u2009\mu 1|\Phi ,T2(1)|HF\u2009$ |

$\epsilon \mu 2t\mu 2(1)$ | = | −⟨ μ_{2}|Φ|HF⟩ |

$\epsilon \mu 2t\mu 2(2)$ | = | $\u2212\u2009\mu 2|\Phi ,T2(1)|HF\u2009$ |

$J\mu i\nu j(0)$ | = | $J\mu i\nu jCCS\delta i1\delta j1+\epsilon \nu j\delta \mu i\nu j\delta i2\delta j2$ |

$J\mu i\nu j(1)$ | = | $\u2009\mu i|\Phi ,\theta \nu j|HF\u2009\delta i1\delta j2+\u2009\mu i|\Phi ,\theta \nu j|HF\u2009\delta i2\delta j1$ |

$\u2009+\u2009\mu i|\Phi ,\theta \nu j|HF\u2009\delta i2\delta j2$ | ||

$J\mu i\nu j(2)$ | = | $\u2009\mu i|\Phi ,T2(1),\theta \nu j|HF\u2009$ |

$J\mu i\nu j(3)$ | = | $\u2009\mu i|\Phi ,T(2),\theta \nu j|HF\u2009$ |

#### 2. Arbitrary-order excitation energy corrections

The excitation energy *ω*_{x} and right excitation vector **R**_{x} in Eq. (1a) may be expanded in orders of the fluctuation potential,

where the zeroth-order excitation energy is equal to $\omega xCCS$ and where the zeroth-order right excitation vector is given in Eq. (19). We will assume that **R**_{x} is intermediate normalized against $Lx(0)=Rx(0)\u2020$, i.e.,

Eq. (27) implies that

To determine the *k*th-order correction to a CCS excitation energy targeting a CCSD excitation energy, we substitute Eqs. (6), (25), and (26) in Eq. (1a) and collect terms of order *k*, giving

which may be rearranged as

where we have used Eq. (8) and where **I** is an identity matrix. Projecting Eq. (31) against the zeroth-order left eigenvector, $Lx(0)$, we obtain the *k*th-order correction to the CCS excitation energy,

where we have used Eqs. (7b), (7c), and (29). Note that the right-hand side of Eq. (32) depends only on right eigenvectors through order (*k* − 1).

When excitation energy corrections have been determined through order *k*, they may be substituted in Eq. (31) to determine the *k*th-order correction to the right eigenvector. In the two-component form, the *k*th-order right eigenvalue equation may be written as

where to obtain Eq. (33a) we have used that the first-order Jacobian does not have a singles-singles block [see Eq. (16)], and to obtain Eq. (33b) we have used that the doubles component of the zeroth-order right eigenvector vanishes [cf. Eq. (17)]. To obtain Eq. (33b), we have further used that the matrix $\epsilon D\u2212\omega xCCSI$ is non-singular and diagonal. The singles component of the *k*th-order right eigenvector is obtained solving the linear equation in Eq. (33a).

Explicit expressions for the calculation of the CPS(D) excitation energies up to third-order are given in Table II, where the following notation is used:

For explicit expressions of the CPS(D) excitation energies up to sixth-order, see Tables XII–XIV. Note that direct relaxation effects are removed from the single excitations subspace when the *k*th-order response amplitudes are determined, because the *k*th-order singles response amplitudes $RxS(k)$ do not depend on $RxS(k\u22121)$ but only on $RxD(k\u22121)$ (see Eq. (33a) and Tables XIII and XIV), and because $RxS(k)$ is determined solving sets of linear equations containing the zeroth-order CCS Jacobian. For comparison, the *k*th-order singles cluster amplitudes $tS(k)$ do not depend on $tS(k\u22121)$ but only on $tD(k\u22121)$, and $tS(k)$ is determined solving sets of linear equations containing the zeroth-order CCS Jacobian [see Eq. (24a) and Table X].

$Rx\nu 1(1)$ | = | 0 |

$\u2211\nu 1J\mu 1\nu 1CCS\u2212\omega xCCS\u2009\delta \mu 1\nu 1Rx\nu 1(2)$ | = | $\omega x(2)Rx\mu 1CCS\u2212\u2009\mu 1|\Phi ,T2(1),RxCCS|HF\u2009\u2212\u2009\mu 1|\Phi ,Rx2(1)|HF\u2009$ |

$\epsilon \mu 2\u2212\omega xCCSRx\mu 2(1)$ | = | $\u2212\u2009\mu 2|\Phi ,RxCCS|HF\u2009$ |

$\epsilon \mu 2\u2212\omega xCCSRx\mu 2(2)$ | = | $\u2212\u2009\mu 2|\Phi ,T2(1),RxCCS|HF\u2009\u2212\u2009\mu 2|\Phi ,Rx2(1)|HF\u2009$ |

$\omega x(1)$ | = | 0 |

$\omega x(2)$ | = | $\u2009LxCCS|\Phi ,T2(1),RxCCS|HF\u2009+\u2009LxCCS|\Phi ,Rx2(1)|HF\u2009$ |

$\omega x(3)$ | = | $\u2009LxCCS|\Phi ,T(2),RxCCS|HF\u2009+\u2009LxCCS|\Phi ,Rx2(2)|HF\u2009$ |

$Rx\nu 1(1)$ | = | 0 |

$\u2211\nu 1J\mu 1\nu 1CCS\u2212\omega xCCS\u2009\delta \mu 1\nu 1Rx\nu 1(2)$ | = | $\omega x(2)Rx\mu 1CCS\u2212\u2009\mu 1|\Phi ,T2(1),RxCCS|HF\u2009\u2212\u2009\mu 1|\Phi ,Rx2(1)|HF\u2009$ |

$\epsilon \mu 2\u2212\omega xCCSRx\mu 2(1)$ | = | $\u2212\u2009\mu 2|\Phi ,RxCCS|HF\u2009$ |

$\epsilon \mu 2\u2212\omega xCCSRx\mu 2(2)$ | = | $\u2212\u2009\mu 2|\Phi ,T2(1),RxCCS|HF\u2009\u2212\u2009\mu 2|\Phi ,Rx2(1)|HF\u2009$ |

$\omega x(1)$ | = | 0 |

$\omega x(2)$ | = | $\u2009LxCCS|\Phi ,T2(1),RxCCS|HF\u2009+\u2009LxCCS|\Phi ,Rx2(1)|HF\u2009$ |

$\omega x(3)$ | = | $\u2009LxCCS|\Phi ,T(2),RxCCS|HF\u2009+\u2009LxCCS|\Phi ,Rx2(2)|HF\u2009$ |

### C. Single-replacement dominated excitation energies

In the CPS(D) series, we calculate approximations to CCSD excitation energies by adding perturbative corrections to CCS excitation energies. For the series to be well behaved and for the results to be trustworthy, it is important that the targeted CCSD state is well described at the CCS level. We refer to such excitations as being single-replacement dominated.

For an excitation to be single-replacement dominated, it is vital that the zeroth-order CCS state is the dominant excitation component in the expansion of **R**_{x} in Eq. (26). Distinguishing between the singles and doubles terms, **R**_{x} may be written as

with

To obtain Eq. (38) we have used that $RxD(0)$ vanishes [see Eq. (19)], and to obtain Eq. (39) we have used Eqs. (9c) and (11). For a single-replacement dominated excitation, the norm of the doubles component |**R**_{xD}| truncated at a given order must be small compared to the norm of the singles component and thus much smaller than $|RxS(0)|=1$. One requirement for an excitation to be single-replacement dominated is therefore that

In addition, for a single-replacement dominated excitation we must require that the correction to the CCS state is much smaller than the CCS state itself, i.e.,

also when the right-hand side of Eq. (41) is truncated at a given order. The sum of the norms $\u2211k=2|RxS(k)|$ truncated at a given order may thus be used to judge if the CPS(D) series converges and how fast it converges. To conclude, the norms of the higher-order singles and doubles components of the response amplitudes are important for estimating the quality of the excitation energy corrections in the CPS(D) series.

In CP theory, like in CC theory, the excitation energies are size-intensive.^{4} Furthermore, excitations in a large molecular system are in general local and with a small interaction between excitations that are located at large distances from each other. When solving the Jacobian eigenvalue equation, the size of the excitation vector therefore does not grow with the system size unlike the size of the ground state vector when a configuration interaction (CI) eigenvalue equation is solved. The size of the single and double excitations component in Eqs. (37) and (38) can therefore be used with confidence also for large molecular systems to test for single-replacement dominance.

### D. Comparison of the computational scaling for a conventional and a perturbative approach

We will now compare the leading-order computational scaling of a CCSD excitation energy calculation using the conventional approach described in Sec. II A with the one using the CPS(D) series, in which excitation energy perturbation corrections are determined through sixth-order. In the conventional approach, we solve the cluster amplitude equations in Eq. (5) and the CCSD right response eigenvalue equation in Eq. (1a). The solution to both equations may be determined using iterative algorithms that require about the same number of iterations *n*_{it}, where *n*_{it} is in the range 10–15.^{35–38} In each iteration of the cluster amplitude equation, we have to provide a trial solution vector,

where **b**_{S} and **b**_{D} are the singles and the doubles components of the trial solution vector, respectively. The singles component of the amplitude equation contains terms with a leading-order computational scaling *N*^{5}, where *N* refers to the size of the molecular system. The leading-order computational scaling for the full CCSD calculation is determined by the doubles component which contains *N*^{6} scaling terms. The doubles component may be expressed as

where

and

is a Hamiltonian with *B*_{1}-transformed integrals.^{32} The construction of $H0B1$ has a leading-order computational scaling *N*^{5}. The first term in Eq. (43), $\u2009\mu 2|H0B1|HF\u2009$, also has a leading-order computational scaling *N*^{5}. The computational scaling for setting up the doubles component is determined by the last two terms in Eq. (43), which both have an *N*^{6} leading-order scaling. We denote the scaling for these two terms as $LD1$ and $LD2$, respectively,

Following the derivation of the CCSD equations in Chap. 13 of Ref. 32, we can express the $LD1$ and $LD2$ sixth-order scaling in terms of the number of occupied (*O*) and virtual (*V*) orbitals,

where we have assumed *B*_{1}-transformed integrals. In more advanced implementations the pre-factor for the dominating scaling contributions can change but the leading-order scaling will in all cases be $LD1\u2243V4O2$ and $LD2\u2243V3O3$. The computational scaling for solving amplitude equations is thus

When solving the CCSD Jacobian eigenvalue equation in Eq. (1a), the leading-order scaling of each iteration is determined by setting up the linear transformation of the CCSD Jacobian on a trial solution, **J b**. The doubles component of this linear transformation determines the leading-order computational scaling and may be written as

The first term is a matrix element between a HF state and a double excited state of a two-electron Hamiltonian,

with some modified integrals. Setting up $H0T1B1$ and evaluating the matrix element has a leading-order computational scaling *N*^{5}. The computational scaling for the rest of the terms in Eq. (52) may be written as

where we have used the structural similarity between the terms in Eqs. (54)–(56) and the terms in Eqs. (47) and (48). The leading-order computational scaling for solving the CCSD Jacobian eigenvalue equation is therefore

For comparison, let us now express the scaling of the lowest order corrections in the CPS(D) series in terms of $LD1$ and $LD2$. The construction of $\omega x(k)$ requires ground state amplitudes, **t**^{(k−1)}, and the right eigenvector, $Rx(k\u22121)$. Explicit expressions for these vectors through fifth-order can be found in Tables X, XIII, and XIV. In addition, the excitation energy corrections involve a projection onto the singles left CCS eigenvector, $LxCCS$. Explicit expressions for $\omega x(k)$ through sixth-order are given in Table XII.

In Table III, we have tabulated the computational scaling for constructing the $\omega x(k)$ correction for *k* = 1, 2, …, 6 as well as the leading-order scaling for the construction of the corresponding **t**^{(k−1)} and $Rx(k\u22121)$ vectors. The leading-order scaling for the vectors is determined by their doubles component and it has been obtained following the same strategy as used previously in this subsection.

Order k
. | t^{(k)}
. | $Rx(k)$ . | $\omega x(k)$ . |
---|---|---|---|

1 | … | … | … |

2 | $1LD1$ | $2LD1$ | |

3 | $1LD1+1LD2$ | $2LD1+1LD2$ | $3LD1$ |

4 | $2LD1+1LD2$ | $5LD1+3LD2$ | $3LD1+2LD2$ |

5 | $3LD1+2LD2$ | $8LD1+5LD2$ | $7LD1+4LD2$ |

6 | $11LD1+7LD2$ |

Order k
. | t^{(k)}
. | $Rx(k)$ . | $\omega x(k)$ . |
---|---|---|---|

1 | … | … | … |

2 | $1LD1$ | $2LD1$ | |

3 | $1LD1+1LD2$ | $2LD1+1LD2$ | $3LD1$ |

4 | $2LD1+1LD2$ | $5LD1+3LD2$ | $3LD1+2LD2$ |

5 | $3LD1+2LD2$ | $8LD1+5LD2$ | $7LD1+4LD2$ |

6 | $11LD1+7LD2$ |

### E. Closed-shell expressions and implementation of CPS(D) excitation energy corrections through third-order

In this section, we will determine the second- and third-order corrections to CCS singlet excitation energies for closed-shell systems. From Table II, we see that the second-order correction requires the determination of the doubles component of the first-order cluster and response amplitudes ($t\mu 2(1)$ and $R\mu 2(1)$), while for the third-order correction, the singles and doubles components of the second-order cluster amplitudes ($t\mu 1(2)$ and $t\mu 2(2)$), as well as the doubles components of the second-order response amplitudes ($R\mu 2(2)$), are needed. We will first determine expressions for the required amplitudes and then use those to obtain working equations for the second- and third-order excitation energy corrections. In this section, we use the notation presented in Table IV.

Occupied canonical MOs: | i, j, k, l |

Virtual canonical MOs: | a, b, c, d |

Canonical MOs of unspecified occupancy: | p, q, r, s |

Auxiliary basis functions (for the RI approximation): | P, Q |

Atomic orbitals (AOs): | α, β, γ, δ |

Fock matrix: | F_{pq} = δ_{pq}ϵ_{p} |

Two-electron integrals: | (pq|rs) = ∑_{αβγδ}C_{αp}C_{βq}C_{γr}C_{δs}(αβ|γδ) |

Occupied canonical MOs: | i, j, k, l |

Virtual canonical MOs: | a, b, c, d |

Canonical MOs of unspecified occupancy: | p, q, r, s |

Auxiliary basis functions (for the RI approximation): | P, Q |

Atomic orbitals (AOs): | α, β, γ, δ |

Fock matrix: | F_{pq} = δ_{pq}ϵ_{p} |

Two-electron integrals: | (pq|rs) = ∑_{αβγδ}C_{αp}C_{βq}C_{γr}C_{δs}(αβ|γδ) |

Response amplitudes for singles and doubles calculations are conveniently obtained using the singlet biorthonormal basis,^{32}

where *E*_{pq} is a singlet excitation operator and

with

For *ai* ≥ *bj* and *ck* ≥ *dl*, we have

In addition, we have

A detailed derivation of the singlet basis presented above can be found in Chap. 13 of Ref. 32.

For the biorthonormal basis in Eq. (58), a double excitations operator for the cluster amplitudes may be expressed as

where *τ*_{aibj} denotes the non-redundant doubles cluster amplitudes in the biorthonormal basis. For the derivation of the working equations, it is convenient to expand the amplitudes to a symmetrized form defined for all (*ai*, *bj*) pairs,

and the double excitations operator then becomes

#### 1. Cluster and response amplitudes

To illustrate the use of the biorthonormal basis and double excitations operators in the derivation of the closed-shell equations, we consider the case of the first-order doubles cluster amplitudes in more details. The derivation presented here follows the strategy presented in Chap. 13 of Ref. 32 for the derivation of the CCSD closed-shell amplitude equations. Let us start with the first-order amplitude equation in the singlet biorthonormal basis (see Table I),

and introduce the symmetrized form of the double excitations cluster operator from Eq. (67),

Proceeding with the strategy established in Ref. 32 for the evaluation of the matrix elements, we get

The non-redundant amplitudes can then be obtained from Eq. (66)

Following the same strategy, the first-order response amplitudes in Table II can be expressed in the singlet basis and in the symmetrized form,

where the *barred* integrals are calculated as

Note that here and in the following, we are suppressing the *x* subscript used to denote a given excitation. Closed-shell expressions for the second-order doubles cluster and response amplitudes can be obtained in the same way,

$Xaibjt=Xaibjt.A+Xaibjt.B+Xaibjt.C+Xaibjt.D+Xaibjt.E$ |

$Xaibjt.A=\u2212\u2211dktaidk(1)(bd|kj)$ |

$Xaibjt.B=\u2212\u2211dktakdj(1)(bd|ki)$ |

$Xaibjt.C=12\u2211dctdjci(1)(ca|bd)$ |

$Xaibjt.D=12\u2211kltakbl(1)(ik|jl)$ |

$Xaibjt.E=\u2211ckt\u0303aick(1)(ck|bj)$ |

$Xaibjt=Xaibjt.A+Xaibjt.B+Xaibjt.C+Xaibjt.D+Xaibjt.E$ |

$Xaibjt.A=\u2212\u2211dktaidk(1)(bd|kj)$ |

$Xaibjt.B=\u2212\u2211dktakdj(1)(bd|ki)$ |

$Xaibjt.C=12\u2211dctdjci(1)(ca|bd)$ |

$Xaibjt.D=12\u2211kltakbl(1)(ik|jl)$ |

$Xaibjt.E=\u2211ckt\u0303aick(1)(ck|bj)$ |

$XaibjR.1$ | $XaibjR.2$ |

$XaibjR.1.A=\u2212\u2211dkRaidk(1)(bd|kj)$ | $XaibjR.2.A=\u2212\u2211dktaidk(1)(b\xafd|jk)+(bd|j\xafk)$ |

$XaibjR.1.B=\u2212\u2211dkRakdj(1)(bd|ki)$ | $XaibjR.2.B=\u2212\u2211dktakdj(1)(b\xafd|ik)+(bd|\u012bk)$ |

$XaibjR.1.C=12\u2211dcRdjci(1)(ca|bd)$ | $XaibjR.2.C=\u2211dctdjci(1)(ca|b\xafd)$ |

$XaibjR.1.D=12\u2211klRakbl(1)(ik|jl)$ | $XaibjR.2.D=\u2211kltakbl(1)(ik|j\xafl)$ |

$XaibjR.1.E=\u2211ckR\u0303aick(1)(ck|bj)$ | $XaibjR.2.E=\u2211ckt\u0303aick(1)(ck|b\xafj)+(ck|bj\xaf)$ |

$XaibjR.2.F=\u2211ctaicj(1)F\xafbc\u2032\u2212\u2211ktaibk(1)F\xafkj\u2032$ | |

$F\xafbc\u2032=\u2211dl2(bc|ld)\u2212(bd|lc)RdlCCS$ | $F\xafkj\u2032=\u2211dl2(kj|ld)\u2212(kd|lj)RdlCCS$ |

$XaibjR.1$ | $XaibjR.2$ |

$XaibjR.1.A=\u2212\u2211dkRaidk(1)(bd|kj)$ | $XaibjR.2.A=\u2212\u2211dktaidk(1)(b\xafd|jk)+(bd|j\xafk)$ |

$XaibjR.1.B=\u2212\u2211dkRakdj(1)(bd|ki)$ | $XaibjR.2.B=\u2212\u2211dktakdj(1)(b\xafd|ik)+(bd|\u012bk)$ |

$XaibjR.1.C=12\u2211dcRdjci(1)(ca|bd)$ | $XaibjR.2.C=\u2211dctdjci(1)(ca|b\xafd)$ |

$XaibjR.1.D=12\u2211klRakbl(1)(ik|jl)$ | $XaibjR.2.D=\u2211kltakbl(1)(ik|j\xafl)$ |

$XaibjR.1.E=\u2211ckR\u0303aick(1)(ck|bj)$ | $XaibjR.2.E=\u2211ckt\u0303aick(1)(ck|b\xafj)+(ck|bj\xaf)$ |

$XaibjR.2.F=\u2211ctaicj(1)F\xafbc\u2032\u2212\u2211ktaibk(1)F\xafkj\u2032$ | |

$F\xafbc\u2032=\u2211dl2(bc|ld)\u2212(bd|lc)RdlCCS$ | $F\xafkj\u2032=\u2211dl2(kj|ld)\u2212(kd|lj)RdlCCS$ |

Finally, the second-order singles ground state amplitudes $tai(2)$ given in Table I can be obtained by solving a set of linear equations for single excitations,

where

and where the linear transformation of a trial vector $bbjCCS$ with the CCS Jacobian is given by

#### 2. Second-order excitation energy correction

The second-order correction to a CCS excitation energy is given in Table II. In the singlet basis, the second-order correction reads

Using the derivation strategy detailed in the preceding part of this subsection (see also Ref. 32), the following expression is obtained:

where *L*_{pqrs} = 2(*pq*|*rs*) − (*ps*|*rq*) and

After some rearrangement, we obtain the following working equations:

with

As discussed in Sec. II C, the norm of higher-order singles and doubles response amplitudes can be used to judge the quality of the CPS(D) series. For the singles space, we consider the norm of the second-order singles response amplitudes which can be obtained from a set of linear equations involving the CCS Jacobian and intermediate quantities already used in the calculation of the CPS(D-2) excitation energy correction,

#### 3. Third-order excitation energy correction

The third-order correction to a CCS excitation energy is given in Table II. In the singlet basis, the third-order correction reads

and the following expression is obtained:

After some rearrangement, we obtain the final working equation

where

and

#### 4. CPS(D-3) implementation

The closed-shell CPS(D-3) equations for excitation energies have been implemented in the LSDALTON^{39,40} program, in the spirit of the RI-CC2 algorithm by Hättig and Weigend in Ref. 15. The two-electron integrals are calculated using the RI approximation as

The formal structure of the second- and third-order excitation energy corrections in Eqs. (85) and (94) allows us to calculate all four-index quantities (amplitudes and integrals) in batches which can be contracted on-the-fly, eliminating the need to store four-index quantities in core-memory or on disk. When the RI approximation is used for the integrals, we may also avoid the recalculation of two-electron integrals since they can be computed in advance and stored as the 3-index fitting coefficients $BpqP$.

A schematic representation of the algorithm used to calculate the third-order correction is given in Algorithm 1. The two contributions to the third-order correction in Eq. (94) are calculated in two different sections which loop over batches of virtual indices. The two loops are built in the same way and the most expensive step in Algorithm 1 is the calculation of integrals of the type (*ai*|*bj*) which scales as *O*^{2}*V*^{2}*N*_{α} when using the RI approximation (see steps 7 and 15). *N*_{α} denotes the number of auxiliary functions, whereas *O* and *V* denote the number of occupied and virtual orbitals, respectively. However, the calculation of the intermediates *X*^{t} and *X*^{R} (steps 6 and 14 in Algorithm 1) involves *N*^{6} scaling steps, as can been seen from Algorithm 2. Indeed, the contraction of doubles amplitudes with fully virtual two-electron integrals (step 11 in Algorithm 2) scales as *O*^{2}*V*^{4}, which represents the leading-order scaling term in the $LD1$ transformation introduced in Sec. II D. As shown in Table III, we do need one $LD1$ transformation to obtain the *X*^{t} intermediate, which comes from the second-order cluster amplitudes, and two $LD1$ transformations for the *X*^{R} intermediate, which arise in the second-order response amplitudes. This should be compared with a standard CCSD calculation of excitation energies in which one and two $LD1$ transformations are required in each iteration of the cluster amplitude and the response eigenvalue equations, respectively (see Sec. II D).

1: Get transformation matrices, $X\xaf$, $Y\xaf$, $X\u0306$, $Y\u0306$. |

2: Get integrals fitting coefficients, $BaiP$, $BijP$, $B\u012daP$, $Bi\u0103P$. |

3: Get transformed Fock matrices, $F\xafia$, $F\xafij\u2032$, $F\xafab\u2032$. |

4: Calculate $BabP$, and $B\u0101bP$ in batches over b and store on disk. |

5: for A (batch over virtual index) do |

6: Get $XAibjt$ (see Algorithm 2). |

7: $(iA|\u0306jb)=\u2211PB\u012dAPBjbP+Bi\u0102PBjbP+BiAPBj\u0306bP+BiAPBjb\u0306P$. |

8: Get $t\u0303Aibj(2)=(2XAibjt\u2212XAjbit)/(\u03f5i\u2212\u03f5a+\u03f5j\u2212\u03f5b)$ |

9: $\omega (3)=\omega (3)+\u2211a\u2208A\u2211ibjt\u0303aibj(2)P^ijabRaiCCSF\xafjb\u2212(ia|\u0306jb)$ |

10: end for A |

11: Discard $B\u012daP$ and $Bi\u0103P$. |

12: Get integrals fitting coefficients, $B\u0101iP$, $Ba\u012bP$, $B\u012bjP$. |

13: for A (batch over virtual index) do |

14: Get $XAibjR=XAibjR.1+XAibjR.2$ Algorithm (see 2). |

15: $(Ai|\xafbj)=\u2211PB\u0100iPBbjP+BA\u012bPBbjP+BAiPBb\xafjP+BAiPBbj\xafP$. |

16: Get $R\u0303Aibj(2)=(2XAibjR\u2212XAjbiR)/(\u03f5i\u2212\u03f5a+\u03f5j\u2212\u03f5b+\omega CCS)$ |

17: $\omega (3)=\omega (3)+\u2211a\u2208A\u2211ibjR\u0303aibj(2)+2RaiCCStbj(2)\u2212RajCCStbi(2)(ai|\xafbj)$ |

18: end for A |

1: Get transformation matrices, $X\xaf$, $Y\xaf$, $X\u0306$, $Y\u0306$. |

2: Get integrals fitting coefficients, $BaiP$, $BijP$, $B\u012daP$, $Bi\u0103P$. |

3: Get transformed Fock matrices, $F\xafia$, $F\xafij\u2032$, $F\xafab\u2032$. |

4: Calculate $BabP$, and $B\u0101bP$ in batches over b and store on disk. |

5: for A (batch over virtual index) do |

6: Get $XAibjt$ (see Algorithm 2). |

7: $(iA|\u0306jb)=\u2211PB\u012dAPBjbP+Bi\u0102PBjbP+BiAPBj\u0306bP+BiAPBjb\u0306P$. |

8: Get $t\u0303Aibj(2)=(2XAibjt\u2212XAjbit)/(\u03f5i\u2212\u03f5a+\u03f5j\u2212\u03f5b)$ |

9: $\omega (3)=\omega (3)+\u2211a\u2208A\u2211ibjt\u0303aibj(2)P^ijabRaiCCSF\xafjb\u2212(ia|\u0306jb)$ |

10: end for A |

11: Discard $B\u012daP$ and $Bi\u0103P$. |

12: Get integrals fitting coefficients, $B\u0101iP$, $Ba\u012bP$, $B\u012bjP$. |

13: for A (batch over virtual index) do |

14: Get $XAibjR=XAibjR.1+XAibjR.2$ Algorithm (see 2). |

15: $(Ai|\xafbj)=\u2211PB\u0100iPBbjP+BA\u012bPBbjP+BAiPBb\xafjP+BAiPBbj\xafP$. |

16: Get $R\u0303Aibj(2)=(2XAibjR\u2212XAjbiR)/(\u03f5i\u2212\u03f5a+\u03f5j\u2212\u03f5b+\omega CCS)$ |

17: $\omega (3)=\omega (3)+\u2211a\u2208A\u2211ibjR\u0303aibj(2)+2RaiCCStbj(2)\u2212RajCCStbi(2)(ai|\xafbj)$ |

18: end for A |

1: Calculate $tAidk(1)$ |

2: Read $BcAP$ from disk |

3: for D (batch over virtual index) do |

4: Extract $tAiDk(1)$ from $tAidk(1)$ |

5: Read $BbDP$ from disk |

6: $(bD|kj)=\u2211PBbDPBkjP$ |

7: $XAibjt=XAibjt\u2212\u2211DktAiDk(1)(bD|kj)$ |

8: $XAibjt=XAibjt\u2212\u2211DktAkDj(1)(bD|ki)$ |

9: $(cA|bD)=\u2211PBcAPBbDP$ |

10: Calculate $tDjci(1)$ |

11: $XAibjt=XAibjt+12\u2211DctDjci(1)(cA|bD)$ |

12: end for D |

13: for L (batch over occupied index) do |

14: Extract $tAkbL(1)$ from $tAkbl(1)$ |

15: $(ik|jL)=\u2211PBikPBjLP$ |

16: $XAibjt=XAibjt+12\u2211kLtAkbL(1)(ik|jL)$ |

17: end for L |

18: $t\u0303Aick(1)=2tAick(1)\u2212tAkci(1)$ |

19: $YAiP=\u2211ckt\u0303Aick(1)BckP$ |

20: $XAibjt=XAibjt+\u2211PYAiPBbjP$ |

1: Calculate $tAidk(1)$ |

2: Read $BcAP$ from disk |

3: for D (batch over virtual index) do |

4: Extract $tAiDk(1)$ from $tAidk(1)$ |

5: Read $BbDP$ from disk |

6: $(bD|kj)=\u2211PBbDPBkjP$ |

7: $XAibjt=XAibjt\u2212\u2211DktAiDk(1)(bD|kj)$ |

8: $XAibjt=XAibjt\u2212\u2211DktAkDj(1)(bD|ki)$ |

9: $(cA|bD)=\u2211PBcAPBbDP$ |

10: Calculate $tDjci(1)$ |

11: $XAibjt=XAibjt+12\u2211DctDjci(1)(cA|bD)$ |

12: end for D |

13: for L (batch over occupied index) do |

14: Extract $tAkbL(1)$ from $tAkbl(1)$ |

15: $(ik|jL)=\u2211PBikPBjLP$ |

16: $XAibjt=XAibjt+12\u2211kLtAkbL(1)(ik|jL)$ |

17: end for L |

18: $t\u0303Aick(1)=2tAick(1)\u2212tAkci(1)$ |

19: $YAiP=\u2211ckt\u0303Aick(1)BckP$ |

20: $XAibjt=XAibjt+\u2211PYAiPBbjP$ |

Regarding the memory requirements, we see from Algorithm 1 that for typical calculations the largest quantities that are kept in memory are the three-index fitting coefficients of the type $BaiP$ and of size *N*_{α}*OV*. The storage of the four-index quantities is avoided by batching over one virtual index. Similarly, to avoid keeping in memory the fully virtual fitting coefficients $BabP$ of size *N*_{α}*V*^{2}, we construct them in batches over *b* and store them on disk (step 4 in Algorithm 1). The virtual fitting coefficients are then read in batches from disk (steps 2 and 5 in Algorithm 2) and contracted on-the-fly.

We note that the calculation of the third-order correction, as described in Algorithms 1 and 2, may be straightforwardly parallelized by distributing the loops in steps 5 and 13 of Algorithm 1 over several nodes. Each node would then calculate a subset of fitting coefficients as well as independent contributions to the third-order corrections. The fully virtual fitting coefficients could then be stored on disk, as it is done in the current version, or distributed over the nodes’ memory. In the latter case, this would imply node-to-node communication of batches of the $BbdP$ fitting coefficients in step 5 of Algorithm 2.

To summarize, the calculation of the CPS(D-3) excitation energies starts by solving the CCS Jacobian eigenvalue equation in Eq. (9a) using an iterative algorithm. The solution of the CCS problem formally scales as *N*^{4} with the system size and can be applied to large molecular systems.^{41} We then proceed with the calculation of the second-order correction from the CPS(D-2) model. The $\eta \xaf(2)$ vector in Eq. (86) is first calculated using a modified version of the RI-CC2 algorithm presented in Ref. 15. The second-order correction is then directly obtained by contracting $\eta \xaf(2)$ with the CCS excitation vector, as in Eq. (85). The $\eta \xaf(2)$ vector is also used to calculate the right-hand side of Eq. (91) which is solved iteratively to provide the $RS(2)$ vector needed to calculate the $|RS(2)|$ norm for diagnostic purposes. We note that while the calculation of the $\eta \xaf(2)$ vector is an *N*^{5} scaling operation, the solution of Eq. (91) using iterative algorithm requires only *N*^{4} scaling operations. Another set of linear equation has to be solved to obtain the second-order ground state singles vectors $t1(2)$ in Eq. (78). Similarly, the right-hand side (*η*^{(2)}) is computed with a modified RI-CC2 algorithm, while the iterative solution involves only CCS-like transformations. Finally, the third-order corrections in Eq. (94) are calculated using Algorithm 1, where all quantities are either singles quantities calculated in advance, or doubles quantities that can be computed and contracted on-the-fly (see Algorithm 2).

## III. RESULTS

We now report calculation of excitation energies for the CPS(D) series. The local convergence of the CPS(D) series has been investigated in Paper II.^{5} In this paper we look at the convergence of the series from a cost-benefit point of view.

### A. Convergence of the CPS(D) series for excitation energies

In this section, the convergence of the CPS(D) series is examined from a computational cost perspective where the evaluation of lower order CPS(D) excitation energy corrections are compared with a conventional evaluation of a CCSD excitation energy. In Sec. II C we have shown that the CPS(D) series of excitation energy corrections can only be applied with confidence for excitations that are single-replacement dominated and we discussed how the norms of the singles and doubles excitation vectors can be used for identifying excitations that are single-replacement dominated. In this section, we show how the norms of the singles excitation vectors may be used as a practical diagnostic tool for identifying when an excitation is single-replacement dominated. When CCSD excitation energies are calculated, the effects of triples and higher excitations are neglected. Therefore, excitation energies in the CPS(D) series are said to have CCSD quality when the deviation of a CPS(D) excitation energy from the CCSD excitation energy is of the same size as the effects of triple excitations. We show in this section that for single-replacement dominated excitations, the CPS(D-3) model gives excitation energies of CCSD quality.

As a prototype example of the convergence of the CPS(D) series of excitation energy corrections, we report in Table VII the deviation from the CCSD excitation energy *ω*_{k} − *ω*^{CCSD} for three low lying triplet excitation energies of CO, where

The excitation energy corrections *ω*^{(k)} for orders *k* = 0, 1, …, 10 are also reported together with the excitation vector norms $|RS(k)|$ and$|RD(k)|$. The supplementary material to this paper contains the corresponding tables for three low lying singlet states of CO, for three low lying singlet and triplet states of H_{2}O and for three low lying triplet–triplet transitions of O_{2}. The calculations on O_{2}, CO and H_{2}O were carried out at the geometries: *R*_{CO} = 1.129 Å, *R*_{OO} = 1.208 Å, *R*_{OH} = 0.957 Å, *α*_{HOH} = 104.2° and all used the aug-cc-pVDZ basis.^{42,43} For O_{2} we used an unrestricted HF reference. The CPS(D) corrections through an arbitrary order have been implemented using a Python interface to the PSI4 program (Psithon)^{44} and employing the Numerical Python Library (Numpy).^{45} This implementation should not be confused with the efficient RI-based implementation described in Sec. II E 4 which is restricted to the CPS(D-2) and CPS(D-3) models for singlet transitions of closed-shell system. The CCSD and CC3 excitation energy calculations were performed with the DALTON^{39,46} (for CO and H_{2}O) and PSI4^{44} (for O_{2}) programs.

. | ^{3}Π
. | ^{3}Σ^{+}
. | ^{3}Δ
. | |||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

k
. | ω_{k} − ω^{CCSD}^{a}
. | ω^{(k)}
. | $|RS(k)|$ . | $|RD(k)|$ . | ω_{k} − ω^{CCSD}^{b}
. | ω^{(k)}
. | $|RS(k)|$ . | $|RD(k)|$ . | ω_{k} − ω^{CCSD}^{c}
. | ω^{(k)}
. | $|RS(k)|$ . | $|RD(k)|$ . |

0 | −0.5514 | 5.8593 | 1.0000 | 0.0000 | −0.6489 | 7.7728 | 1.0000 | 0.0000 | −0.6589 | 8.7240 | 1.0000 | 0.0000 |

1 | −0.5514 | 0.0000 | 0.0000 | 0.2061 | −0.6489 | 0.0000 | 0.0000 | 0.2870 | −0.6589 | 0.0000 | 0.0000 | 0.2553 |

2 | 0.1311 | 0.6824 | 0.0488 | 0.0469 | 0.2533 | 0.9022 | 0.0639 | 0.0489 | 0.1599 | 0.8188 | 0.0717 | 0.0559 |

3 | 0.0275 | −0.1036 | 0.0152 | 0.0346 | −0.0102 | −0.2635 | 0.0140 | 0.0232 | −0.0372 | −0.1971 | 0.0140 | 0.0258 |

4 | 0.0632 | 0.0357 | 0.0123 | 0.0161 | 0.0255 | 0.0357 | 0.0081 | 0.0127 | 0.0098 | 0.0470 | 0.0067 | 0.0123 |

5 | −0.0047 | −0.0679 | 0.0116 | 0.0106 | −0.0119 | −0.0374 | 0.0059 | 0.0061 | −0.0119 | −0.0217 | 0.0060 | 0.0047 |

6 | 0.0128 | 0.0175 | 0.0069 | 0.0076 | 0.0087 | 0.0207 | 0.0024 | 0.0039 | 0.0059 | 0.0178 | 0.0031 | 0.0029 |

7 | −0.0050 | −0.0178 | 0.0069 | 0.0050 | −0.0036 | −0.0123 | 0.0021 | 0.0020 | −0.0034 | −0.0093 | 0.0019 | 0.0017 |

8 | 0.0042 | 0.0092 | 0.0047 | 0.0041 | 0.0041 | 0.0077 | 0.0011 | 0.0012 | 0.0027 | 0.0061 | 0.0015 | 0.0010 |

9 | −0.0015 | −0.0057 | 0.0041 | 0.0030 | −0.0017 | −0.0058 | 0.0011 | 0.0008 | −0.0017 | −0.0044 | 0.0010 | 0.0007 |

10 | 0.0016 | 0.0030 | 0.0031 | 0.0024 | 0.0018 | 0.0036 | 0.0006 | 0.0007 | 0.0013 | 0.0030 | 0.0007 | 0.0005 |

ω^{CC3} − ω^{CCSD} | −0.0686^{d} | 0.1125^{e} | 0.0779^{f} |

. | ^{3}Π
. | ^{3}Σ^{+}
. | ^{3}Δ
. | |||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

k
. | ω_{k} − ω^{CCSD}^{a}
. | ω^{(k)}
. | $|RS(k)|$ . | $|RD(k)|$ . | ω_{k} − ω^{CCSD}^{b}
. | ω^{(k)}
. | $|RS(k)|$ . | $|RD(k)|$ . | ω_{k} − ω^{CCSD}^{c}
. | ω^{(k)}
. | $|RS(k)|$ . | $|RD(k)|$ . |

0 | −0.5514 | 5.8593 | 1.0000 | 0.0000 | −0.6489 | 7.7728 | 1.0000 | 0.0000 | −0.6589 | 8.7240 | 1.0000 | 0.0000 |

1 | −0.5514 | 0.0000 | 0.0000 | 0.2061 | −0.6489 | 0.0000 | 0.0000 | 0.2870 | −0.6589 | 0.0000 | 0.0000 | 0.2553 |

2 | 0.1311 | 0.6824 | 0.0488 | 0.0469 | 0.2533 | 0.9022 | 0.0639 | 0.0489 | 0.1599 | 0.8188 | 0.0717 | 0.0559 |

3 | 0.0275 | −0.1036 | 0.0152 | 0.0346 | −0.0102 | −0.2635 | 0.0140 | 0.0232 | −0.0372 | −0.1971 | 0.0140 | 0.0258 |

4 | 0.0632 | 0.0357 | 0.0123 | 0.0161 | 0.0255 | 0.0357 | 0.0081 | 0.0127 | 0.0098 | 0.0470 | 0.0067 | 0.0123 |

5 | −0.0047 | −0.0679 | 0.0116 | 0.0106 | −0.0119 | −0.0374 | 0.0059 | 0.0061 | −0.0119 | −0.0217 | 0.0060 | 0.0047 |

6 | 0.0128 | 0.0175 | 0.0069 | 0.0076 | 0.0087 | 0.0207 | 0.0024 | 0.0039 | 0.0059 | 0.0178 | 0.0031 | 0.0029 |

7 | −0.0050 | −0.0178 | 0.0069 | 0.0050 | −0.0036 | −0.0123 | 0.0021 | 0.0020 | −0.0034 | −0.0093 | 0.0019 | 0.0017 |

8 | 0.0042 | 0.0092 | 0.0047 | 0.0041 | 0.0041 | 0.0077 | 0.0011 | 0.0012 | 0.0027 | 0.0061 | 0.0015 | 0.0010 |

9 | −0.0015 | −0.0057 | 0.0041 | 0.0030 | −0.0017 | −0.0058 | 0.0011 | 0.0008 | −0.0017 | −0.0044 | 0.0010 | 0.0007 |

10 | 0.0016 | 0.0030 | 0.0031 | 0.0024 | 0.0018 | 0.0036 | 0.0006 | 0.0007 | 0.0013 | 0.0030 | 0.0007 | 0.0005 |

ω^{CC3} − ω^{CCSD} | −0.0686^{d} | 0.1125^{e} | 0.0779^{f} |

^{a}

*ω*^{CCSD} = 6.4107 eV.

^{b}

*ω*^{CCSD} = 8.4217 eV.

^{c}

*ω*^{CCSD} = 9.3830 eV.

^{d}

*ω*^{CC3} = 6.3421 eV.

^{e}

*ω*^{CC3} = 8.5342 eV.

^{f}

*ω*^{CC3} = 9.4609 eV.

We first consider the local convergence of the excitation energy series. From Table VII and the tables in the supplementary material, we see that the second- and third-order corrections most often have opposite sign and that the third-order corrections typically are about one third of the numerical value of the second-order correction. Exceptions occur for example, for the A^{1}Π state of CO, where the second-order correction is −0.1398 eV and the third-order correction is −0.1515 eV. In this case, the third-order correction has the same sign as the second-order correction and is numerically larger than the second-order correction. Also for the B ^{1}Σ^{+} state of CO, the second- and third-order corrections have the same sign. For the ^{3}Π_{g} state of O_{2}, the second- and third-order corrections are 0.5896 eV and −0.4217 eV respectively, and are thus numerically large and similar in magnitude.

The numerical value of the second-order excitation energy corrections covers a large interval with the largest value being the second-order excitation energy correction of 2.0104 eV for the C ^{3}Δ_{u} state of O_{2} and with the smallest value being −0.1398 eV for the A ^{1}Π state of CO. For the third-order corrections, the largest corrections are of the size ∼0.5 eV and are obtained for several excitations. For example for the C ^{3}Δ_{u} and A ^{3}$\Sigma u+$ states of O_{2}. The smallest third-order correction is −0.1036 eV and is obtained for the ^{3}Π state of CO. The largest second-order excitation energy error compared to CCSD is −0.4963 eV for the ^{1}A_{2} state of H_{2}O. The largest third-order error is 0.1511 eV for the B ^{1}Σ^{+} state of CO. For this state $|RS(2)|$ has a value 0.56 which is much larger than the $|RS(2)|$ values for the other considered states where 0.31 is the largest value that is obtained for the ^{3}A_{2} state of H_{2}O. For all states except the B ^{1}Σ^{+} state of CO the third-order excitation energy errors are all smaller than 0.1 eV. In general there is a tendency that small third-order errors are obtained for small $|RS(2)|$ values.

The doubles first-order excitation vector norms $|RD(1)|$ cover an interval of 0.2-0.4. The doubles second-order excitation vector norms $|RD(2)|$ are typically one third to one fourth of the first-order correction. The magnitude of the doubles excitation norms does not appear to have any direct influence on the excitation energy errors for the considered interval of doubles excitation norms.

The fourth-order excitation energy corrections are numerically smaller than the third-order excitation energy corrections and in general of opposite sign. The largest fourth-order excitation energy correction is −0.3079 eV for the ^{3}A_{2} state of H_{2}O giving the largest fourth-order excitation energy error of −0.2256 eV. The fourth- and fifth-order excitation energy corrections generally have opposite sign. The largest fifth-order excitation energy correction is 0.2794 eV for the ^{3}A_{2} state of H_{2}O giving one of the largest fifth-order errors of 0.0538 eV. The fifth-order errors are in general smaller than the third-order errors and are typically in the range 0.05–0.01 eV. In sixth-order the errors in the excitation energies are typically of the order of 0.01 eV or smaller.

The effect of excitations higher than doubles is usually dominated by the effect of triples which can be well represented by Δ*ω*^{CC3} = *ω*^{CC3} − *ω*^{CCSD} since CC3 excitation energies accurately describe full CCSDT excitation energies.^{29,47,48} In Table VII (and the tables in the supplementary material), we report the effect of triples as Δ*ω*^{CC3}. From the tables we see that the third-order excitation energy errors for all excitations are smaller and in most cases much smaller than Δ*ω*^{CC3}. The third-order excitation energies can therefore be considered to have CCSD quality. Further, for most of the excitation energies the third-order excitation energies are closer than the CCSD excitation energies to the CC3 excitation energy.

We now consider the global convergence of the CPS(D) excitation energy series. For fifth- and higher-orders the largest excitation energy errors are found for the B ^{1}Σ^{+} state of CO. The singles excitation norms for this state indicate that the excitation energy series for this state diverge rather than converge. The $|RS(2)|$ norm is 0.56 and the higher-order norms $|RS(k)|,\u2009k=2,3,\u2026,10$ fluctuate between values in the range 0.06–0.33 with no sign of convergence. Furthermore, the corrections show no definite converging trend. Calculations at higher-order also show that the CPS(D) series diverge for the B ^{1}Σ^{+} state. For the other states the $|RS(2)|$ norm are much smaller with the largest $|RS(2)|$ norm being 0.31 for the ^{3}A_{2} state of H_{2}O. For these states the singles and doubles higher-order norms typically decrease for increasing orders, and the excitation energies are converging toward the CCSD excitation energies. The $|RS(2)|$ norm of 0.56 for the B ^{1}Σ^{+} state also indicate that this excitation is not single-replacement dominated. Even though the deviation from the CCSD excitation energies for the B ^{1}Σ^{+} state is only 0.15 eV at the CPS(D-3) level it is still significantly larger than 0.10 eV that is the largest deviation from the CCSD excitation energies for the other states in the tables. The other excitations are clearly single-replacement dominated since the largest $|RS(2)|$ norm is 0.31 for these states.

The leading order scaling for the excitation energy corrections in the CPS(D-k) series is *N*^{6} (for k > 2). It corresponds to the $LD1$ and $LD2$ terms given in Sec. II D. Using Eqs. (49) and (50) and Table III we see that the leading order scalings for *ω*^{(k)}, *k* = 3, 4, 5, 6 are 3 *V*^{4}*O*^{2}, 3*V*^{4}*O*^{2}, 7*V*^{4}*O*^{2}, and 11*V*^{4}*O*^{2}. From the analysis in Sec. II D, the leading order scaling for evaluating a CCSD excitation energy becomes 3*n*_{it}*V*^{4}*O*^{2} where *n*_{it} is the number of iterations for solving the cluster amplitude and the response eigenvalue equation. In Fig. 1 we have plotted excitation energy errors for the CPS(D-k) series for *k* = 0, 1, …, 6 for the C ^{3}Δ_{u} state of O_{2}. For each order we have in parentheses given the leading order computational effort relative to the leading order effort for a conventional CCSD excitation energy calculation assuming *n*_{it} = 15. The numbers in parentheses are obtained using that a conventional CCSD excitation energy calculation has a leading order scaling of 3*n*_{it}*V*^{4}*O*^{2} = 45*V*^{4}*O*^{2}. For the third-order correction the leading order scaling is 3*V*^{4}*O*^{2} and the relative effort for a third-order calculation is thus (3/45) as displayed in Fig. 1. The relative effort for orders 4, 5, and 6 are obtained in the same way. In second-order, the leading order scaling is *N*^{5} and the computational effort relative to a CCSD calculation is thus vanishing. For that reason, we have marked the relative effort for the second-order correction as (≪1/1000). In zeroth-order, the CCS excitation energies may be calculated using a linear-scaling algorithm and we have marked the corresponding relative effort as (≪1/10 000) in Fig. 1.

Considering the fast local convergence of the CPS(D) series for excitation energies it thus becomes computationally very tractable to use low-order corrections in the CPS(D) series to obtain approximate CCSD excitation energies. The second-order correction which is equal to the CIS(D) correction gives a fast estimate of the doubles correction to CCS excitation energies, while the third-order correction gives excitation energies of CCSD quality. For excitation energies where $|RS(2)|$ is smaller than 0.3, the CP excitation energy series converges, and the largest error in the third-order excitation energies compared to CCSD is smaller than 0.1 eV. Further, the third-order correction can be calculated at a small fraction of the cost of a conventional CCSD excitation energy calculation. In Sec. III B we will further benchmark the quality of the excitation energies obtained using the CPS(D-2) and CPS(D-3) models.

### B. Accuracy and comparison of the CPS(D-2), CPS(D-3), and CC2 models

The calculations presented in Sec. III A indicate that the CPS(D-3) model is a very good compromise in terms of accuracy and computational cost. In this section we thus focus on the CPS(D-3) model and perform a statistical analysis to further test the accuracy of the second- and third-order excitation energies compared to CCSD excitation energies and also compared to CC3 excitation energies. For comparison, CC2 excitation energies are also reported. We also consider and test the norm of the second-order singles excitation vector $|RS(2)|$, as a diagnostics for identifying single-replacement dominated excitations. For the statistical analysis, we have used the test set of Goings *et al.*,^{28} which was first used to investigate different low-scaling approximations to CCSD excitation energies. This test set contains a total of 69 excited states, 30 of which can be characterized as valence states, while the remaining 39 states are Rydberg states. In order to enable a direct comparison to the results of Goings *et al.*, we have used the 6-311(3+,3+)G^{**} basis set^{49} where all atoms have been augmented with three diffuse functions.

The CCS, CPS(D-2), and CPS(D-3) excitation energies have been calculated as described in Sec. II E, using the LSDALTON program.^{39,40} The CC2 excitation energies are directly obtained from Ref. 28, while the CCSD and CC3 numbers have been obtained with the DALTON program. In order to compare the results between the different models we have followed the assignment of the CCS and CC2 states given in Ref. 28. By construction the CPS(D-2) and CPS(D-3) states are then assigned to the parent CCS states. It is however necessary to associate a CCSD target state to all the considered CCS parent states. For each of the calculated CCSD states we have calculated the overlap with the CCS states and made the assignment based on the largest overlap. The mapping between the CCSD and CC3 states is straightforward. In the supplementary material we report excitation energies for all the states and models considered here. The assignment and character of each transition is also reported.

Let us start by comparing the CC2, CPS(D-2), and CPS(D-3) models as approximations to the CCSD model for excitation energies. In Fig. 2 we have plotted the CPS(D-3) errors against the CPS(D-2) errors taking CCSD excitation energies as a reference. The three plots in Fig. 2 correspond to the full set of states as well as subsets of states with $|RS(2)|<1.0$ and $|RS(2)|<0.3$. Similarly, we report the CPS(D-3) errors against the CC2 errors in Fig. 3 and the CC2 errors against the CPS(D-2) errors in Fig. 4. From Figs. 2 and 3, we see that the CPS(D-3) models outperforms the second-order models which perform particularly poorly for Rydberg states. Both figures also reveal that cancellation of errors can happen and that the second-order models occasionally perform better than the CPS(D-3) model (points above the diagonal). However, those cases are rare and can often be identified since they are generally associated with a large value of the $|RS(2)|$ diagnostic which indicates that these excitation energies cannot be trusted. Indeed, we see that the $|RS(2)|$ measure is a valuable tool since it identifies the largest errors from the CPS(D-3) results. In Fig. 4 we see that the CPS(D-2) and CC2 models have a relatively similar accuracy with respect to the reference CCSD excitation energies. Both models perform particularly well for valence states and gives significantly larger errors for Rydberg states. The $|RS(2)|$ diagnostic may also be used for the CPS(D-2) and CC2 results but since the errors for these models are much larger than the errors of the CPS(D-3) model, all the large errors can in general not be scanned away using the $|RS(2)|$ diagnostic.

In Table VIII, we present a statistical analysis of the errors in the excitation energies of the CC2, CCS, CPS(D-2), and CPS(D-3) models compared to the reference CCSD numbers. The results are also divided into the contributions from valence and Rydberg states and subsets of states with $|RS(2)|<1.0$ and $|RS(2)|<0.3$. From Table VIII, we see that the CCS model performs very poorly with an average absolute error of 1.0 eV for all states compared to the CCSD results and that no significant difference is observed between valence and Rydberg transitions. However, when keeping only singles dominated transitions based on the $|RS(2)|$ diagnostic ($|RS(2)|<1.0$ and $|RS(2)|<0.3$) the CCS average absolute error drops to 0.84 eV and 0.79 eV, respectively. The second-order models (CC2 and CPS(D-2)) perform particularly well for valence states as already noted in Ref. 28. For example, for the CPS(D-2) model we have $\Delta \xafabs=0.15$ eV for valence states only, while for Rydberg states it increases to $\Delta \xafabs=0.45$ eV. The results are similar for the CC2 model. Finally, Table VIII confirms that the new CPS(D-3) model performs significantly better than the second-order models, especially for Rydberg states which results in a much more balanced description of the two types of transitions ($\Delta \xafabs=0.09$ eV and $\Delta \xafabs=0.10$ eV for valence and Rydberg states, respectively). Here the impact of the $RS(2)$ diagnostic measure can be seen from the maximum error, dropping from 0.52 eV for all states to 0.14 eV for states with $|RS(2)|<0.3$.

All states . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . |
---|---|---|---|---|

$\Delta \xafabs$ | 0.37 | 1.00 | 0.32 | 0.09 |

Δ_{max} | 0.82 | 2.64 | 0.96 | 0.52 |

$\Delta \xaf$ | −0.36 | 0.85 | −0.26 | −0.02 |

Δ_{std} | 0.28 | 0.86 | 0.33 | 0.13 |

All states . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . |
---|---|---|---|---|

$\Delta \xafabs$ | 0.37 | 1.00 | 0.32 | 0.09 |

Δ_{max} | 0.82 | 2.64 | 0.96 | 0.52 |

$\Delta \xaf$ | −0.36 | 0.85 | −0.26 | −0.02 |

Δ_{std} | 0.28 | 0.86 | 0.33 | 0.13 |

Valence states . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . |
---|---|---|---|---|

$\Delta \xafabs$ | 0.23 | 1.03 | 0.15 | 0.09 |

Δ_{max} | 0.80 | 2.64 | 0.67 | 0.41 |

$\Delta \xaf$ | −0.21 | 0.92 | −0.02 | −0.01 |

Δ_{std} | 0.21 | 0.76 | 0.20 | 0.13 |

Valence states . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . |
---|---|---|---|---|

$\Delta \xafabs$ | 0.23 | 1.03 | 0.15 | 0.09 |

Δ_{max} | 0.80 | 2.64 | 0.67 | 0.41 |

$\Delta \xaf$ | −0.21 | 0.92 | −0.02 | −0.01 |

Δ_{std} | 0.21 | 0.76 | 0.20 | 0.13 |

Rydberg states . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . |
---|---|---|---|---|

$\Delta \xafabs$ | 0.48 | 0.98 | 0.45 | 0.10 |

Δ_{max} | 0.82 | 1.86 | 0.96 | 0.52 |

$\Delta \xaf$ | −0.48 | 0.80 | −0.45 | −0.02 |

Δ_{std} | 0.27 | 0.93 | 0.30 | 0.13 |

Rydberg states . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . |
---|---|---|---|---|

$\Delta \xafabs$ | 0.48 | 0.98 | 0.45 | 0.10 |

Δ_{max} | 0.82 | 1.86 | 0.96 | 0.52 |

$\Delta \xaf$ | −0.48 | 0.80 | −0.45 | −0.02 |

Δ_{std} | 0.27 | 0.93 | 0.30 | 0.13 |

$|RS(2)|<1.0$ . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . |
---|---|---|---|---|

$\Delta \xafabs$ | 0.31 | 0.84 | 0.28 | 0.08 |

Δ_{max} | 0.81 | 1.86 | 0.96 | 0.33 |

$\Delta \xaf$ | −0.30 | 0.68 | −0.22 | −0.04 |

Δ_{std} | 0.25 | 0.80 | 0.33 | 0.09 |

$|RS(2)|<1.0$ . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . |
---|---|---|---|---|

$\Delta \xafabs$ | 0.31 | 0.84 | 0.28 | 0.08 |

Δ_{max} | 0.81 | 1.86 | 0.96 | 0.33 |

$\Delta \xaf$ | −0.30 | 0.68 | −0.22 | −0.04 |

Δ_{std} | 0.25 | 0.80 | 0.33 | 0.09 |

$|RS(2)|<0.3$ . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . |
---|---|---|---|---|

$\Delta \xafabs$ | 0.32 | 0.79 | 0.30 | 0.07 |

Δ_{max} | 0.81 | 1.86 | 0.96 | 0.14 |

$\Delta \xaf$ | −0.30 | 0.60 | −0.24 | −0.05 |

Δ_{std} | 0.27 | 0.79 | 0.34 | 0.07 |

$|RS(2)|<0.3$ . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . |
---|---|---|---|---|

$\Delta \xafabs$ | 0.32 | 0.79 | 0.30 | 0.07 |

Δ_{max} | 0.81 | 1.86 | 0.96 | 0.14 |

$\Delta \xaf$ | −0.30 | 0.60 | −0.24 | −0.05 |

Δ_{std} | 0.27 | 0.79 | 0.34 | 0.07 |

For a better visualization of the results, we report in Fig. 5 the normal distribution of the errors taking CCSD excitation energies as reference. In addition, we have used the $|RS(2)|$ diagnostic to plot only results with $|RS(2)|<1.0$ and $|RS(2)|<0.3$ for the CPS(D-3) model. In Fig. 5, the normal distribution of the errors of the different models confirms the clear superiority of the CPS(D-3) model over the second-order models in terms of accuracy. The reliability of the $|RS(2)|$ diagnostic is also nicely exemplified since the CPS(D-3) error distributions become more narrow as we remove states with large $|RS(2)|$ values.

In order to put our investigation in a broader perspective, we compare our results with CC3 excitation energies that accurately describe the effect of triple excitations in excitation energy calculations.^{30} Table IX contains the statistical measures for the CC2, CCS, CPS(D-2), CPS(D-3) and CCSD models compared to CC3 results, while the normal distribution of the errors are presented in Fig. 6. First we note that for single-replacement dominated excitations with $|RS(2)|<0.3$ we have statistical errors comparing CPS(D-3) excitation energies with CCSD excitation energies (see Table VIII) of the same size as the errors in CCSD excitation energies compared to CC3 excitation energies (see Table IX), in particular the CPS(D-3) maximum error of 0.14 eV is much smaller than maximum error of 0.32 eV for CCSD excitation energies compared to CC3 excitation energies. The CPS(D-3) excitation energies therefore have CCSD quality. From Table IX, we also note the remarkable performance of the CC2 model for valence states ($\Delta \xafabs=0.12$ eV) which has to be seen in perspective with the relatively poor performance for Rydberg states ($\Delta \xafabs=0.44$ eV). On the other hand, the CPS(D-3) model provides a more balanced behaviour with $\Delta \xafabs=0.17$ eV and $\Delta \xafabs=0.08$ eV for valence and Rydberg states, respectively. More importantly, the results from Table IX show that the performance of the CPS(D-3) and CCSD excitation energies compared to the CC3 reference are very similar. For the whole test set we have $\Delta \xafabs=0.12$ eV for CPS(D-3) and $\Delta \xafabs=0.11$ eV for CCSD and when considering only states with $|RS(2)|<0.3$, the CPS(D-3) model outperforms the CCSD results with $\Delta \xaf=0.04$ eV for CPS(D-3) and $\Delta \xaf=0.08$ eV for CCSD. However, the standard deviations and maximum errors are larger for the CPS(D-3) than for the CCSD model. This is confirmed by the normal distributions in Fig. 6 where we also see that all the lower-scaling models benefit from cancellation of errors with the CPS(D-3) model being the most reliable one, i.e. the one with the most narrow distribution.

All states . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . | CCSD . |
---|---|---|---|---|---|

$\Delta \xafabs$ | 0.30 | 1.07 | 0.31 | 0.12 | 0.11 |

Δ_{max} | 0.82 | 3.04 | 0.97 | 0.59 | 0.39 |

$\Delta \xaf$ | −0.26 | 0.96 | −0.16 | 0.08 | 0.10 |

Δ_{std} | 0.31 | 0.89 | 0.38 | 0.15 | 0.10 |

All states . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . | CCSD . |
---|---|---|---|---|---|

$\Delta \xafabs$ | 0.30 | 1.07 | 0.31 | 0.12 | 0.11 |

Δ_{max} | 0.82 | 3.04 | 0.97 | 0.59 | 0.39 |

$\Delta \xaf$ | −0.26 | 0.96 | −0.16 | 0.08 | 0.10 |

Δ_{std} | 0.31 | 0.89 | 0.38 | 0.15 | 0.10 |

Valence states . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . | CCSD . |
---|---|---|---|---|---|

$\Delta \xafabs$ | 0.12 | 1.18 | 0.19 | 0.17 | 0.18 |

Δ_{max} | 0.57 | 3.04 | 0.57 | 0.59 | 0.39 |

$\Delta \xaf$ | −0.03 | 1.10 | 0.16 | 0.17 | 0.18 |

Δ_{std} | 0.18 | 0.83 | 0.18 | 0.14 | 0.09 |

Valence states . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . | CCSD . |
---|---|---|---|---|---|

$\Delta \xafabs$ | 0.12 | 1.18 | 0.19 | 0.17 | 0.18 |

Δ_{max} | 0.57 | 3.04 | 0.57 | 0.59 | 0.39 |

$\Delta \xaf$ | −0.03 | 1.10 | 0.16 | 0.17 | 0.18 |

Δ_{std} | 0.18 | 0.83 | 0.18 | 0.14 | 0.09 |

Rydberg states . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . | CCSD . |
---|---|---|---|---|---|

$\Delta \xafabs$ | 0.44 | 0.99 | 0.41 | 0.08 | 0.05 |

Δ_{max} | 0.82 | 1.98 | 0.97 | 0.52 | 0.12 |

$\Delta \xaf$ | −0.44 | 0.84 | −0.40 | 0.02 | 0.04 |

Δ_{std} | 0.27 | 0.93 | 0.30 | 0.13 | 0.03 |

Rydberg states . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . | CCSD . |
---|---|---|---|---|---|

$\Delta \xafabs$ | 0.44 | 0.99 | 0.41 | 0.08 | 0.05 |

Δ_{max} | 0.82 | 1.98 | 0.97 | 0.52 | 0.12 |

$\Delta \xaf$ | −0.44 | 0.84 | −0.40 | 0.02 | 0.04 |

Δ_{std} | 0.27 | 0.93 | 0.30 | 0.13 | 0.03 |

$|RS(2)|<1.0$ . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . | CCSD . |
---|---|---|---|---|---|

$\Delta \xafabs$ | 0.25 | 0.90 | 0.29 | 0.10 | 0.10 |

Δ_{max} | 0.82 | 2.02 | 0.97 | 0.59 | 0.32 |

$\Delta \xaf$ | −0.20 | 0.78 | −0.12 | 0.06 | 0.10 |

Δ_{std} | 0.28 | 0.83 | 0.38 | 0.14 | 0.09 |

$|RS(2)|<1.0$ . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . | CCSD . |
---|---|---|---|---|---|

$\Delta \xafabs$ | 0.25 | 0.90 | 0.29 | 0.10 | 0.10 |

Δ_{max} | 0.82 | 2.02 | 0.97 | 0.59 | 0.32 |

$\Delta \xaf$ | −0.20 | 0.78 | −0.12 | 0.06 | 0.10 |

Δ_{std} | 0.28 | 0.83 | 0.38 | 0.14 | 0.09 |

$|RS(2)|<0.3$ . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . | CCSD . |
---|---|---|---|---|---|

$\Delta \xafabs$ | 0.27 | 0.82 | 0.31 | 0.08 | 0.09 |

Δ_{max} | 0.82 | 2.02 | 0.97 | 0.43 | 0.32 |

$\Delta \xaf$ | −0.22 | 0.69 | −0.16 | 0.04 | 0.08 |

Δ_{std} | 0.30 | 0.81 | 0.38 | 0.11 | 0.08 |

$|RS(2)|<0.3$ . | CC2 . | CCS . | CPS(D-2) . | CPS(D-3) . | CCSD . |
---|---|---|---|---|---|

$\Delta \xafabs$ | 0.27 | 0.82 | 0.31 | 0.08 | 0.09 |

Δ_{max} | 0.82 | 2.02 | 0.97 | 0.43 | 0.32 |

$\Delta \xaf$ | −0.22 | 0.69 | −0.16 | 0.04 | 0.08 |

Δ_{std} | 0.30 | 0.81 | 0.38 | 0.11 | 0.08 |

The analysis performed here shows and confirms that the CPS(D-3) model provides excitation energies of CCSD quality that outperforms the CC2 and CPS(D-2) models in terms of accuracy. The convergence of the series can be analysed by looking at the intermediate CCS and CPS(D-2) results. More importantly, outliers in the CPS(D-3) results due to non-single-replacement dominated excitations can be identified by computing the norm of the second-order singles response amplitudes ($|RS(2)|$). Comparison with the more accurate CC3 model shows that the CPS(D-3) excitation energies benefit from a small cancellation of errors. The CPS(D-3) model is thus a very attractive alternative to the CCSD model for the calculation of excitation energies.

### C. Performance investigations

In order to evaluate the performance of the CPS(D-3) model in terms of computational cost, we have performed a series of calculations on saturated fatty acids of increasing length, CH_{3}(CH_{2})_{n}COOH, where *n* = 0, 2, 6, 10, 14, 18, 22. For all molecules, the lowest excitation energy has been calculated at the CC2, CPS(D-2), CPS(D-3), and CCSD level using the aug-cc-pVDZ′ (Fig. 7) and aug-cc-pVTZ′ (Fig. 8) basis sets. The prime in the basis set notation indicates that diffuse functions were not included for hydrogen atoms. The CCSD calculations have been performed with the Turbomole program package,^{50} while for all other calculations the LSDALTON program was used.^{39,40} All calculations were performed on a single Lenovo nx360 M5 node with 28 cores @ 2.4 GHz and 256 GB memory. We note that due to the high computational requirements of the CCSD model, the largest systems (with *n* > 10) could not be treated at that level of theory.

From Figs. 7 and 8, we see that for large systems, as expected, the CPS(D-3) model is much more expensive than the second-order models. For example, for the largest fatty acid, CH_{3}(CH_{2})_{22}COOH in the aug-cc-pVTZ′ basis (1868 basis functions), the CPS(D-3) calculation takes 213 h, while only 24 min and 7 h are required for the CPS(D-2) and CC2 models, respectively. This difference is of course due to the difference in the computational scaling of the methods [$O(N6)$ versus $O(N5)$].

However, the CCSD timings put things in perspective and we see that even though the CPS(D-3) and CCSD models formally have the same scaling, the difference in the pre-factor is of great importance in practice. The CPS(D-3) model can thus be applied to much larger systems than the CCSD model and with excitation energies of similar quality (see Secs. III A and III B). This is confirmed in Fig. 8 where a CPS(D-3)/aug-cc-pVTZ′ calculation for the Lignoceric acid (1868 basis functions) is reported.

## IV. SUMMARY AND CONCLUSION

Excitation energies have been determined using the cluster perturbation series CPS(D), where a series of excitation energy corrections is determined in orders of the fluctuation potential that added to a CCS excitation energy formally converges to a CCSD excitation energy. In particular, we have investigated the performance of the low-order models of the CPS(D) series from a cost benefit point of view, where the computational cost and the convergence of the individual models have been taken into account.

The CPS(D) series of excitation energies is only well behaved if the targeted CCSD excited state is well described at the CCS level. We denote such excitations as single-replacement dominated excitations and have described how single-replacement excitations can be identified. We introduced the $|RS(2)|$ diagnostic for their identification.

For single-replacement dominated excitations, we have confirmed that the second-order model CPS(D-2) gives a good estimate of the doubles correction to the CCS excitation energies. Further, by carrying out calculations for a large variety of excitation energies, we have found that the third-order model CPS(D-3) gives excitation energies of CCSD quality in the sense that the difference between CPS(D-3) and CCSD excitation energies is of the same size or smaller than the effect of adding triples corrections to CCSD excitation energies. By comparing the accuracy of CPS(D-3) and CCSD excitation energies relative to excitation energy calculations where triple excitations are considered, we have also found that the performance of the CPS(D-3) and CCSD models is very similar. The CPS(D-3) model thus constitutes a very attractive alternative to the CCSD model for determining excitation energies. For single-replacement dominated excitations, we thus recommend that CCSD calculations are replaced by CPS(D-3) calculations if triple excitation energy corrections are not calculated explicitly.

CPS(D-3) excitation energies can also be determined for system sizes that are far beyond what can be considered in conventional CCSD excitation energy calculations (i.e., when the CCSD amplitude and CCSD response eigenvalue equations are solved explicitly). This has several reasons. Both models have the same leading-order computational scaling, however, the CPS(D-3) calculations have a pre-factor that is ∼3/45 of the one of a conventional CCSD calculation. Equally important, the equations for the CPS(D-3) model, in contrast to the equations for the CCSD model, contain only a product of a first-order amplitude and a two-electron integral. This reduces dramatically the storage requirement in the CPS(D-3) excitation calculations compared to a conventional CCSD calculation, in particular, when the resolution of the identity (RI) approximation is used for the integrals.

We have described an efficient implementation of the CPS(D-3) model for calculating excitation energies, where we have used the RI approximation for integrals and we have further discussed how this implementation may straightforwardly be parallelized. This will dramatically extend the system size for which CPS(D-3) excitation energies can be determined.

In this paper, we have demonstrated the potential of using CP theory to calculate excitation energies in terms of accuracy and computational efficiency. Future work may involve using the CP series CPSD(T) to efficiently calculate triples corrections to CCSD excitation energies and to use CP theory to calculate other molecular properties. An obvious extension of this work would be to use the CPS(D) model for the calculation of transition strength and excited state molecular properties.

## SUPPLEMENTARY MATERIAL

See supplementary material for the convergence of the CPS(D) excitation energy series for the three lowest singlet and triplet states of H_{2}O, the three lowest singlet states of CO, and three lowest triplet–triplet transitions of O_{2}. In addition, the supplementary material contains individual excitation energies at the CC2, CCS, CPS(D-2), CPS(D-3), CCSD and CC3 level used in Sec. III B. The assignment and the character of each transition is also reported.

## ACKNOWLEDGMENTS

P.B., F.P., D.B., and K.K. acknowledge qLEAP Center for Theoretical Chemistry, Department of Chemistry, Aarhus University, Langelandsgade 140, DK-8000 Aarhus C, Denmark, where they were formerly employed and where research leading to this article was carried out. P.B., F.P., D.B., K.K., and P.J. acknowledge support from The European Research Council under the European Union’s (EU) Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement No. 291371. D.B. acknowledges the Marie Curie Individual Fellowship funding for “DECOS,” Project No. 657514. J.O. acknowledges support from the Danish Council for Independent Research, Grant No. DFF-4181-00537. This research used resources of the Centre for Scientific Computing, Aarhus, Denmark (http://phys.au.dk/forskning/cscaa/) and resources provided by Wroclaw Centre for Networking and Supercomputing, Poland (http://www.wcss.pl/en/), Grant No. 277. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.

### APPENDIX: CPS(D) EQUATIONS FOR EXCITATION ENERGIES THROUGH SIXTH-ORDER

In this appendix we gather the equations necessary to calculate excitation energies with the CPS(D) series up to sixth-order.

The ground-state singles and doubles cluster amplitudes through fifth-order are given in Table X.

The Jacobian matrix elements through sixth-order are given in Table XI.

The CPS(D) excitation energy expressions through sixth-order are given in Table XII.

The right singles excitation vectors through fifth-order are given in Table XIII.

The right doubles excitation vectors through fifth-order are given in Table XIV.

$t\mu 1(1)$ | = | 0 |

$\u2211\nu 1J\mu 1\nu 1CCSt\nu 1(2)$ | = | $\u2212\u2009\mu 1|\Phi ,T2(1)|HF\u2009$ |

$\u2211\nu 1J\mu 1\nu 1CCSt\nu 1(3)$ | = | $\u2212\u2009\mu 1|\Phi ,T2(2)|HF\u2009$ |

$\u2211\nu 1J\mu 1\nu 1CCSt\nu 1(4)$ | = | $\u2212\u2009\mu 1|\Phi ,T2(3)|HF\u2009\u2212\u2009\mu 1|\Phi ,T1(2),T2(1)|HF\u2009$ |

$\u2211\nu 1J\mu 1\nu 1CCSt\nu 1(5)$ | = | $\u2212\u2009\mu 1|\Phi ,T2(4)|HF\u2009\u2212\u2009\mu 1|\Phi ,T1(3),T2(1)|HF\u2009\u2212\u2009\mu 1|\Phi ,T1(2),T2(2)|HF\u2009$ |

$\u221212\u2009\mu 1|\Phi ,T1(2),T1(2)|HF\u2009$ | ||

$\epsilon \mu 2t\mu 2(1)$ | = | −⟨μ_{2}|Φ|HF⟩ |

$\epsilon \mu 2t\mu 2(2)$ | = | $\u2212\u2009\mu 2|\Phi ,T2(1)|HF\u2009$ |

$\epsilon \mu 2t\mu 2(3)$ | = | $\u2212\u2009\mu 2|\Phi ,T1(2)+T2(2)|HF\u2009\u221212\u2009\mu 2|\Phi ,T2(1),T2(1)|HF\u2009$ |

$\epsilon \mu 2t\mu 2(4)$ | = | $\u2212\u2009\mu 2|\Phi ,T1(3)+T2(3)|HF\u2009\u2212\u2009\mu 2|\Phi ,T2(1),T1(2)+T2(2)|HF\u2009$ |

$\epsilon \mu 2t\mu 2(5)$ | = | $\u2212\u2009\mu 2|\Phi ,T1(4)+T2(4)|HF\u2009\u2212\u2009\mu 2|\Phi ,T2(1),T1(3)+T2(3)|HF\u2009$ |

$\u221212\u2009\mu 2|\Phi ,T1(2)+T2(2),T1(2)+T2(2)|HF\u2009$ |

$t\mu 1(1)$ | = | 0 |

$\u2211\nu 1J\mu 1\nu 1CCSt\nu 1(2)$ | = | $\u2212\u2009\mu 1|\Phi ,T2(1)|HF\u2009$ |

$\u2211\nu 1J\mu 1\nu 1CCSt\nu 1(3)$ | = | $\u2212\u2009\mu 1|\Phi ,T2(2)|HF\u2009$ |

$\u2211\nu 1J\mu 1\nu 1CCSt\nu 1(4)$ | = | $\u2212\u2009\mu 1|\Phi ,T2(3)|HF\u2009\u2212\u2009\mu 1|\Phi ,T1(2),T2(1)|HF\u2009$ |

$\u2211\nu 1J\mu 1\nu 1CCSt\nu 1(5)$ | = | $\u2212\u2009\mu 1|\Phi ,T2(4)|HF\u2009\u2212\u2009\mu 1|\Phi ,T1(3),T2(1)|HF\u2009\u2212\u2009\mu 1|\Phi ,T1(2),T2(2)|HF\u2009$ |

$\u221212\u2009\mu 1|\Phi ,T1(2),T1(2)|HF\u2009$ | ||

$\epsilon \mu 2t\mu 2(1)$ | = | −⟨μ_{2}|Φ|HF⟩ |

$\epsilon \mu 2t\mu 2(2)$ | = | $\u2212\u2009\mu 2|\Phi ,T2(1)|HF\u2009$ |

$\epsilon \mu 2t\mu 2(3)$ | = | $\u2212\u2009\mu 2|\Phi ,T1(2)+T2(2)|HF\u2009\u221212\u2009\mu 2|\Phi ,T2(1),T2(1)|HF\u2009$ |

$\epsilon \mu 2t\mu 2(4)$ | = | $\u2212\u2009\mu 2|\Phi ,T1(3)+T2(3)|HF\u2009\u2212\u2009\mu 2|\Phi ,T2(1),T1(2)+T2(2)|HF\u2009$ |

$\epsilon \mu 2t\mu 2(5)$ | = | $\u2212\u2009\mu 2|\Phi ,T1(4)+T2(4)|HF\u2009\u2212\u2009\mu 2|\Phi ,T2(1),T1(3)+T2(3)|HF\u2009$ |

$\u221212\u2009\mu 2|\Phi ,T1(2)+T2(2),T1(2)+T2(2)|HF\u2009$ |

$J\mu i\nu j(0)$ | = | $J\mu i\nu jCCS\delta i1\delta j1+\epsilon \nu j\delta \mu i\nu j\delta i2\delta j2$ |

$J\mu i\nu j(1)$ | = | $\u2009\mu i|\Phi ,\theta \nu j|HF\u2009\delta i1\delta j2+\u2009\mu i|\Phi ,\theta \nu j|HF\u2009\delta i2\delta j1+\u2009\mu i|\Phi ,\theta \nu j|HF\u2009\delta i2\delta j2$ |

$J\mu i\nu j(2)$ | = | $\u2009\mu i|\Phi ,T2(1),\theta \nu j|HF\u2009$ |

$J\mu i\nu j(3)$ | = | $\u2009\mu i|\Phi ,T(2),\theta \nu j|HF\u2009$ |

$J\mu i\nu j(4)$ | = | $\u2009\mu i|\Phi ,T(3),\theta \nu j|HF\u2009+\u2009\mu i|\Phi ,T2(1),T(2),\theta \nu j|HF\u2009\delta i2$ |

$J\mu i\nu j(5)$ | = | $\u2009\mu i|\Phi ,T(4),\theta \nu j|HF\u2009\u2009+\u2009\u2009\mu i|\Phi ,T2(1),T(3),\theta \nu j|HF\u2009\delta i2\u2009+\u200912\u2009\mu i|\Phi ,T(2),T(2),\theta \nu j|HF\u2009$ |

$J\mu i\nu j(6)$ | = | $\u2009\mu i|\Phi ,T(5),\theta \nu j|HF\u2009\u2009+\u2009\u2009\mu i|\Phi ,T2(1),T(4),\theta \nu j|HF\u2009\delta i2\u2009+\u2009\u2009\mu i|\Phi ,T(2),T(3),\theta \nu j|HF\u2009$ |

$J\mu i\nu j(0)$ | = | $J\mu i\nu jCCS\delta i1\delta j1+\epsilon \nu j\delta \mu i\nu j\delta i2\delta j2$ |

$J\mu i\nu j(1)$ | = | $\u2009\mu i|\Phi ,\theta \nu j|HF\u2009\delta i1\delta j2+\u2009\mu i|\Phi ,\theta \nu j|HF\u2009\delta i2\delta j1+\u2009\mu i|\Phi ,\theta \nu j|HF\u2009\delta i2\delta j2$ |

$J\mu i\nu j(2)$ | = | $\u2009\mu i|\Phi ,T2(1),\theta \nu j|HF\u2009$ |

$J\mu i\nu j(3)$ | = | $\u2009\mu i|\Phi ,T(2),\theta \nu j|HF\u2009$ |

$J\mu i\nu j(4)$ | = | $\u2009\mu i|\Phi ,T(3),\theta \nu j|HF\u2009+\u2009\mu i|\Phi ,T2(1),T(2),\theta \nu j|HF\u2009\delta i2$ |

$J\mu i\nu j(5)$ | = | $\u2009\mu i|\Phi ,T(4),\theta \nu j|HF\u2009\u2009+\u2009\u2009\mu i|\Phi ,T2(1),T(3),\theta \nu j|HF\u2009\delta i2\u2009+\u200912\u2009\mu i|\Phi ,T(2),T(2),\theta \nu j|HF\u2009$ |

$J\mu i\nu j(6)$ | = | $\u2009\mu i|\Phi ,T(5),\theta \nu j|HF\u2009\u2009+\u2009\u2009\mu i|\Phi ,T2(1),T(4),\theta \nu j|HF\u2009\delta i2\u2009+\u2009\u2009\mu i|\Phi ,T(2),T(3),\theta \nu j|HF\u2009$ |

$\omega x(1)$ | = | 0 |

$\omega x(2)$ | = | $\u2009LxCCS|\Phi ,T2(1),RxCCS|HF\u2009+\u2009LxCCS|\Phi ,Rx2(1)|HF\u2009$ |

$\omega x(3)$ | = | $\u2009LxCCS|\Phi ,T(2),RxCCS|HF\u2009+\u2009LxCCS|\Phi ,Rx2(2)|HF\u2009$ |

$\omega x(4)$ | = | $\u2009LxCCS|\Phi ,T(3),RxCCS|HF\u2009+\u2009LxCCS|\Phi ,Rx2(3)|HF\u2009+\u2009LxCCS|\Phi ,T2(1),Rx1(2)|HF\u2009$ |

$+\u2009\u2009LxCCS|\Phi ,T(2),Rx2(1)|HF\u2009$ | ||

$\omega x(5)$ | = | $\u2009LxCCS|\Phi ,T(4),RxCCS|HF\u2009+12\u2009LxCCS|\Phi ,T1(2),T1(2),RxCCS|HF\u2009+\u2009LxCCS|\Phi ,Rx2(4)|HF\u2009$ |

$+\u2009\u2009LxCCS|\Phi ,T2(1),Rx1(3)|HF\u2009+\u2009LxCCS|\Phi ,T(2),Rx1(2)|HF\u2009+\u2009LxCCS|\Phi ,T1(2),Rx2(2)|HF\u2009$ | ||

$+\u2009\u2009LxCCS|\Phi ,T1(3),Rx2(1)|HF\u2009$ | ||

$\omega x(6)$ | = | $\u2009LxCCS|\Phi ,T(5),RxCCS|HF\u2009+\u2009LxCCS|\Phi ,T1(2),T1(3),RxCCS|HF\u2009+\u2009LxCCS|\Phi ,Rx2(5)|HF\u2009$ |

$+\u2009\u2009LxCCS|\Phi ,T2(1),Rx1(4)|HF\u2009+\u2009LxCCS|\Phi ,T(2),Rx1(3)|HF\u2009+\u2009LxCCS|\Phi ,T1(2),Rx2(3)|HF\u2009$ | ||

$+\u2009\u2009LxCCS|\Phi ,T(3),Rx1(2)|HF\u2009+\u2009LxCCS|\Phi ,T1(3),Rx2(2)|HF\u2009+\u2009LxCCS|\Phi ,T1(4),Rx2(1)|HF\u2009$ |

$\omega x(1)$ | = | 0 |

$\omega x(2)$ | = | $\u2009LxCCS|\Phi ,T2(1),RxCCS|HF\u2009+\u2009LxCCS|\Phi ,Rx2(1)|HF\u2009$ |

$\omega x(3)$ | = | $\u2009LxCCS|\Phi ,T(2),RxCCS|HF\u2009+\u2009LxCCS|\Phi ,Rx2(2)|HF\u2009$ |

$\omega x(4)$ | = | $\u2009LxCCS|\Phi ,T(3),RxCCS|HF\u2009+\u2009LxCCS|\Phi ,Rx2(3)|HF\u2009+\u2009LxCCS|\Phi ,T2(1),Rx1(2)|HF\u2009$ |

$+\u2009\u2009LxCCS|\Phi ,T(2),Rx2(1)|HF\u2009$ | ||

$\omega x(5)$ | = | $\u2009LxCCS|\Phi ,T(4),RxCCS|HF\u2009+12\u2009LxCCS|\Phi ,T1(2),T1(2),RxCCS|HF\u2009+\u2009LxCCS|\Phi ,Rx2(4)|HF\u2009$ |

$+\u2009\u2009LxCCS|\Phi ,T2(1),Rx1(3)|HF\u2009+\u2009LxCCS|\Phi ,T(2),Rx1(2)|HF\u2009+\u2009LxCCS|\Phi ,T1(2),Rx2(2)|HF\u2009$ | ||

$+\u2009\u2009LxCCS|\Phi ,T1(3),Rx2(1)|HF\u2009$ | ||

$\omega x(6)$ | = | $\u2009LxCCS|\Phi ,T(5),RxCCS|HF\u2009+\u2009LxCCS|\Phi ,T1(2),T1(3),RxCCS|HF\u2009+\u2009LxCCS|\Phi ,Rx2(5)|HF\u2009$ |

$+\u2009\u2009LxCCS|\Phi ,T2(1),Rx1(4)|HF\u2009+\u2009LxCCS|\Phi ,T(2),Rx1(3)|HF\u2009+\u2009LxCCS|\Phi ,T1(2),Rx2(3)|HF\u2009$ | ||

$+\u2009\u2009LxCCS|\Phi ,T(3),Rx1(2)|HF\u2009+\u2009LxCCS|\Phi ,T1(3),Rx2(2)|HF\u2009+\u2009LxCCS|\Phi ,T1(4),Rx2(1)|HF\u2009$ |

$Rx\nu 1(1)$ | = | 0 |

$\u2211\nu 1J\mu 1\nu 1CCS\u2212\omega xCCS\u2009\delta \mu 1\nu 1Rx\nu 1(2)$ | = | $\omega x(2)Rx\mu 1CCS\u2212\u2009\mu 1|\Phi ,T2(1),RxCCS|HF\u2009\u2212\u2009\mu 1|\Phi ,Rx2(1)|HF\u2009$ |

$\u2211\nu 1J\mu 1\nu 1CCS\u2212\omega xCCS\u2009\delta \mu 1\nu 1Rx\nu 1(3)$ | = | $\omega x(3)Rx\mu 1CCS\u2212\u2009\mu 1|\Phi ,T(2),RxCCS|HF\u2009\u2212\u2009\mu 1|\Phi ,Rx2(2)|HF\u2009$ |

$\u2211\nu 1J\mu 1\nu 1CCS\u2212\omega xCCS\u2009\delta \mu 1\nu 1Rx\nu 1(4)$ | = | $\omega x(4)Rx\mu 1CCS+\omega x(2)Rx\mu 1(2)\u2212\u2009\mu 1|\Phi ,T(3),RxCCS|HF\u2009$ |

$\u2212\u2009\u2009\mu 1|\Phi ,T2(1),Rx1(2)|HF\u2009\u2212\u2009\mu 1|\Phi ,Rx2(3)|HF\u2009$ | ||

$\u2212\u2009\u2009\mu 1|\Phi ,T(2),Rx2(1)|HF\u2009$ | ||

$\u2211\nu 1J\mu 1\nu 1CCS\u2212\omega xCCS\u2009\delta \mu 1\nu 1Rx\nu 1(5)$ | = | $\omega x(5)Rx\mu 1CCS+\omega x(2)Rx\mu 1(3)+\omega x(3)Rx\mu 1(2)\u2212\u2009\mu 1|\Phi ,T(4),RxCCS|HF\u2009$ |

$\u2212\u200912\u2009\mu 1|\Phi ,T1(2),T1(2),RxCCS|HF\u2009\u2212\u2009\mu 1|\Phi ,T2(1),Rx1(3)|HF\u2009$ | ||

$\u2212\u2009\u2009\mu 1|\Phi ,T(2),Rx1(2)|HF\u2009\u2212\u2009\mu 1|\Phi ,Rx2(4)|HF\u2009\u2212\u2009\mu 1|\Phi ,Rx2(4)|HF\u2009$ | ||

$\u2212\u2009\u2009\mu 1|\Phi ,T1(2),Rx2(2)|HF\u2009\u2212\u2009\mu 1|\Phi ,T1(3),Rx2(1)|HF\u2009$ |

$Rx\nu 1(1)$ | = | 0 |

$\u2211\nu 1J\mu 1\nu 1CCS\u2212\omega xCCS\u2009\delta \mu 1\nu 1Rx\nu 1(2)$ | = | |

$\u2211\nu 1J\mu 1\nu 1CCS\u2212\omega xCCS\u2009\delta \mu 1\nu 1Rx\nu 1(3)$ | = | $\omega x(3)Rx\mu 1CCS\u2212\u2009\mu 1|\Phi ,T(2),RxCCS|HF\u2009\u2212\u2009\mu 1|\Phi ,Rx2(2)|HF\u2009$ |

$\u2211\nu 1J\mu 1\nu 1CCS\u2212\omega xCCS\u2009\delta \mu 1\nu 1Rx\nu 1(4)$ | = | $\omega x(4)Rx\mu 1CCS+\omega x(2)Rx\mu 1(2)\u2212\u2009\mu 1|\Phi ,T(3),RxCCS|HF\u2009$ |

$\u2212\u2009\u2009\mu 1|\Phi ,T2(1),Rx1(2)|HF\u2009\u2212\u2009\mu 1|\Phi ,Rx2(3)|HF\u2009$ | ||

$\u2212\u2009\u2009\mu 1|\Phi ,T(2),Rx2(1)|HF\u2009$ | ||

$\u2211\nu 1J\mu 1\nu 1CCS\u2212\omega xCCS\u2009\delta \mu 1\nu 1Rx\nu 1(5)$ | = | $\omega x(5)Rx\mu 1CCS+\omega x(2)Rx\mu 1(3)+\omega x(3)Rx\mu 1(2)\u2212\u2009\mu 1|\Phi ,T(4),RxCCS|HF\u2009$ |

$\u2212\u200912\u2009\mu 1|\Phi ,T1(2),T1(2),RxCCS|HF\u2009\u2212\u2009\mu 1|\Phi ,T2(1),Rx1(3)|HF\u2009$ | ||

$\u2212\u2009\u2009\mu 1|\Phi ,T(2),Rx1(2)|HF\u2009\u2212\u2009\mu 1|\Phi ,Rx2(4)|HF\u2009\u2212\u2009\mu 1|\Phi ,Rx2(4)|HF\u2009$ | ||

$\u2212\u2009\u2009\mu 1|\Phi ,T1(2),Rx2(2)|HF\u2009\u2212\u2009\mu 1|\Phi ,T1(3),Rx2(1)|HF\u2009$ |

$\epsilon \mu 2\u2212\omega xCCSRx\mu 2(1)$ | = | $\u2212\u2009\mu 2|\Phi ,RxCCS|HF\u2009$ |

$\epsilon \mu 2\u2212\omega xCCSRx\mu 2(2)$ | = | $\u2212\u2009\mu 2|\Phi ,T2(1),RxCCS|HF\u2009\u2212\u2009\mu 2|\Phi ,Rx2(1)|HF\u2009$ |

$\epsilon \mu 2\u2212\omega xCCSRx\mu 2(3)$ | = | $\omega x(2)Rx\mu 2(1)\u2212\u2009\mu 2|\Phi ,T(2),RxCCS|HF\u2009\u2212\u2009\mu 2|\Phi ,Rx1(2)|HF\u2009\u2212\u2009\mu 2|\Phi ,Rx2(2)|HF\u2009$ |

$\u2212\u2009\u2009\mu 2|\Phi ,T2(1),Rx2(1)|HF\u2009$ | ||

$\epsilon \mu 2\u2212\omega xCCSRx\mu 2(4)$ | = | $\omega x(2)Rx\mu 2(2)+\omega x(3)Rx\mu 2(1)\u2212\u2009\mu 2|\Phi ,T(3),RxCCS|HF\u2009$ |

$\u2212\u2009\u2009\mu 2|\Phi ,T2(1),T(2),RxCCS|HF\u2009\u2212\u2009\mu 2|\Phi ,Rx1(3)|HF\u2009$ | ||

$\u2212\u2009\u2009\mu 2|\Phi ,T2(1),Rx1(2)|HF\u2009\u2212\u2009\mu 2|\Phi ,Rx2(3)|HF\u2009$ | ||

$\u2212\u2009\u2009\mu 2|\Phi ,T2(1),Rx2(2)|HF\u2009\u2212\u2009\mu 2|\Phi ,T(2),Rx2(1)|HF\u2009$ | ||

$\epsilon \mu 2\u2212\omega xCCSRx\mu 2(5)$ | = | $\omega x(2)Rx\mu 2(3)+\omega x(3)Rx\mu 2(2)+\omega x(4)Rx\mu 2(1)\u2212\u2009\mu 2|\Phi ,T(4),RxCCS|HF\u2009$ |

$\u2212\u2009\u2009\mu 2|\Phi ,T2(1),T1(3),RxCCS|HF\u2009\u221212\u2009\mu 2|\Phi ,T(2),T(2),RxCCS|HF\u2009$ | ||

$\u2212\u2009\u2009\mu 2|\Phi ,Rx1(4)|HF\u2009\u2212\u2009\mu 2|\Phi ,T2(1),Rx1(3)|HF\u2009\u2212\u2009\mu 2|\Phi ,T(2),Rx1(2)|HF\u2009$ | ||

$\u2212\u2009\u2009\mu 2|\Phi ,Rx2(4)|HF\u2009\u2212\u2009\mu 2|\Phi ,T2(1),Rx2(3)|HF\u2009\u2212\u2009\mu 2|\Phi ,T(2),Rx2(2)|HF\u2009$ | ||

$\u2212\u2009\u2009\mu 2|\Phi ,T(3),Rx2(1)|HF\u2009$ |

$\epsilon \mu 2\u2212\omega xCCSRx\mu 2(1)$ | = | $\u2212\u2009\mu 2|\Phi ,RxCCS|HF\u2009$ |

$\epsilon \mu 2\u2212\omega xCCSRx\mu 2(2)$ | = | $\u2212\u2009\mu 2|\Phi ,T2(1),RxCCS|HF\u2009\u2212\u2009\mu 2|\Phi ,Rx2(1)|HF\u2009$ |

$\epsilon \mu 2\u2212\omega xCCSRx\mu 2(3)$ | = | $\omega x(2)Rx\mu 2(1)\u2212\u2009\mu 2|\Phi ,T(2),RxCCS|HF\u2009\u2212\u2009\mu 2|\Phi ,Rx1(2)|HF\u2009\u2212\u2009\mu 2|\Phi ,Rx2(2)|HF\u2009$ |

$\u2212\u2009\u2009\mu 2|\Phi ,T2(1),Rx2(1)|HF\u2009$ | ||

$\epsilon \mu 2\u2212\omega xCCSRx\mu 2(4)$ | = | $\omega x(2)Rx\mu 2(2)+\omega x(3)Rx\mu 2(1)\u2212\u2009\mu 2|\Phi ,T(3),RxCCS|HF\u2009$ |

$\u2212\u2009\u2009\mu 2|\Phi ,T2(1),T(2),RxCCS|HF\u2009\u2212\u2009\mu 2|\Phi ,Rx1(3)|HF\u2009$ | ||

$\u2212\u2009\u2009\mu 2|\Phi ,T2(1),Rx1(2)|HF\u2009\u2212\u2009\mu 2|\Phi ,Rx2(3)|HF\u2009$ | ||

$\u2212\u2009\u2009\mu 2|\Phi ,T2(1),Rx2(2)|HF\u2009\u2212\u2009\mu 2|\Phi ,T(2),Rx2(1)|HF\u2009$ | ||

$\epsilon \mu 2\u2212\omega xCCSRx\mu 2(5)$ | = | $\omega x(2)Rx\mu 2(3)+\omega x(3)Rx\mu 2(2)+\omega x(4)Rx\mu 2(1)\u2212\u2009\mu 2|\Phi ,T(4),RxCCS|HF\u2009$ |

$\u2212\u2009\u2009\mu 2|\Phi ,T2(1),T1(3),RxCCS|HF\u2009\u221212\u2009\mu 2|\Phi ,T(2),T(2),RxCCS|HF\u2009$ | ||

$\u2212\u2009\u2009\mu 2|\Phi ,Rx1(4)|HF\u2009\u2212\u2009\mu 2|\Phi ,T2(1),Rx1(3)|HF\u2009\u2212\u2009\mu 2|\Phi ,T(2),Rx1(2)|HF\u2009$ | ||

$\u2212\u2009\u2009\mu 2|\Phi ,Rx2(4)|HF\u2009\u2212\u2009\mu 2|\Phi ,T2(1),Rx2(3)|HF\u2009\u2212\u2009\mu 2|\Phi ,T(2),Rx2(2)|HF\u2009$ | ||

$\u2212\u2009\u2009\mu 2|\Phi ,T(3),Rx2(1)|HF\u2009$ |