Higher-order diagrammatic vibrational coupled-cluster theory

Faucheaux, Jacob A.; Hirata, So

doi:10.1063/1.4931472

Diagrammatically size-consistent and basis-set-free vibrational coupled-cluster (XVCC) theory for both zero-point energies and transition frequencies of a molecule, the latter through the equation-of-motion (EOM) formalism, is defined for an nth-order Taylor-series potential energy surface (PES). Quantum-field-theoretical tools (the rules of normal-ordered second quantization and Feynman–Goldstone diagrams) for deriving their working equations are established. The equations of XVCC and EOM-XVCC including up to the mth-order excitation operators are derived and implemented with the aid of computer algebra in the range of 1 ≤ m ≤ 8. Algorithm optimizations known as strength reduction, intermediate reuse, and factorization are carried out before code generation, reducing the cost scaling of the mth-order XVCC and EOM-XVCC in an nth-order Taylor-series PES (m ≥ n) to the optimal value of O(N^m+⌊n/2⌋), where N is the number of modes. The calculated zero-point energies and frequencies of fundamentals, overtones, and combinations as well as Fermi-resonant modes display rapid and nearly monotonic convergence with m towards the exact values for the PES. The theory with the same excitation rank as the truncation order of the Taylor-series PES (m = n) seems to strike the best cost-accuracy balance, achieving the accuracy of a few tenths of cm⁻¹ for transitions involving (m − 3) modes and of a few cm⁻¹ for those involving (m − 2) modes. The relationships between XVCC and the vibrational coupled-cluster theories of Prasad and coworkers and of Christiansen and coworkers as well as the size-extensive vibrational self-consistent-field and many-body perturbation theories are also elucidated.

I. INTRODUCTION

Coupled-cluster (CC) theory^1,2 was introduced by Coester and Kümmel^3,4 in nuclear physics and soon adopted by Čížek and coworkers^5,6 in the field of quantum chemistry, where it was fully developed into the most accurate, versatile, and widely applicable electronic structure theory.^7–9 Being rigorously (i.e., diagrammatically) size consistent,^10,11 it was applied to problems in solid-state physics.^12–18 It also witnessed a revival in nuclear physics,^19,20 as the Hamiltonian for nucleons became increasingly more accurate. Another area where the applicability of this theory is being explored is quantum nuclear dynamics in a bound potential energy surface (PES), i.e., anharmonic vibrational structure theory.^21,22

Following its applications to model PES’s,^23–25 vibrational CC (VCC) theory was pioneered and extensively developed by Prasad and coworkers for general molecules with truncated Taylor-series PES’s.^26–34 They considered both zero-point energies and transition frequencies using the equation-of-motion (EOM)^27,32 or multireference^29,34 formalism. They also calculated spectra and band intensities in either the time-dependent^26,28,30,31 or independent³³ framework. These methods are based on a single Hartree product of harmonic-oscillator wave functions as a reference, of which the mode frequencies and coordinate center are varied to minimize its energy. Their working equations are derived using the second-quantized ladder operators, which are valid for harmonic oscillators. The reference wave function adopted can be identified as that of the size-extensive vibrational self-consistent-field (XVSCF) method with an anharmonic geometry correction, i.e., XVSCF[n] (n is the rank of the Taylor-series PES).³⁵

Christiansen and coworkers^36–42 proposed alternative formulations of VCC theory and its excited-state analog, VCC linear response (VCCLR). They may be viewed as a direct extension of vibrational correlation theories pioneered by Bowman and coworkers and by Ratner and coworkers, such as vibrational configuration-interaction (VCI)⁴³ and vibrational Møller–Plesset perturbation (VMP) theories,⁴⁴ and share similar design philosophy. They all rely on a single-reference Hartree product of one-mode functions (“modals”) along normal coordinates furnished by vibrational self-consistent-field (VSCF) theory.^45–47 These modals are determined so as to minimize the VSCF energy and are numerically defined on a grid or a basis set. These vibrational correlation theories tend to be implemented in symbolic (as opposed to matrix-algebraic) algorithms, in which excited Hartree products are manipulated during a calculation (not during formula derivation) with the Hamiltonian matrix elements in modals evaluated by quadrature.

In the meantime, it has been shown that VSCF theory is not diagrammatically size consistent.⁴⁸ Unlike a truncated CI for electrons, which gives zero electron-correlation energy per electron in the thermodynamic limit, VSCF yields a finite and thus physically meaningful anharmonic correction in this limit. Nevertheless, this correction is due to a minute portion of a large number of calculated contributions, most of which are shown to vanish in the same limit as they correspond to disconnected diagrams. This renders VSCF an exceedingly inefficient method for large molecules and solids. When such vanishing contributions are erased from the VSCF equation, we arrive at a streamlined, rigorously size-consistent, but less accurate variant of VSCF, which we call XVSCF(n) or XVSCF[n],^35,49,50 depending on whether an anharmonic geometry correction is suppressed (in the former) or not (in the latter). XVSCF is no longer exact for any one-dimensional potential, which VSCF is, in the infinite basis-set limit. Instead, each modal of XVSCF (or that of VSCF in the thermodynamic limit) is shown to be a harmonic-oscillator wave function with its frequency and center dressed with mean-field anharmonic effects,^48,49 which is consistent with Makri’s theorem.⁵¹ Since the Hamiltonian matrix elements in harmonic-oscillator modals are expressed simply with force constants, XVSCF and vibrational correlation theories based on it can be defined compactly with force constants and excitation amplitudes. Grids, basis sets, and quadrature are no longer necessary.

Insofar as the defining feature of CC theory is its size consistency, which also ensures its uniformly high accuracy across system sizes, its formalism should be stipulated with many-body-theoretical techniques such as normal-ordered second quantization and Feynman–Goldstone diagrams and be based on a rigorously size-consistent reference wave function. A normal-ordered second-quantized form of the vibrational Hamiltonian has been introduced recently.⁵² Remarkably, its constant term and excitation amplitudes are identified as the energy, gradient, and frequency expressions of XVSCF and of the closely related self-consistent-phonon (SCP) method^53–57 (but not of VSCF). This is reminiscent of the normal-ordered electronic Hamiltonian, in which the expressions of the Hartree–Fock (HF) energy and Fock operator naturally emerge.² This observation highlights the central place XVSCF and SCP occupy in anharmonic vibrational structure theory for large molecules and solids and suggests that they provide an ideal reference wave function for vibrational correlation theories for such systems. The establishment of XVSCF,^35,49,50 along with the elucidation of their relationship to the normal-ordered vibrational Hamiltonian,⁵² therefore, sets the stage for a more complete development of a VCC theory.

In this article, we present a rigorously (i.e., diagrammatically) size-consistent VCC theory for the ground and excited states, the latter through the EOM formalism. We call them XVCCm and EOM-XVCCm theories, respectively, where m stands for the highest connected excitation rank considered. We propose a rational and expedient method of deriving their algebraic definitions using the rules of normal-ordered second quantization and Feynman–Goldstone diagrams. We demonstrate its utility through the derivation of low-order XVCC and EOM-XVCC equations. The formalisms thus derived are valid for any single-reference Hartree product of harmonic-oscillator modals along delocalized, rectilinear coordinates, such as XVSCF(n),⁴⁹ XVSCF[n],³⁵ and SCP (Refs. 53–57) as well as the harmonic approximation. Some diagrammatic contributions are shown to vanish in the XVSCF and SCP reference, making them more favorable, although higher-order XVCC and EOM-XVCC are increasingly insensitive to the reference choice. XVCC and EOM-XVCC with the XVSCF[n] reference are found to be equivalent to the VCC formalism first developed by Prasad and coworkers,^27,32 shedding light on the relationship between the two VCC theories and XVCC theory: the VCC theory of Prasad and coworkers and XVCC theory may be diagrammatically size-consistent and infinite-basis counterparts of the VCC theory of Christiansen and coworkers in a truncated Taylor-series PES, if no other approximations are made.

We furthermore derive and implement high-order XVCCm and EOM-XVCCm equations up to m = 8 or octuple excitations, with the aid of computerized derivation, optimization, and code synthesis. We apply these methods to the zero-point energies and frequencies of the fundamentals, overtones, and combinations as well as Fermi-resonant modes of the water and formaldehyde molecules, using quartic force fields (QFF’s).⁵⁸ A comparative analysis is performed for the results of these methods and of XVSCF(n), XVSCF[n], size-extensive second-order vibrational many-body perturbation theory (XVMP2),⁵⁹ and full VCI. It shows that XVCCm and EOM-XVCCm with m equal to the truncation rank of the Taylor-series PES are uniformly accurate and stable for both resonant and nonresonant modes and are likely the methods with the greatest practical utility of all considered.

II. THEORY

A. XVCC formalism

We express the mth-order XVCC (XVCCm) wave function for the zero-point (ground) state (labeled by subscript “0”) as

| Ψ_{0} 〉 = e^{\hat{T}} | Φ_{0} 〉,

(1)

where |Φ₀〉 is a reference wave function. The latter, in turn, is a Hartree product of harmonic-oscillator modals (which include but are not limited to those from the harmonic approximation to the PES),

| Φ_{0} 〉 = \prod_{i} | 0_{i} 〉,

(2)

where |n_i〉 denotes the harmonic-oscillator wave function with quantum number n_i and frequency ω_i along the ith normal, first-order Dyson,⁵⁰ or other delocalized, rectilinear coordinate. Operator $\hat{T}$ is the sum of the single, double, through m-fold cluster excitation operators,

\hat{T} = {\hat{T}}_{1} + {\hat{T}}_{2} + \dots + {\hat{T}}_{m},

(3)

with

{\hat{T}}_{n} = \frac{1}{n!} \sum_{\underset{1}{i}, i_{2}, \dots, i_{n}} τ_{i_{1} i_{2} \dots i_{n}} {\hat{a}}_{i_{1}}^{†} {\hat{a}}_{i_{2}}^{†} \dots {\hat{a}}_{i_{n}}^{†}

(4)

= \frac{1}{n!} \sum_{\underset{1}{i}, i_{2}, \dots, i_{n}} τ_{i_{1} i_{2} \dots i_{n}} {{\hat{a}}_{i_{1}}^{†} {\hat{a}}_{i_{2}}^{†} \dots {\hat{a}}_{i_{n}}^{†}},

(5)

where τ_{i₁i₂⋯i_n} is an excitation amplitude, the value of which is to be determined by solving the amplitude equations (see below). Operator ${\hat{a}}_{i}^{†}$ and its adjoint ${\hat{a}}_{i}$ are the ladder operators for the harmonic-oscillator modals, i.e.,

{\hat{a}}_{i} | n_{i} 〉 = n_{i}^{1 / 2} | n_{i} - 1 〉,

(6)

{\hat{a}}_{i}^{†} | n_{i} 〉 = {(n_{i} + 1)}^{1 / 2} | n_{i} + 1 〉 .

(7)

The braces {⋯} in Eq. (5) bring the operators enclosed by them into a normal order.⁵² Two or more of the mode indices (i₁, i₂, …, i_n) in Eq. (4) are permitted to be coincident.

We require that the wave function of Eq. (1) satisfy the vibrational Schrödinger equation in the space of Hartree products reachable by $\hat{T}$ from |Φ₀〉,

〈 Φ_{0} | \hat{H} e^{\hat{T}} | Φ_{0} 〉 = 〈 Φ_{0} | E_{0}^{(m)} e^{\hat{T}} | Φ_{0} 〉

(8)

and

〈 Φ_{ν_{1} ν_{2} \dots ν_{n}} | \hat{H} e^{\hat{T}} | Φ_{0} 〉 = 〈 Φ_{ν_{1} ν_{2} \dots ν_{n}} | E_{0}^{(m)} e^{\hat{T}} | Φ_{0} 〉,

(9)

for all ν₁, …, ν_n and each n (1 ≤ n ≤ m). Here, $E_{0}^{(m)}$ is the XVCCm total energy and

| Φ_{ν_{1} ν_{2} \dots ν_{n}} 〉 = {\hat{a}}_{ν_{1}}^{†} {\hat{a}}_{ν_{2}}^{†} \dots {\hat{a}}_{ν_{n}}^{†} | Φ_{0} 〉 .

(10)

The number of the amplitude equations expressed collectively by Eq. (9) is equal to the number of unknown excitation amplitudes {τ} in Eq. (4). Furthermore, it can be shown below that $E_{0}^{(m)}$ ⁠, which is also unknown, is eliminated from Eq. (9). Hence, by solving Eq. (9), we obtain {τ}. Substituting the latter in the energy equation [Eq. (8)], we determine $E_{0}^{(m)}$ ⁠.

The reference wave function |Φ₀〉 can be furnished by any of the following: XVSCF(n),⁴⁹ XVSCF[n],³⁵ SCP,^53–57 or the harmonic approximation. They differ from one another in the coordinates, their center, and the modal frequencies,⁵⁰ but they all share the same form: the Hartree product of harmonic-oscillator modals in delocalized, rectilinear coordinates. The XVCCm formalisms are essentially unchanged by the choice of the references, except that some terms may vanish in some references. The foregoing definition of XVCCm and the working equations derived in the following are valid for any of them.

The rather restricted form of reference wave functions is dictated by the fact that every modal of a Hartree product in delocalized, rectilinear coordinates (such as normal coordinates) reduces to a harmonic oscillator in the thermodynamic limit, as shown both analytically⁴⁸ and numerically⁴⁹ in the case of VSCF and also in a more general sense by Makri.⁵¹ This makes it convenient to express the formalisms of vibrational correlation methods in terms of the ladder operators and force constants (as opposed to molecular integrals) and also renders the methods free of numerically defined basis functions, as discussed below. Its drawback is a lack of choices of coordinates^60–62 and the poor performance of the reference methods for strongly anharmonic potentials such as a Morse potential and a double-well potential, for overtones and combinations, and for anharmonic resonances. As shown below, some of these issues are ameliorated by XVCC, but others (such as a double-well problem) are unlikely to be resolved.

One may notice at this stage a striking difference between electronic CC and XVCC. Since the one-mode basis functions of XVCC are harmonic-oscillator wave functions, which need not be stored numerically, XVCC can take account of an infinite number of them implicitly. This means that, even when the Hamiltonian is truncated after the n-body term (n ≥ 3) for a molecule with N vibrational degrees of freedom, neither XVCCn nor XVCCN is exact mathematically, though XVCCm is shown below to be rapidly convergent with m in practice. In other words, XVCC is an infinite-basis or basis-set-free method despite the finiteness of its equations. This is where XVCC and Christiansen’s VCC (Ref. 36) also differ from each other fundamentally. The excitation space reached by $\hat{T}$ is limited by the number of basis functions included in VCC, whereas it is not in XVCC.

The VCC ansatz of Prasad and coworkers^27,32 is the same as above, but it is predicated on a reference Hartree product of Gaussian modals along normal coordinates, whose centers and widths are varied to minimize its energy. This reference method is identified as the XVSCF[n] method.³⁵ Therefore, XVCCm with the XVSCF[n] reference corresponds to the VCC theory of Prasad and coworkers,^27,32 although XVCCm is defined for more general reference Hartree products, including XVSCF(n), SCP, and the harmonic approximation.

B. Algebraic derivation of XVCC

Here, we illustrate an algebraic derivation of the working equations of XVCCm using the normal-ordered second quantization and Wick’s theorem.⁵² We use, as an example, XVCC1 or XVCC with single excitations (XVCCS) for a vibrational Hamiltonian with a QFF.

This Hamiltonian can be written in the normal-ordered second-quantized form⁵² as

\hat{H} = E_{0}^{(0)} + {\hat{H}}_{1 a} + {\hat{H}}_{1 b} + {\hat{H}}_{2 a} + {\hat{H}}_{2 b} + {\hat{H}}_{2 c} + {\hat{H}}_{3 a} + {\hat{H}}_{3 b} + {\hat{H}}_{3 c} + {\hat{H}}_{3 d} + {\hat{H}}_{4 a} + {\hat{H}}_{4 b} + {\hat{H}}_{4 c} + {\hat{H}}_{4 d} + {\hat{H}}_{4 e},

(11)

with

{\hat{H}}_{1 a} + {\hat{H}}_{1 b} = \sum_{i} W_{i} {\hat{a_{i}}} + \sum_{i} W_{i} {{\hat{a}}_{i}^{†}},

(12)

{\hat{H}}_{2 a} + {\hat{H}}_{2 b} + {\hat{H}}_{2 c} = \frac{1}{2!} \sum_{i, j} W_{i j}^{-} {{\hat{a}}_{i} {\hat{a}}_{j}} + \frac{1}{2!} \sum_{i, j} W_{i j}^{-} {{\hat{a}}_{i}^{†} {\hat{a}}_{j}^{†}} + \sum_{i, j} W_{i j}^{+} {{\hat{a}}_{i}^{†} {\hat{a}}_{j}},

(13)

{\hat{H}}_{3 a} + {\hat{H}}_{3 b} + {\hat{H}}_{3 c} + {\hat{H}}_{3 d} = \frac{1}{3!} \sum_{i, j, k} W_{i j k} {{\hat{a}}_{i} {\hat{a}}_{j} {\hat{a}}_{k}} + \frac{1}{3!} \sum_{i, j, k} W_{i j k} {{\hat{a}}_{i}^{†} {\hat{a}}_{j}^{†} {\hat{a}}_{k}^{†}} + \frac{1}{2!} \sum_{i, j, k} W_{i j k} {{\hat{a}}_{i}^{†} {\hat{a}}_{j} {\hat{a}}_{k}} + \frac{1}{2!} \sum_{i, j, k} W_{i j k} {{\hat{a}}_{i}^{†} {\hat{a}}_{j}^{†} {\hat{a}}_{k}},

(14)

{\hat{H}}_{4 a} + {\hat{H}}_{4 b} + {\hat{H}}_{4 c} + {\hat{H}}_{4 d} + {\hat{H}}_{4 e} = \frac{1}{4!} \sum_{i, j, k, l} W_{i j k l} {{\hat{a}}_{i} {\hat{a}}_{j} {\hat{a}}_{k} {\hat{a}}_{l}} + \frac{1}{4!} \sum_{i, j, k, l} W_{i j k l} {{\hat{a}}_{i}^{†} {\hat{a}}_{j}^{†} {\hat{a}}_{k}^{†} {\hat{a}}_{l}^{†}} + \frac{1}{2!} \frac{1}{2!} \sum_{i, j, k, l} W_{i j k l} {{\hat{a}}_{i}^{†} {\hat{a}}_{j}^{†} {\hat{a}}_{k} {\hat{a}}_{l}} + \frac{1}{3!} \sum_{i, j, k, l} W_{i j k l} {{\hat{a}}_{i}^{†} {\hat{a}}_{j} {\hat{a}}_{k} {\hat{a}}_{l}} + \frac{1}{3!} \sum_{i, j, k, l} W_{i j k l} {{\hat{a}}_{i}^{†} {\hat{a}}_{j}^{†} {\hat{a}}_{k}^{†} {\hat{a}}_{l}},

(15)

where it should be understood that there is term-by-term correspondence between the left- and right-hand sides of each equation.

The constant term in Eq. (11), $E_{0}^{(0)}$ ⁠, is the total ground-state energy expression of XVSCF,^35,49,52 regardless of the actual reference method adopted in the calculation, and is given by

E_{0}^{(0)} = V_{0} + \frac{1}{4} \sum_{i} ω_{i} + \frac{1}{2} \sum_{i} \frac{F_{i i}}{(2 ω_{i})} + \frac{1}{2! 2^{2}} \sum_{i, j} \frac{F_{i i j j}}{(2 ω_{i}) (2 ω_{j})} + \dots,

(16)

where V₀ is the value of the PES at the minimum, ω_i is again the frequency of the ith modal, and F’s are the force constants. The Hamiltonian amplitudes, {W}, are given by

W_{i} = \frac{F_{i}}{{(2 ω_{i})}^{1 / 2}} + \frac{1}{2} \sum_{j} \frac{F_{i j j}}{{(2 ω_{i})}^{1 / 2} (2 ω_{j})},

(17)

W_{i j}^{\pm} = \pm \frac{1}{2} δ_{i j} ω_{i} + \frac{F_{i j}}{{(2 ω_{i})}^{1 / 2} {(2 ω_{j})}^{1 / 2}} + \frac{1}{2} \sum_{k} \frac{F_{i j k k}}{{(2 ω_{i})}^{1 / 2} {(2 ω_{j})}^{1 / 2} (2 ω_{k})},

(18)

W_{i j k} = \frac{F_{i j k}}{{(2 ω_{i})}^{1 / 2} {(2 ω_{j})}^{1 / 2} {(2 ω_{k})}^{1 / 2}},

(19)

W_{i j k l} = \frac{F_{i j k l}}{{(2 ω_{i})}^{1 / 2} {(2 ω_{j})}^{1 / 2} {(2 ω_{k})}^{1 / 2} {(2 ω_{l})}^{1 / 2}},

(20)

in a QFF. These amplitudes are invariant to any permutation of indices. With the XVSCF(4) reference, $W_{i i}^{+} = ω_{i}$ holds. In the XVSCF[4] reference, additionally, W_i = 0 is true. Furthermore, $W_{i j}^{+} = 0$ for i≠j with the SCP reference. The fact that these key quantities of the XVSCF methods naturally emerge in the expression of the normal-ordered Hamiltonian underscores their fundamental significance in many-body vibrational theory.⁵²

From Eq. (8), we find the XVCCS energy equation to be

E_{0}^{(1)} = 〈 Φ_{0} | \hat{H} (1 + {\hat{T}}_{1} + \frac{1}{2!} {\hat{T}}_{1}^{2} + \frac{1}{3!} {\hat{T}}_{1}^{3} + \frac{1}{4!} {\hat{T}}_{1}^{4}) | Φ_{0} 〉,

(21)

where $E_{0}^{(1)}$ is the XVCCS total energy. The Taylor series of the exponential operator is truncated after the ${\hat{T}}_{1}^{4}$ term because the higher-order terms cannot contribute in a QFF (see below). The XVCCS amplitude equation, Eq. (9), similarly reduces to

E_{0}^{(1)} τ_{ν} = 〈 Φ_{0} | {{\hat{a}}_{ν}} \hat{H} \times (1 + {\hat{T}}_{1} + \frac{1}{2!} {\hat{T}}_{1}^{2} + \frac{1}{3!} {\hat{T}}_{1}^{3} + \frac{1}{4!} {\hat{T}}_{1}^{4} + \frac{1}{5!} {\hat{T}}_{1}^{5}) | Φ_{0} 〉,

(22)

in a QFF. Here, the Taylor series can be terminated after the ${\hat{T}}_{1}^{5}$ term (see below).

Wick’s theorem states that the vacuum expectation value of a product of two or more normal-ordered series of operators is the sum of all full contractions with no internal contractions.^2,52,63 It also states that the vacuum expectation value of just one normal-ordered series is zero.^2,52,63

Using this theorem, one can evaluate the first term of Eq. (21) as

〈 Φ_{0} | \hat{H} | Φ_{0} 〉 = E_{0}^{(0)} .

(23)

The second term has only one nonzero full contraction,

(24)

The third term has two nonzero full contractions,

(25)

The fourth and fifth terms are evaluated similarly, yielding

\frac{1}{3!} 〈 Φ_{0} | \hat{H} {\hat{T}}_{1}^{3} | Φ_{0} 〉 = \frac{1}{3!} \sum_{i, j, k} W_{i j k} τ_{i} τ_{j} τ_{k}

(26)

and

\frac{1}{4!} 〈 Φ_{0} | \hat{H} {\hat{T}}_{1}^{4} | Φ_{0} 〉 = \frac{1}{4!} \sum_{i, j, k, l} W_{i j k l} τ_{i} τ_{j} τ_{k} τ_{l} .

(27)

Since the Hamiltonian with a QFF has only up to four-body operators, the ${\hat{T}}_{1}^{5}$ term and higher in the exponential operator cannot produce a full contraction.

Together, the XVCCS energy is given by

E_{0}^{(1)} = E_{0}^{(0)} + \sum_{i} W_{i} τ_{i} + \frac{1}{2!} \sum_{i, j} W_{i j}^{-} τ_{i} τ_{j} + \frac{1}{3!} \sum_{i, j, k} W_{i j k} τ_{i} τ_{j} τ_{k} + \frac{1}{4!} \sum_{i, j, k, l} W_{i j k l} τ_{i} τ_{j} τ_{k} τ_{l} .

(28)

It may be noticed that each of the second and subsequent terms is connected in the sense that every factor has at least one common summation index with another factor. This ensures the extensivity of $E_{0}^{(1)}$ ⁠.^10,11

The connectedness of the energy equation is a universal feature of XVCCm for all m. For future convenience, therefore, we introduce an XVCCm energy equation alternative to Eq. (8), which reads

Δ E_{0}^{(m)} = 〈 Φ_{0} | {({\hat{H}}_{N} e^{\hat{T}})}_{C} | Φ_{0} 〉,

(29)

where ${\hat{H}}_{N} = \hat{H} - E_{0}^{(0)}$ ⁠, $Δ E_{0}^{(m)} = E_{0}^{(m)} - E_{0}^{(0)}$ ⁠, and (⋯)_C deletes terms that are disconnected. We call this the energy equation of XVCCm.

The XVCCS amplitude equation can likewise be derived by Wick’s theorem. The first term of Eq. (22) has one nonzero full contraction,

(30)

The second term has two nonzero full contractions,

(31)

The third term has four such contractions,

(32)

The remainders are evaluated similarly, yielding

\frac{1}{3!} 〈 Φ_{0} | {{\hat{a}}_{ν}} \hat{H} {\hat{T}}_{1}^{3} | Φ_{0} 〉 = \frac{1}{2!} \sum_{i, j} W_{i j}^{-} τ_{i} τ_{j} τ_{ν} + \frac{1}{3!} \sum_{i, j, k} W_{ν i j k} τ_{i} τ_{j} τ_{k},

(33)

\frac{1}{4!} 〈 Φ_{0} | {{\hat{a}}_{ν}} \hat{H} {\hat{T}}_{1}^{4} | Φ_{0} 〉 = \frac{1}{3!} \sum_{i, j, k} W_{i j k} τ_{i} τ_{j} τ_{k} τ_{ν},

(34)

\frac{1}{5!} 〈 Φ_{0} | {{\hat{a}}_{ν}} \hat{H} {\hat{T}}_{1}^{5} | Φ_{0} 〉 = \frac{1}{4!} \sum_{i, j, k, l} W_{i j k l} τ_{i} τ_{j} τ_{k} τ_{l} τ_{ν} .

(35)

The higher-order terms in the exponential operator cannot produce a full contraction and are zero.

Combining these, we find the XVCCS amplitude equation to read

E_{0}^{(1)} τ_{ν} = E_{0}^{(0)} τ_{ν} + \sum_{i} W_{i} τ_{i} τ_{ν} + \frac{1}{2!} \sum_{i, j} W_{i j}^{-} τ_{i} τ_{j} τ_{ν} + \frac{1}{3!} \sum_{i, j, k} W_{i j k} τ_{i} τ_{j} τ_{k} τ_{ν} + \frac{1}{4!} \sum_{i, j, k, l} W_{i j k l} τ_{i} τ_{j} τ_{k} τ_{l} τ_{ν} + W_{ν} + \sum_{i} W_{ν i}^{+} τ_{i} + \frac{1}{2!} \sum_{i, j} W_{ν i j} τ_{i} τ_{j} + \frac{1}{3!} \sum_{i, j, k} W_{ν i j k} τ_{i} τ_{j} τ_{k} .

(36)

Note that the first five terms in the right-hand side are disconnected in the sense that there is a factor (τ_ν) that has no common summation index with any other factor; each of them is, therefore, a simple product of two factors.

Factoring these disconnected terms, we can rewrite the above equation as

E_{0}^{(1)} τ_{ν} = (E_{0}^{(0)} + \sum_{i} W_{i} τ_{i} + \frac{1}{2!} \sum_{i, j} W_{i j}^{-} τ_{i} τ_{j}) + (\frac{1}{3!} \sum_{i, j, k} W_{i j k} τ_{i} τ_{j} τ_{k} + \frac{1}{4!} \sum_{i, j, k, l} W_{i j k l} τ_{i} τ_{j} τ_{k} τ_{l}) τ_{ν} + W_{ν} + \sum_{i} W_{ν i}^{+} τ_{i} + \frac{1}{2!} \sum_{i, j} W_{ν i j} τ_{i} τ_{j} + \frac{1}{3!} \sum_{i, j, k} W_{ν i j k} τ_{i} τ_{j} τ_{k}

(37)

= E_{0}^{(1)} τ_{ν} + W_{ν} + \sum_{i} W_{ν i}^{+} τ_{i} + \frac{1}{2!} \sum_{i, j} W_{ν i j} τ_{i} τ_{j} + \frac{1}{3!} \sum_{i, j, k} W_{ν i j k} τ_{i} τ_{j} τ_{k},

(38)

where we use the fact that the sum in the parenthesis is equal to $E_{0}^{(1)}$ as per Eq. (28). We can, therefore, further simplify the amplitude equation into

T_{1}^{(1)} \equiv W_{ν} + \sum_{i} W_{ν i}^{+} τ_{i} + \frac{1}{2!} \sum_{i, j} W_{ν i j} τ_{i} τ_{j} + \frac{1}{3!} \sum_{i, j, k} W_{ν i j k} τ_{i} τ_{j} τ_{k} = 0 .

(39)

Notice that this equation consists of connected terms only.

In fact, we can recast the amplitude equations of XVCCm in general into a form that is equivalent to, but more convenient than Eq. (9),

〈 Φ_{ν_{1} ν_{2} \dots ν_{n}} | {({\hat{H}}_{N} e^{\hat{T}})}_{C} | Φ_{0} 〉 = 0,

(40)

for all ν₁, …, ν_n and all n (1 ≤ n ≤ m). We call this the T_n amplitude equation of XVCCm. The diagrammatic derivations of XVCC and the EOM-XVCC formulations to be discussed below are based on Eqs. (29) and (40).

In the foregoing derivation, the connected structure of the XVCC energy and amplitude equations is only illustrated with the XVCCS example. In Appendix A, we document an alternative derivation that proves more rigorously the connectedness of these equations at any order.^1,2

C. Diagrammatic derivation of XVCC

The equations of XVCC can be derived even more expediently with Feynman–Goldstone diagrams, which is equivalent to the normal-ordered second-quantized derivation; a diagram corresponds to all topologically related Wick’s contractions. Here, we introduce the rules of drawing closed diagrams for XVCC energy equations and open diagrams for XVCC amplitude equations as well as the rules of translating them into algebraic expressions. These rules differ from those for XVMP energies and self-energies.⁶⁴ We illustrate their application to XVCCm with m = 1 (XVCCS) and m = 2 (XVCCSD) in this subsection.

We recast each term in the normal-ordered Hamiltonian, ${\hat{H}}_{N}$ ⁠, into diagrammatic forms, as shown in Fig. 1. Each diagram represents one of the normal-ordered operators defined by Eqs. (12)–(15). The cluster excitation operators, ${{\hat{T}}_{n}}$ ⁠, are expressed diagrammatically in Fig. 2.

FIG. 1.

View large Download slide

Diagrammatic representation of the normal-ordered Hamiltonian ${\hat{H}}_{N}$ with a QFF.

FIG. 2.

View large Download slide

Diagrammatic representation of the cluster excitation operators.

Table I lists the rules to generate diagrams of the XVCCm energy and amplitude equations in an nth-order Taylor-series PES. Applying these rules, we find the energy diagrams of XVCCS in a QFF (n = 4) to be those in Fig. 3. The energy diagrams of XVCCSD in a QFF are the union of Figs. 3 and 4. Hence,

Δ E_{0}^{(1)} = \sum_{i = 1}^{4} E_{i},

(41)

Δ E_{0}^{(2)} = \sum_{i = 1}^{8} E_{i},

(42)

where the summands are the algebraic interpretation of the diagrams. Since each ${\hat{T}}_{k}$ vertex in an energy diagram has only upgoing lines (signifying the fact that they describe an excitation process), all of its lines have no choice but to connect with the ${\hat{H}}_{N}$ vertex, making the whole diagram always connected when closed, guaranteeing its thermodynamic extensivity as per the extensive diagram theorem.^10,11,65

TABLE I.

Rules to draw all closed diagrams in the energy equation and all open diagrams in the T_l amplitude equation of XVCCm in an nth-order Taylor-series PES.

(1)

Draw an

{\hat{H}}_{N}

vertex having up to n lines.

(2)

Draw 0 to

n {\hat{T}}_{k}

vertices (1 ≤ k ≤ m) beneath the

{\hat{H}}_{N}

vertex.

(3)

Connect the lines of the

{\hat{T}}_{k}

vertices with the lines of the

{\hat{H}}_{N}

vertex without altering (upgoing or downgoing) line directions,

leaving l external (i.e., dangling) lines to form a connected

diagram, which is either closed (l = 0) or open (l > 0).

(4)

Repeat steps (1)–(3) to enumerate all topologically distinct

diagrams. If permutation of

\hat{T}

vertices allows two diagrams

to overlap, they are deemed topologically equivalent.

FIG. 3.

View large Download slide

Diagrams in the energy equation of XVCCS in a QFF.

FIG. 4.

View large Download slide

Additional diagrams in the energy equation of XVCCSD in a QFF.

Table II provides a set of rules to translate these diagrams into algebraic expressions. For example, diagram E₄ of Fig. 3 can be labeled as shown in Fig. 5 using rules 1 through 5 of the table. Following the remaining rules, we find the algebraic expression of the whole diagram to be

E_{4} \equiv \frac{1}{4!} \sum_{i, j, k, l} W_{i j k l} τ_{i} τ_{j} τ_{k} τ_{l},

(43)

where the 1/4! factor comes from rule 8 because there are four equivalent ${\hat{T}}_{1}$ vertices. It is straightforward to confirm that the algebraic interpretations of E₁, E₂, and E₃ obtained with Table II agree exactly with the second, third, and fourth terms of the XVCCS energy expression [Eq. (28)], respectively.

TABLE II.

Rules to translate a closed or open diagram of XVCC and EOM-XVCC into an algebraic expression.

(1)

Label each internal line with a mode index (i, j, k, …).

An internal line is a line terminated by vertices at both ends.

(2)

Label each external line with a mode index (ν, μ, λ, …).

An external line is a line terminated by a vertex at only one end.

(3)

Associate each

{\hat{H}}_{N}

vertex with a Hamiltonian amplitude W.

The indices of W are given by the labels of the lines connected

to that vertex. A two-line

{\hat{H}}_{N}

vertex is associated with W⁻

if both lines are upgoing or downgoing; with W⁺ otherwise.

(4)

Associate each

{\hat{T}}_{k}

vertex with an excitation amplitude τ.

The indices of τ are given by the labels of the lines connected

to that vertex.

(5)

Associate each

{\hat{R}}_{k}

vertex with an excitation amplitude r.

The indices of r are given by the labels of the lines connected

to that vertex.

(6)

Sum over all internal line indices (i, j, k, …).

(7)

Sum over all permutations of external line indices (ν, μ, λ, …).

(8)

For each set of n equivalent vertices, multiply by a factor of

(1/n!). If permutation of two vertices leaves the diagram

unchanged topologically, the two vertices are equivalent.

(9)

For each set of m equivalent lines, multiply by a factor of

(1/m!). If permutation of two lines leaves the diagram

unchanged topologically, the two lines are equivalent.

FIG. 5.

View large Download slide

Diagram E₄ with labels.

Likewise, the application of the rules in Table II to the diagrams in Fig. 4 leads to the XVCCSD energy equation,

Δ E_{0}^{(2)} = Δ E_{0}^{(1)} + E_{5} + E_{6} + E_{7} + E_{8}

(44)

= Δ E_{0}^{(1)} + \frac{1}{2!} \sum_{i, j} W_{i j}^{-} τ_{i j} + \frac{1}{2!} \sum_{i, j, k} W_{i j k} τ_{i} τ_{j k} + \frac{1}{2!} \frac{1}{2!} \frac{1}{2!} \sum_{i, j, k, l} W_{i j k l} τ_{i j} τ_{k l} + \frac{1}{2!} \frac{1}{2!} \sum_{i, j, k, l} W_{i j k l} τ_{i} τ_{j} τ_{k l} .

(45)

The factor of (1/2!)³ in the penultimate term (diagram E₇) arises from the fact that there are one pair of equivalent ${\hat{T}}_{2}$ vertices and two pairs of equivalent lines.

Table I also outlines rules to enumerate all diagrams in the XVCCm amplitude equations. Unlike energy diagrams, a diagram in the T_l amplitude equation is open with l external lines. Such diagrams can be disconnected, in which case they must be deleted as per rule 3.

Figure 6 lists all four connected diagrams of the XVCCS T₁ amplitude equation in a QFF, which can thus be written as

0 = \sum_{i = 1}^{4} S_{i} \equiv T_{1}^{(1)},

(46)

where S_i is the algebraic interpretation of the corresponding diagram. The T₁ amplitude equation of XVCCSD in a QFF is diagrammatically expressed by Figs. 6 and 7. The T₂ amplitude equation of the same method is defined with the diagrams in Fig. 8. The T₁ and T₂ amplitude equations of XVCCSD are thus formally written as

0 = \sum_{i = 1}^{12} S_{i} \equiv T_{1}^{(2)},

(47)

0 = \sum_{i = 1}^{12} D_{i} \equiv T_{2}^{(2)} .

(48)

FIG. 6.

View large Download slide

Diagrams in the T₁ amplitude equation of XVCCS in a QFF.

FIG. 7.

View large Download slide

Additional diagrams in the T₁ amplitude equation of XVCCSD in a QFF.

FIG. 8.

View large Download slide

Diagrams in the T₂ amplitude equation of XVCCSD in a QFF.

The interpretation rules in Table II apply to open diagrams in these amplitude equations also. For example, lines and vertices in diagram S₄ can be labeled as shown in Fig. 9. Its algebraic expression then reads

S_{4} = \frac{1}{3!} \sum_{i, j, k} W_{ν j k l} τ_{i} τ_{j} τ_{k} .

(49)

Applying the same rules to diagrams S₁, S₂, and S₃ and substituting the results in Eq. (46) reproduces the T₁ amplitude equation [Eq. (39)] derived algebraically.

FIG. 9.

View large Download slide

Diagram S₄ with labels.

The T₁ amplitude equation of XVCCSD in a QFF is likewise obtained directly from Figs. 6 and 7,

0 = T_{1}^{(2)}

(50)

= T_{1}^{(1)} + \sum_{i} W_{i} τ_{i ν} + \frac{1}{2!} \sum_{i, j} W_{ν i j} τ_{i j} + \sum_{i, j} W_{i j}^{-} τ_{i} τ_{j ν} + \frac{1}{2!} \sum_{i, j, k} W_{ν i j k} τ_{i} τ_{j k} + \frac{1}{2!} \sum_{i, j, k} W_{i j k} τ_{i j} τ_{k ν} + \frac{1}{2!} \sum_{i, j, k} W_{i j k} τ_{i} τ_{j} τ_{k ν} + \frac{1}{2!} \sum_{i, j, k, l} W_{i j k l} τ_{i} τ_{j k} τ_{l ν} + \frac{1}{3!} \sum_{i, j, k, l} W_{i j k l} τ_{i} τ_{j} τ_{k} τ_{l ν},

(51)

where $T_{1}^{(1)}$ is given algebraically in Eq. (39).

Similarly, the T₂ amplitude equation of the same method can directly result from Fig. 8 and reads

0 = T_{2}^{(2)}

(52)

= \frac{1}{2!} \hat{P} (ν μ) W_{ν μ}^{-} + \frac{1}{2!} \hat{P} (ν μ) \sum_{i} W_{ν μ i} τ_{i} + \hat{P} (ν μ) \sum_{i} W_{ν i}^{+} τ_{i μ} + \frac{1}{2!} \frac{1}{2!} \hat{P} (ν μ) \sum_{i, j} W_{ν μ i j} τ_{i j} + \frac{1}{2!} \frac{1}{2!} \hat{P} (ν μ) \sum_{i, j} W_{ν μ i j} τ_{i} τ_{j} + \hat{P} (ν μ) \sum_{i, j} W_{ν i j} τ_{i} τ_{j μ} + \frac{1}{2!} \hat{P} (ν μ) \sum_{i, j} W_{i j}^{-} τ_{i ν} τ_{j μ} + \frac{1}{2!} \hat{P} (ν μ) \sum_{i, j, k} W_{ν i j k} τ_{i j} τ_{k μ} + \frac{1}{2!} \hat{P} (ν μ) \sum_{i, j, k} W_{ν i j k} τ_{i} τ_{j} τ_{k μ} + \frac{1}{2!} \hat{P} (ν μ) \sum_{i, j, k} W_{i j k} τ_{i} τ_{j ν} τ_{k μ} + \frac{1}{2!} \frac{1}{2!} \hat{P} (ν μ) \sum_{i, j, k, l} W_{i j k l} τ_{i j} τ_{k ν} τ_{l μ} + \frac{1}{2!} \frac{1}{2!} \hat{P} (ν μ) \sum_{i, j, k, l} W_{i j k l} τ_{i} τ_{j} τ_{k ν} τ_{l μ},

(53)

where $\hat{P} (ν_{1} \dots ν_{n})$ yields the sum of n! permutations of indices {ν₁, …, ν_n} in the term it acts on (cf. rule 7 of Table II). For instance,

\hat{P} (ν μ) \sum_{i} W_{ν i}^{+} τ_{i μ} = \sum_{i} W_{ν i}^{+} τ_{i μ} + \sum_{i} W_{μ i}^{+} τ_{i ν} .

(54)

In the same way, the diagrammatic and algebraic energy and amplitude equations of XVCCm with m = 3 (XVCCSDT) in a QFF are derived. They are documented in Appendix B. For 4 ≤ m ≤ 8, diagrams and equations were generated with the aid of computer algebra based on the rules provided in Tables I and II and are documented in the supplementary material.⁶⁶

D. EOM-XVCC formalism

In a vibrational analysis, it is usually the frequencies (excitation energies) that are of greater interest than the zero-point energies (total energies in the ground state). Here, we stipulate and implement an excited-state XVCC theory using the EOM formalism,^67–72 which is equivalent to CC linear response^73–78 for excitation energies.

The EOM-XVCCm theory can be viewed as a CI for excited states using the XVCC effective Hamiltonian ${\bar{H}}_{N}$ defined as

{\bar{H}}_{N} = {({\hat{H}}_{N} e^{\hat{T}})}_{C},

(55)

where $\hat{T}$ is truncated after the m-fold cluster excitation operator [Eq. (3)] and all of its amplitudes should have been determined by a preceding XVCCm calculation.

Let $\hat{R}$ be a linear excitation operator truncated after the same rank,

\hat{R} = r_{0} + {\hat{R}}_{1} + {\hat{R}}_{2} + \dots + {\hat{R}}_{m},

(56)

with

{\hat{R}}_{n} = \frac{1}{n!} \sum_{\underset{1}{i}, i_{2}, \dots, i_{n}} r_{i_{1} i_{2} \dots i_{n}} {\hat{a}}_{i_{1}}^{†} {\hat{a}}_{i_{2}}^{†} \dots {\hat{a}}_{i_{n}}^{†}

(57)

= \frac{1}{n!} \sum_{\underset{1}{i}, i_{2}, \dots, i_{n}} r_{i_{1} i_{2} \dots i_{n}} {{\hat{a}}_{i_{1}}^{†} {\hat{a}}_{i_{2}}^{†} \dots {\hat{a}}_{i_{n}}^{†}},

(58)

where r_{i₁i₂⋯i_n} is an excitation amplitude. We require that the excited-state wave function $\hat{R} | Φ_{0} 〉$ satisfy the Schrödinger equation with the modified Hamiltonian ${\bar{H}}_{N}$ in the space of Hartree products reachable by $\hat{R}$ from |Φ₀〉. That is,

〈 Φ_{0} | {\bar{H}}_{N} \hat{R} | Φ_{0} 〉 = 〈 Φ_{0} | Δ E_{k}^{(m)} \hat{R} | Φ_{0} 〉,

(59)

〈 Φ_{ν_{1} ν_{2} \dots ν_{n}} | {\bar{H}}_{N} \hat{R} | Φ_{0} 〉 = 〈 Φ_{ν_{1} ν_{2} \dots ν_{n}} | Δ E_{k}^{(m)} \hat{R} | Φ_{0} 〉,

(60)

for all ν₁, …, ν_n and all n (1 ≤ n ≤ m), where $Δ E_{k}^{(m)} = E_{k}^{(m)} - E_{0}^{(0)}$ and $E_{k}^{(m)}$ is the total EOM-XVCCm energy of the kth excited state. The constant term r₀ is formally necessary in Eq. (56) for excited states of the same symmetry as the ground state. Although $\hat{R}$ and its amplitudes vary with k and should thus carry a k identifier, it is suppressed for brevity.

Equations (59) and (60) define EOM-XVCCm as an eigenvalue problem with a non-Hermitian matrix of ${\bar{H}}_{N}$ ⁠, which thus has both left and right eigenvectors with $\hat{R}$ corresponding to the latter. They can, however, be further simplified by using the fact that T amplitudes in ${\bar{H}}_{N}$ satisfy the XVCCm equations. Equations (29) and (40) indicate that |Φ₀〉 satisfies the Schrödinger equation with ${\bar{H}}_{N}$ in the same Hartree-product space. This implies

〈 Φ_{0} | \hat{R} {\bar{H}}_{N} | Φ_{0} 〉 = 〈 Φ_{0} | \hat{R} Δ E_{0}^{(m)} | Φ_{0} 〉,

(61)

〈 Φ_{ν_{1} ν_{2} \dots ν_{n}} | \hat{R} {\bar{H}}_{N} | Φ_{0} 〉 = 〈 Φ_{ν_{1} ν_{2} \dots ν_{n}} | \hat{R} Δ E_{0}^{(m)} | Φ_{0} 〉,

(62)

because ${\hat{R}}^{†} | Φ_{0} 〉 \cup {{\hat{R}}^{†} | Φ_{ν_{1} ν_{2} \dots ν_{n}} 〉}$ spans the space no greater than |Φ₀〉∪{|Φ_{ν₁ν₂⋯ν_n}〉}. Subtracting these from Eqs. (59) and (60), we obtain

〈 Φ_{0} | [{\bar{H}}_{N}, \hat{R}] | Φ_{0} 〉 = ω_{k}^{(m)} r_{0},

(63)

〈 Φ_{ν_{1} ν_{2} \dots ν_{n}} | [{\bar{H}}_{N}, \hat{R}] | Φ_{0} 〉 = ω_{k}^{(m)} r_{ν_{1} ν_{2} \dots ν_{n}},

(64)

where $ω_{k}^{(m)} = Δ E_{k}^{(m)} - Δ E_{0}^{(m)}$ and is the EOM-XVCCm transition energy (frequency) to the kth excited state.

Equation (64) is the ansatz for EOM-XVCCm for frequencies, which can also be written as

〈 Φ_{ν_{1} ν_{2} \dots ν_{n}} | {({({\hat{H}}_{N} e^{\hat{T}})}_{C} \hat{R})}_{C} | Φ_{0} 〉 = ω_{k}^{(m)} r_{ν_{1} ν_{2} \dots ν_{n}},

(65)

for all ν₁, …, ν_n and all n (1 ≤ n ≤ m). We call this the R_n amplitude equation of EOM-XVCCm. The nested connectedness requirement is always met when the term is simply connected because both $\hat{T}$ and $\hat{R}$ are excitation operators and can only directly connect with ${\hat{H}}_{N}$ ⁠, but not with each other. The connectedness requirement, however, deletes all terms involving r₀ because the latter has no index (or a line in its diagrammatic representation) to be shared with other factors.

Equation (65), therefore, refers to a set of equations whose number is equal to the number of unknowns (excluding r₀) minus one (because $ω_{k}^{(m)}$ is also unknown) and, with an appropriate normalization condition, can determine $ω_{k}^{(m)}$ and all the amplitudes in $\hat{R}$ except r₀. The value of r₀, in turn, can be determined by Eq. (63), though it is not necessary if only frequencies are sought. Equation (65) satisfies the intensive diagram theorem,^10,11 ensuring the size consistency of EOM-XVCCm for total excited-state energies and frequencies, despite its CI-like structure.

Just like the electronic counterpart, EOM-XVCCm can also be derived by applying the first-order time-dependent perturbation (linear-response) theory to XVCCm. This derivation (not repeated here because it is identical to the electronic case⁷⁸) elegantly explains why the truncation ranks of $\hat{T}$ and $\hat{R}$ should be the same. EOM-XVCCm with the XVSCF[n] reference is the same as excited-state VCC of Prasad and coworkers.^27,32

E. Algebraic derivation of EOM-XVCC

Here, we illustrate the normal-ordered second-quantized derivation of the amplitude equations of EOM-XVCCm using m = 1 (EOM-XVCCS) in a QFF as an example. The working equations of EOM-XVCCm with m = 2 (EOM-XVCCSD) are documented in Appendix C.

The left-hand side of Eq. (65) in the case of m = 1 reads

〈 Φ_{0} | {{\hat{a}}_{ν}} {({({\hat{H}}_{N} e^{{\hat{T}}_{1}})}_{C} {\hat{R}}_{1})}_{C} | Φ_{0} 〉 = 〈 Φ_{0} | {{\hat{a}}_{ν}} {({\hat{H}}_{N} {\hat{R}}_{1})}_{C} | Φ_{0} 〉 + 〈 Φ_{0} | {{\hat{a}}_{ν}} {({\hat{H}}_{N} {\hat{T}}_{1} {\hat{R}}_{1})}_{C} | Φ_{0} 〉 + \frac{1}{2!} 〈 Φ_{0} | {{\hat{a}}_{ν}} {({\hat{H}}_{N} {\hat{T}}_{1}^{2} {\hat{R}}_{1})}_{C} | Φ_{0} 〉,

(66)

where the nested connectedness means simple connectedness in this case and terms containing higher powers of ${\hat{T}}_{1}$ cannot yield any connected contributions.

The first term of Eq. (66) has only one full contraction,

(67)

which is connected. The next term has two connected full contractions,

(68)

The last term yields the connected full contraction,

\frac{1}{2!} 〈 Φ_{0} | {{\hat{a}}_{ν}} {({\hat{H}}_{N} {\hat{T}}_{1}^{2} {\hat{R}}_{1})}_{C} | Φ_{0} 〉 = \frac{1}{2!} \sum_{i, j, k} W_{ν i j k} τ_{i} τ_{j} r_{k} .

(69)

Together, we obtain the R₁ amplitude equation of EOM-XVCCS in a QFF as

R_{1}^{(1)} \equiv \sum_{i} W_{ν i}^{+} r_{i} + \sum_{i, j} W_{ν i j} τ_{i} r_{j} + \frac{1}{2!} \sum_{i, j, k} W_{ν i j k} τ_{i} τ_{j} r_{k} = ω_{k}^{(1)} r_{ν} .

(70)

F. Diagrammatic derivation of EOM-XVCC

The diagrammatic forms of the ${\hat{R}}_{n}$ operators are given in Fig. 10. They are isomorphic to the diagrams of the ${\hat{T}}_{n}$ operators with the only difference being the colors of the vertices. The rules to generate all diagrams entering the amplitude equations of EOM-XVCCm are given in Table III. The rules for interpreting these diagrams algebraically are provided in Table II.

FIG. 10.

View large Download slide

Diagrammatic representation of the linear excitation operators.

TABLE III.

Rules to draw all open diagrams in the R_l amplitude equation of EOM-XVCCm in an nth-order Taylor-series PES.

(1)

Draw an

{\hat{H}}_{N}

vertex having up to n lines.

(2)

Draw a

{\hat{R}}_{h}

vertex (1 ≤ h ≤ m) and zero to

(n - 1) {\hat{T}}_{k}

vertices

(1 ≤ k ≤ m) beneath the

{\hat{H}}_{N}

vertex.

(3)

Connect the lines of the

{\hat{R}}_{h}

and

{\hat{T}}_{k}

vertices with the lines of

the

{\hat{H}}_{N}

vertex without altering (upgoing or downgoing) line

directions, leaving l external (i.e., dangling) lines to form

a connected diagram.

(4)

Repeat steps (1)–(3) to enumerate all topologically distinct

diagrams. If permutation of

\hat{T}

vertices allows two

diagrams to overlap, they are deemed topologically equivalent.

As an illustration, we derive the R₁ amplitude equation of EOM-XVCCS in a QFF diagrammatically. Diagrams in this equation are drawn in Fig. 11. The R₁ amplitude equation can thus be formally written as

R_{1}^{(1)} \equiv S_{1}^{'} + S_{2}^{'} + S_{3}^{'} = ω_{k}^{(1)} r_{ν} .

(71)

Diagram $S_{3}^{'}$ may be labeled as shown in Fig. 12 according to the rules in Table II. The corresponding algebraic translation is

S_{3}^{'} \equiv \frac{1}{2!} \sum_{i, j, k} W_{ν i j k} τ_{j} τ_{k} r_{i},

(72)

where the 1/2! factor is due to a pair of equivalent ${\hat{T}}_{1}$ vertices. The other two diagrams can be interpreted similarly, yielding the same R₁ amplitude equation as Eq. (70).

FIG. 11.

View large Download slide

Diagrams in the R₁ amplitude equation of EOM-XVCCS in a QFF.

FIG. 12.

View large Download slide

Diagram $S_{3}^{'}$ with labels.

The diagrammatic and algebraic equations of EOM- XVCCm for m = 2 through 8 were generated with the aid of computer algebra based on the rules provided in Tables II and III. The EOM-XVCCSD equations are documented in Appendix C, while the EOM-XVCCm equations for 3 ≤ m ≤ 8 are documented in the supplementary material.⁶⁶

III. AUTOMATED COMPUTER IMPLEMENTATION

The XVCCm and EOM-XVCCm methods with m = 1 through 8 in a QFF were implemented. The Hamiltonian employed was the pure vibrational one, neglecting any coupling between vibrations and rotations. The XVSCF(4), XVSCF[4], or the harmonic reference wave function was used (“4” standing for the QFF). The choice of the reference wave function is indicated by a suffix as in XVCCm (4), XVCCm[4], or XVCCm{4}, respectively; “{n}” suffix indicates the use of a harmonic reference with the nth-order PES.

The T_n amplitude equations of XVCCm were solved with a standard Newton–Krylov algorithm,⁷⁹ implemented in igmres, available from scipy. The initial values of τ were zero. The R_n amplitude equations of EOM-XVCCm with m ≥ 2 were solved with the iterative Arnoldi method, a non-Hermitian generalization of the Lanczos method,^80–82 available in arpack.⁸³ The eigenvalue convergence criterion was set to machine precision. The Arnoldi method does not require the matrix elements of ${\bar{H}}_{N}$ stored in memory; instead, it evaluates the left-hand sides of the R_n amplitude equations [Eq. (65)] directly for several trial r vectors. This allows a small subset of low-lying excited-state roots to be determined in a memory-efficient manner. On the other hand, the EOM-XVCCS problem was solved by diagonalization of the ${\bar{H}}_{N}$ matrix, solving Eqs. (59) and (60) for all excited-state energies instead.

The XVCCm and EOM-XVCCm diagrammatic equations were derived by a symbolic algebra program using the rules presented in Tables I and III. These diagrammatic equations were then translated by the same symbolic algebra program into the algebraic form according to the interpretation rules given in Table II. Terms that originate from a single diagram and are thus related by a permutation operator were identified and consolidated. When the resulting expressions were evaluated literally, the operation costs of XVCCm and EOM-XVCCm in an nth-order Taylor-series PES would scale as steeply as O(N^m+n), where N is the number of modes. To reduce the cost, three types of algorithm optimizations were performed computationally before the automatic synthesis of c++ codes was carried out.^84,85 The optimizations considered in this work were the strength reduction, intermediate reuse, and factorization.⁸⁴ The details of each optimization procedure are given below.

Every tensor has a high degree of index permutation symmetry such as τ_i₁i₂i₃ = τ_i₁i₃i₂ = τ_i₂i₁i₃ = τ_i₂i₃i₁ = τ_i₃i₁i₂ = τ_i₃i₂i₁. Exploiting it can compress the memory spaces of the m-fold cluster and linear excitation amplitudes and of mth-order force constants by a factor of m! and the operation cost of their contractions to an even greater degree. Nevertheless, the compression factor of m! in XVCCm is much smaller than (m!)² in electronic CC theory with connected m-fold excitations, and our preliminary implementation of the algorithm that takes into account this symmetry did not demonstrate any gain in computational speed owing to an enormous cost overhead associated with complicated data layout. Therefore, in this study, we neglected this symmetry.

A. Strength reduction

The algebraic equations of XVCC and EOM-XVCC are sums of contractions of tensors (force constants W, T amplitudes, and R amplitudes). The operation cost of evaluating a contraction (product) of more than two tensors is greatly reduced by carrying out a series of binary tensor contractions. Suppose a tensor contraction E = ABCD, where A through E are tensors of some dimension (including possibly scalars). This product can be evaluated stepwise as

I_{1} = C D, I_{2} = B I_{1}, E = A I_{2},

(73)

where I₁ and I₂ are the so-called intermediate tensors. Alternatively, the same product may be evaluated as

I_{3} = A B, I_{4} = C I_{3}, E = D I_{4} .

(74)

Both orders of contraction give the same final result (E) but at generally different operation and memory costs. Strength reduction determines the sequential contraction order with the minimal operation cost for each product by blanket search; it neglects parallel contraction such as (AB)(CD).

The result of this procedure is depicted in Fig. 13 for diagram T₁₇, which occurs in the T₃ amplitude equation of XVCCSDT (see Appendix B). It is algebraically defined as

{(T_{17})}_{ν μ λ} = \frac{1}{2!} \frac{1}{2!} \hat{P} (ν μ λ) \sum_{i, j, k, l} W_{i j k l} τ_{l} τ_{j k ν} τ_{i μ λ} .

(75)

A literal evaluation of this equation involves an O(N⁷) arithmetic operations, as it must perform N⁴ multiplications and additions to evaluate each of N³ elements of T₁₇, where N is the number of modes. However, the same result can be obtained at an O(N⁴) cost by breaking the quaternary tensor contraction into three separate binary tensor contractions carried out in the following order. First, the force-constant tensor is multiplied with τ_l, which has no external index,

{(T_{17 a})}_{i j k} = \sum_{l} W_{i j k l} τ_{l} .

(76)

This step involves N⁴ multiplications and additions that define an intermediate, T_17a, of dimension N³. This intermediate is then multiplied by τ_jkν at an O(N⁴) cost to define another intermediate, T_17b, of dimension N²,

{(T_{17 b})}_{ν i} = \sum_{j, k} {(T_{17 a})}_{i j k} τ_{j k ν} .

(77)

Finally, T₁₇ is obtained by contraction of intermediate T_17b with τ_iμλ at an O(N⁴) cost,

{(T_{17})}_{ν μ λ} = \frac{1}{2!} \frac{1}{2!} \hat{P} (ν μ λ) \sum_{i} {(T_{17 b})}_{ν i} τ_{i μ λ} .

(78)

Therefore, the cost of evaluating diagram T₁₇ is reduced from O(N⁷) to O(N⁴).

FIG. 13.

View large Download slide

Strength reduction of diagram T₁₇ in the T₃ amplitude equation of XVCCSDT (see Appendix B).

B. Intermediate reuse

This optimization simply stores an intermediate appearing two or more times in the whole equation and reuses it without recalculating. Consider the following sum-of-product expression of tensors: E = A(BC) + D(BC). If the strength reduction procedure suggests the optimal contraction orders specified by the parentheses, intermediate I₁ = BC is calculated only once and stored for reuse. Hence, the whole sum is evaluated as

I_{1} = B C, E = A I_{1} + D I_{1} .

(79)

The numbers of multiplications and additions as a function of N required to evaluate all the amplitude equations of XVCCm and EOM-XVCCm with no optimization and with strength reduction and intermediate reuse are compared in Table IV. For all methods with m ≥ 2, strength reduction leads to a great reduction in the operation cost from O(N^m+4) to O(N^m+2). Generally, for an nth-order PES, the optimal scaling of XVCCm and EOM-XVCCm is O(N^m+⌊n/2⌋) when m ≥ n and is no greater than O(N^{n+⌊(m−1)/2⌋}) when m < n. A proof of this statement is given in Appendix D.

TABLE IV.

The number of multiplications and additions as a function of the number of modes (N) required in evaluating all the amplitude equations of XVCC or EOM-XVCC implemented with no optimization, with strength reduction and intermediate reuse, or with strength reduction, intermediate reuse, and factorization. The three leading terms only are shown.

Method	No optimization	Strength reduction	Factorization^a
XVCCS	5N⁴ + 4N³ + 2N²	2N⁴ + 4N³ + 8N²	2N⁴ + 6N³ + 4N²
XVCCSD	11N⁶ + 31N⁵ + 34N⁴	13N⁴ + 54N³ + 36N²	13N⁴ + 44N³ + 25N²
XVCCSDT	60N⁷ + 119N⁶ + 110N⁵	11N⁵ + 165N⁴ + 96N³	11N⁵ + 86N⁴ + 70N³
XVCCSDTQ	241N⁸ + 360N⁷ + 282N⁶	24N⁶ + 425N⁵ + 225N⁴	24N⁶ + 172N⁵ + 109N⁴
XVCC5	799N⁹ + 1053N⁸ + 720N⁷	104N⁷ + 1061N⁶ + 488N⁵	82N⁷ + 342N⁶ + 194N⁵
XVCC6	2577N¹⁰ + 3086N⁹ + 1888N⁸	175N⁸ + 3079N⁷ + 1134N⁶	143N⁸ + 1135N⁷ + 366N⁶
XVCC7	8398N¹¹ + 9246N¹⁰ + 5070N⁹	285N⁹ + 8902N⁸ + 3156N⁷	241N⁹ + 4392N⁸ + 1161N⁷
XVCC8	28 099N¹² + 28 380N¹¹ + 13 989N¹⁰	468N¹⁰ + 26 641N⁹ + 8991N⁸	410N¹⁰ + 18 021N⁹ + 4420N⁸
EOM-XVCCS	5N⁴ + 3N³ + 2N²	2N⁴ + 4N³ + 7N²	2N⁴ + 4N³ + 7N²
EOM-XVCCSD	30N⁶ + 70N⁵ + 54N⁴	21N⁴ + 97N³ + 63N²	21N⁴ + 57N³ + 52N²
EOM-XVCCSDT	183N⁷ + 286N⁶ + 220N⁵	15N⁵ + 343N⁴ + 189N³	15N⁵ + 131N⁴ + 99N³
EOM-XVCCSDTQ	729N⁸ + 916N⁷ + 598N⁶	34N⁶ + 965N⁵ + 474N⁴	34N⁶ + 311N⁵ + 170N⁴
EOM-XVCC5	2490N⁹ + 2794N⁸ + 1614N⁷	207N⁷ + 2649N⁶ + 1096N⁵	143N⁷ + 696N⁶ + 356N⁵
EOM-XVCC6	8295N¹⁰ + 8533N⁹ + 4432N⁸	348N⁸ + 8164N⁷ + 2782N⁶	254N⁸ + 2527N⁷ + 731N⁶
EOM-XVCC7	27 975N¹¹ + 26 506N¹⁰ + 12 457N⁹	569N⁹ + 24 962N⁸ + 8309N⁷	439N⁹ + 11 231N⁸ + 2570N⁷
EOM-XVCC8	96 603N¹² + 83 992N¹¹ + 35 934N¹⁰	934N¹⁰ + 78 522N⁹ + 25 129N⁸	762N¹⁰ + 51 289N⁹ + 11 274N⁸

Method	No optimization	Strength reduction	Factorization^a
XVCCS	5N⁴ + 4N³ + 2N²	2N⁴ + 4N³ + 8N²	2N⁴ + 6N³ + 4N²
XVCCSD	11N⁶ + 31N⁵ + 34N⁴	13N⁴ + 54N³ + 36N²	13N⁴ + 44N³ + 25N²
XVCCSDT	60N⁷ + 119N⁶ + 110N⁵	11N⁵ + 165N⁴ + 96N³	11N⁵ + 86N⁴ + 70N³
XVCCSDTQ	241N⁸ + 360N⁷ + 282N⁶	24N⁶ + 425N⁵ + 225N⁴	24N⁶ + 172N⁵ + 109N⁴
XVCC5	799N⁹ + 1053N⁸ + 720N⁷	104N⁷ + 1061N⁶ + 488N⁵	82N⁷ + 342N⁶ + 194N⁵
XVCC6	2577N¹⁰ + 3086N⁹ + 1888N⁸	175N⁸ + 3079N⁷ + 1134N⁶	143N⁸ + 1135N⁷ + 366N⁶
XVCC7	8398N¹¹ + 9246N¹⁰ + 5070N⁹	285N⁹ + 8902N⁸ + 3156N⁷	241N⁹ + 4392N⁸ + 1161N⁷
XVCC8	28 099N¹² + 28 380N¹¹ + 13 989N¹⁰	468N¹⁰ + 26 641N⁹ + 8991N⁸	410N¹⁰ + 18 021N⁹ + 4420N⁸
EOM-XVCCS	5N⁴ + 3N³ + 2N²	2N⁴ + 4N³ + 7N²	2N⁴ + 4N³ + 7N²
EOM-XVCCSD	30N⁶ + 70N⁵ + 54N⁴	21N⁴ + 97N³ + 63N²	21N⁴ + 57N³ + 52N²
EOM-XVCCSDT	183N⁷ + 286N⁶ + 220N⁵	15N⁵ + 343N⁴ + 189N³	15N⁵ + 131N⁴ + 99N³
EOM-XVCCSDTQ	729N⁸ + 916N⁷ + 598N⁶	34N⁶ + 965N⁵ + 474N⁴	34N⁶ + 311N⁵ + 170N⁴
EOM-XVCC5	2490N⁹ + 2794N⁸ + 1614N⁷	207N⁷ + 2649N⁶ + 1096N⁵	143N⁷ + 696N⁶ + 356N⁵
EOM-XVCC6	8295N¹⁰ + 8533N⁹ + 4432N⁸	348N⁸ + 8164N⁷ + 2782N⁶	254N⁸ + 2527N⁷ + 731N⁶
EOM-XVCC7	27 975N¹¹ + 26 506N¹⁰ + 12 457N⁹	569N⁹ + 24 962N⁸ + 8309N⁷	439N⁹ + 11 231N⁸ + 2570N⁷
EOM-XVCC8	96 603N¹² + 83 992N¹¹ + 35 934N¹⁰	934N¹⁰ + 78 522N⁹ + 25 129N⁸	762N¹⁰ + 51 289N⁹ + 11 274N⁸

^a

Factorization leads to a cost increase in XVCCS and thus was turned off in its implementation.

C. Factorization

The operation cost of XVCC and EOM-XVCC can be further reduced by factorization. Let us consider the sum-of-product tensor equation, E = A(BC) + A(BD), where the parentheses indicate the optimal order of contraction for each product. This sum can be evaluated as

I_{1} = C + D, I_{2} = B I_{1}, E = A I_{2} .

(80)

Two contractions by the common tensor B in I₂ = BC + BD are factored as I₂ = B(C + D), reducing the multiplication cost by half. Also, two contractions by A are also merged into one, further halving the multiplication cost. This procedure should not be confused with intermediate reuse.

Factorability of equations depends on the contraction order of each product, making strength reduction and factorization a coupled optimization problem. Since a complete solution of this coupled problem is exponentially complex, we approach it as two separate, sequential optimizations: strength reduction is performed first for each product, whereupon factorization is carried out by blanket search without altering the order of contraction. Hence, a sum of intermediates multiplied by the same tensor and the same permutation operator is factored. The resulting new intermediates are sums of tensors and they are recursively factored. In this work, the permutation operator itself is not decomposed into a product of smaller permutation operators, unlike in previous work for electronic CC.⁸⁴

An example of the procedure is given below and drawn in Fig. 14,

D_{6} = \hat{P} (ν μ) \sum_{i j} τ_{i ν} (τ_{j} W_{μ i j}),

(81)

D_{8} = \frac{1}{2!} \hat{P} (ν μ) \sum_{i j k} τ_{i ν} (τ_{j k} W_{μ i j k}),

(82)

D_{9} = \frac{1}{2!} \hat{P} (ν μ) \sum_{i j k} τ_{i ν} (τ_{j} (τ_{k} W_{μ i j k})),

(83)

D_{17} = \frac{1}{2!} \hat{P} (ν μ) \sum_{i j k} τ_{i ν} (τ_{j k μ} W_{i j k}),

(84)

D_{23} = \frac{1}{2!} \hat{P} (ν μ) \sum_{i j k l} τ_{i ν} (τ_{j k μ} (τ_{l} W_{i j k l})),

(85)

where the parentheses specify the optimal order of contraction for each product determined by the preceding strength reduction step. Factorization does not alter these orders. When evaluating the sum of these five terms, the common factor τ_iν is pulled out,

{(D_{sum})}_{ν μ} = \hat{P} (ν μ) \sum_{i} τ_{i ν} {(D_{6 + 8 + 9 + 17 + 23})}_{μ i},

(86)

where the intermediate D_6+8+9+17+23 is the sum of three intermediates,

{(D_{6 + 8 + 9 + 17 + 23})}_{μ i} = {(D_{8})}_{μ i} + {(D_{6 + 9})}_{μ i} + {(D_{17 + 23})}_{μ i}

(87)

with the subscripts indicating the diagrammatic origin of the intermediates. The first intermediate does not lend itself to any further factorization and is given by

{(D_{8})}_{μ i} = \frac{1}{2!} \sum_{j k} τ_{j k} W_{μ i j k} .

(88)

The second intermediate is further factored by the common tensor τ_j and is evaluated as

{(D_{6 + 9})}_{μ i} = \sum_{j} τ_{j} {(D_{6 + 9})}_{μ i j},

(89)

{(D_{6 + 9})}_{μ i j} = W_{μ i j} + \frac{1}{2!} \sum_{k} τ_{k} W_{μ i j k} .

(90)

The third intermediate is also further factored as

{(D_{17 + 23})}_{μ i} = \sum_{j k} τ_{j k μ} {(D_{17 + 23})}_{i j k},

(91)

{(D_{17 + 23})}_{i j k} = \frac{1}{2!} W_{i j k} + \frac{1}{2!} \sum_{l} τ_{l} W_{i j k l} .

(92)

FIG. 14.

View large Download slide

Factorable diagrams in the T₂ amplitude equation of XVCCSDT (see Appendix B).

The operation costs of XVCCm and EOM-XVCCm after factorization (in addition to strength reduction and intermediate reuse) are listed in Table IV. It should be understood that intermediate reuse (combined with factorization) applies to common intermediates appearing two or more times that are either a product or a sum-of-product of tensors. Since factorization is more likely for the final contraction step in our implementation, it is especially effective when the final step is the most costly. Factorization does not alter the overall scaling of the cost function, but generally reduces the prefactor multiplying the leading-order term, while it may increase or decrease prefactors of sub-leading-order terms.

Factorization is more important for EOM-XVCC and for higher values of m. In particular, for m > 4, factorization lowers the prefactors on the leading-order term of the cost functions. For m = 3 or 4, only the prefactors of the second and third leading-order terms of the cost functions are reduced. For m = 1 and 2, the cost is dominated by the contraction of the force-constant tensor with a tensor of T amplitudes, which occurs first in each product (see Appendix D); therefore, factorization is found least effective there.

D. Observed computational cost

In Fig. 15, the CPU time spent in one cycle of the iterative solution of the XVCCm or EOM-XVCCm equation (3 ≤ m ≤ 5) in a model QFF is plotted as a function of the number of modes N. Each plot falls accurately on the asymptotic cost-scaling line of O(N^m+2) for N ≥ 20, suggesting that the leading-order terms in the cost function dominate the cost. For instance, the XVCCSDTQ and EOM-XVCCSDTQ methods with the aforementioned optimizations are remarkably efficient, taking only on the order of ten seconds per iteration for N = 35, despite the fact that their algorithms do not exploit the index permutation symmetry. These costs are expected to be negligible as compared with the cost associated with generating all higher-order force constants at the electronic structure level appropriate for these accurate vibrational methods.

FIG. 15.

View large Download slide

CPU time (in seconds) spent in one cycle of the iterative solution of the XVCCm or EOM-XVCCm equation for a system with N modes in a QFF.

IV. NUMERICAL TESTS

The XVCCm and EOM-XVCCm methods were applied to the vibrational zero-point energies and frequencies of the water and formaldehyde molecules. Their equilibrium structure, normal modes, and QFF’s were calculated at the MP2/aug-cc-pVTZ level using nwchem ⁸⁶ and sindo.⁸⁷ The results were compared with those obtained from full VCI calculations with mavi ⁸⁸ using the 20 lowest-lying harmonic-oscillator wave functions of each mode as the basis set. The VCI results were converged to within 0.1 cm⁻¹ of the exact solutions for the QFF. Our results were also compared with the XVMP2(4) and XVMP2[4] results,⁵⁹ obtained either in the frequency-independent approximation to the self-energy or by solving the Dyson equation with the frequency-dependent self-energy.

A. Water

Table V compiles the errors in the XVCCm zero-point energy (E₀) and EOM-XVCCm frequencies (ν) from the VCI results for the water molecule. These errors decrease nearly monotonically with the excitation rank m, reaching the accuracy of a few cm⁻¹ at m = 3 (CCSDT) for both E₀ and ν. At m = 4 (CCSDTQ), the results are generally within tenths of cm⁻¹ of the VCI results with a few exceptions. The results for m ≥ 5 are given in Table VI and are converged with accuracy higher than necessary, e.g., within 0.001 cm⁻¹ at m = 8.

TABLE V.

The calculated zero-point energy (E₀) and frequencies (ν_i) of the fundamental transitions (in cm⁻¹) of the water molecule. Except for VCI, the errors from VCI are shown.

Method^a	E₀	ν₁	ν₂	ν₃
VCI	4642.4	3682.4	1556.6	3791.5
Harmonic	108.4	139.7	71.7	156.4
XVSCF(4)	104.5	280.2	−17.9	289.4
XVSCF[4]	21.5	105.3	11.7	116.1
XVCCS{4}	21.9	196.5	−1.2	206.7
XVCCS(4)	26.6	197.2	−2.4	207.0
XVCCS[4]	21.5	105.3	11.7	116.1
XVCCSD{4}	22.1	33.6	7.9	34.3
XVCCSD(4)	23.8	38.3	8.5	38.3
XVCCSD[4]	21.5	26.0	6.8	26.6
XVCCSDT{4}	0.2	−1.7	2.2	−2.0
XVCCSDT(4)	0.1	−2.4	2.2	−2.7
XVCCSDT[4]	0.4	−1.7	2.4	−1.9
XVCCSDTQ{4}	0.0	−0.7	0.2	−0.7
XVCCSDTQ(4)	−0.0	−1.0	0.1	−1.0
XVCCSDTQ[4]	0.0	−0.1	0.1	−0.1
XVMP2(4)^b	4.5	−24.6	3.8	−28.3
XVMP2[4]^b	−1.3	−18.9	2.4	−20.3
XVMP2(4)^c	4.5	−14.3	3.8	−20.8
XVMP2[4]^c	−1.3	−14.8	2.5	−17.0

Method^a	E₀	ν₁	ν₂	ν₃
VCI	4642.4	3682.4	1556.6	3791.5
Harmonic	108.4	139.7	71.7	156.4
XVSCF(4)	104.5	280.2	−17.9	289.4
XVSCF[4]	21.5	105.3	11.7	116.1
XVCCS{4}	21.9	196.5	−1.2	206.7
XVCCS(4)	26.6	197.2	−2.4	207.0
XVCCS[4]	21.5	105.3	11.7	116.1
XVCCSD{4}	22.1	33.6	7.9	34.3
XVCCSD(4)	23.8	38.3	8.5	38.3
XVCCSD[4]	21.5	26.0	6.8	26.6
XVCCSDT{4}	0.2	−1.7	2.2	−2.0
XVCCSDT(4)	0.1	−2.4	2.2	−2.7
XVCCSDT[4]	0.4	−1.7	2.4	−1.9
XVCCSDTQ{4}	0.0	−0.7	0.2	−0.7
XVCCSDTQ(4)	−0.0	−1.0	0.1	−1.0
XVCCSDTQ[4]	0.0	−0.1	0.1	−0.1
XVMP2(4)^b	4.5	−24.6	3.8	−28.3
XVMP2[4]^b	−1.3	−18.9	2.4	−20.3
XVMP2(4)^c	4.5	−14.3	3.8	−20.8
XVMP2[4]^c	−1.3	−14.8	2.5	−17.0

^a

XVCC for zero-point energy and EOM-XVCC for frequencies. Suffix “{4}” means the use of a harmonic reference, “(4)” means the use of a XVSCF(4) reference, and “[4]” means the use of a XVSCF[4] reference.

^b

In the frequency-independent approximation.⁵⁹

^c

Solutions of the Dyson equation with the frequency-dependent self-energy.⁵⁹

TABLE VI.

The calculated zero-point energy (E₀) and frequencies (ν_i) of the fundamental transitions (in cm⁻¹) of the water molecule. The errors from XVCC8 with the corresponding reference wave function are shown.

Method^a	E₀	ν₁	ν₂	ν₃
XVCC5{4}	−0.02	−0.05	0.02	−0.06
XVCC5(4)	−0.03	−0.12	0.01	−0.12
XVCC5[4]	−0.03	0.00	0.02	0.00
XVCC6{4}	−0.001	0.004	0.004	0.005
XVCC6(4)	−0.001	0.005	0.002	0.006
XVCC6[4]	−0.001	0.006	0.003	0.007
XVCC7{4}	−0.0004	0.003	0.001	0.003
XVCC7(4)	−0.0002	0.004	0.001	0.005
XVCC7[4]	−0.0005	−0.001	0.002	−0.001
XVCC8{4}	0	0	0	0
XVCC8(4)	0	0	0	0
XVCC8[4]	0	0	0	0

Method^a	E₀	ν₁	ν₂	ν₃
XVCC5{4}	−0.02	−0.05	0.02	−0.06
XVCC5(4)	−0.03	−0.12	0.01	−0.12
XVCC5[4]	−0.03	0.00	0.02	0.00
XVCC6{4}	−0.001	0.004	0.004	0.005
XVCC6(4)	−0.001	0.005	0.002	0.006
XVCC6[4]	−0.001	0.006	0.003	0.007
XVCC7{4}	−0.0004	0.003	0.001	0.003
XVCC7(4)	−0.0002	0.004	0.001	0.005
XVCC7[4]	−0.0005	−0.001	0.002	−0.001
XVCC8{4}	0	0	0	0
XVCC8(4)	0	0	0	0
XVCC8[4]	0	0	0	0

^a

XVCC for zero-point energy and EOM-XVCC for frequencies. Suffix “{4}” means the use of a harmonic reference, “(4)” means the use of a XVSCF(4) reference, and “[4]” means the use of a XVSCF[4] reference.

Note that the VCI results in this table were obtained with 20 basis functions per mode. The number of Hartree products involved in such a calculation grows exponentially with size or as O(b^N), where b is the number of basis functions and N is the number of modes. The numbers of the T and R amplitudes in XVCCm and EOM-XVCCm, on the other hand, increase only polynomially or as O(N^m). The values of b^N for VCI and N^m for XVCCSDTQ (N = 3, b = 20, and m = 4) are 8000 and 81, respectively, implying a much severer limitation in VCI’s applicability to larger molecules. The number of unknowns in Christiansen’s VCC (Ref. 36) grows as O(b^m). Therefore, it too suffers from the consequence of having a finite number of basis functions, unlike XVCCm, which is inherently an infinite-basis method.

The accuracy of E₀ from XVCCm is relatively insensitive to the choice of the reference wave function, despite much larger differences in the results of the reference methods themselves: harmonic approximation, XVSCF(4), and XVSCF[4]. This is in line with a similar observation⁷² made in the electronic CC results for correlation energies, which is understood to be due to the single-excitation operator which corrects any deficiency in the reference wave function. Likewise, the XVCCS(4) and XVCCS[4] results of E₀ are close to each other, while the XVSCF(4) results differ considerably from those of XVSCF[4]. This difference is due to the anharmonic shift in the geometry, neglected in XVSCF(4) but included in XVSCF[4]. Diagrammatically, this contribution is largely represented by E₁ of Fig. 3. The fact that this diagram is found among the XVCCS energy equation means that the effect of anharmonic geometry shift is included in XVCCS regardless of the reference wave function.

XVSCF[4] determines and adopts the first-order Dyson geometry⁵⁰ as the center of coordinates, where W_i = 0 for all i. The T₁ amplitude equation of XVCCS[4], Eq. (39), then has a trivial solution of τ_i = 0 for all i. This renders the XVCCS correction, $E_{0}^{(1)} - E_{0}^{(0)}$ in Eq. (28), zero, making the XVCCS[4] and XVSCF[4] values of E₀ identical to each other. This expected behavior is observed in Table V. Furthermore, substituting τ_i = 0 into the R₁ amplitude equation of EOM-XVCCS [Eq. (70)], we find the EOM-XVCCS[4] frequencies to be the eigenvalues of the $W_{i j}^{+}$ matrix, whose eigenvectors are known to define the first-order Dyson coordinates.^50,52 These eigenvalues are usually, but not always, close to the XVSCF[4] frequencies,⁵⁰ the latter being the diagonal elements (⁠ $W_{i i}^{+}$ ⁠) of the same matrix. This is also observed in Table V. Therefore, XVCCS can take into account the leading-order effects in E₀ and ν from anharmonic corrections both in geometry and coordinates.

It may be noted that XVMP2 performs considerably better than XVCCSD for E₀. This is in striking contrast with the electronic counterparts; CCSD has been observed to consistently outperform MP2.^1,2 This difference can be traced to the difference in the rank of the Hamiltonian in the two cases: the electronic Hamiltonian is a two-body operator, whereas the vibrational one is an n-body operator, where n is the rank of the Taylor-series PES (n = 4 in this calculation). The most important anharmonic correction to E₀ is diagram E₉ in Fig. 16 (in Appendix B), which involves cubic force constants and first appears at XVCCSDT. This contribution is to a great extent accounted for by the isomorphic diagram in XVMP2 labeled 2C′ in Ref. 59. This explains why only at m ≥ 3 does XVCCm outperform XVMP2. It is generally recommended to set the value of m equal to n, if not higher, because the overall cost of the calculation is dominated by the O(Nⁿ) step of the force-constant evaluation.

FIG. 16.

View large Download slide

Additional diagrams in the energy equation of XVCCSDT in a QFF.

The rather small change in E₀ from XVCCS to XVCCSD can also be understood as follows. The chief contribution of XVCCSD to E₀, diagram E₅ of Fig. 4, is nearly zero in the XVSCF reference, where $W_{i i}^{-} = 0$ holds.⁵² When the normal coordinates are close to the first-order Dyson coordinates, which is usually the case, $W_{i j}^{-} = W_{i j}^{+} \approx 0$ is also true for i≠j, making the remaining (i.e., off-diagonal) contributions to E₅ insignificant. Furthermore, the smallness of $W_{i j}^{-}$ implies the smallness of the T₂ amplitudes [cf. Eq. (53)]. Although XVCCSD uses cubic force constants in E₆, they are multiplied by these small T₂ amplitudes and can hardly alter E₀ of XVCCS. The frequencies from EOM-XVCCSD, on the other hand, are considerably better than EOM-XVCCS and are similar to those of XVMP2. This may be ascribed to diagram $S_{5}^{'}$ in Fig. 20 (in Appendix C) included in EOM-XVCCSD, which can describe the same essential cubic corrections accounted for by self-energy diagrams 2e′ and 2f′ (Ref. 59) of XVMP2.

FIG. 20.

View large Download slide

Additional diagrams in the R₁ amplitude equation of EOM-XVCCSD in a QFF.

For the zero-point energies and frequencies that are not affected by any resonance, we observe the following general order of accuracy: harmonic ≈ XVSCF(n) < XVCCS{n} ≈ XVCCS(n) < XVSCF[n] ≈ XVCCS[n] < XVCCSD < XVMP2 < XVCCSDT < XVCCSDTQ ≤ XVCC5 ≤ XVCC6 ≤ XVCC7 ≤ XVCC8, which differs from the order established in electronic-structure theory.^1,2 In a QFF (n = 4), XVCCSDTQ or m = 4 already achieves accuracy high enough for practical purposes and XVCCm with m > n seem to have diminishing returns. Here, XVCC stands for both XVCC and EOM-XVCC and the effects of the reference wave functions and the frequency-independent approximation in XVMP2 are insignificant because there are no resonances.

B. Formaldehyde

The XVCCm and EOM-XVCCm results for the formaldehyde molecule are summarized in Tables VII and VIII. The calculated zero-point energy and frequencies of the fundamentals (except for ν₁ and ν₅) confirm the general observations made above for the water molecule.

TABLE VII.

The calculated zero-point energy (E₀) and frequencies (ν_i) of the fundamental transitions (in cm⁻¹) of the formaldehyde molecule. Except for VCI, the errors from VCI are shown.

Method^a	E₀	ν₁	ν₂	ν₃	ν₄	ν₅	ν₃ν₆	ν₂ν₆	ν₆
VCI	5813.6	2837.2	1719.8	1501.5	1160.7	2873.3	2702.8	2989.2	1237.5
Harmonic	74.3	136.2	33.2	38.6	36.2	174.2	…	…	29.4
XVSCF(4)	73.0	185.0	41.0	10.8	−13.7	216.4	…	…	−3.3
XVSCF[4]	22.9	79.0	20.4	9.8	10.6	106.7	…	…	11.2
XVCCS{4}	24.0	134.8	31.9	9.4	−1.3	164.7	…	…	4.0
XVCCS(4)	25.5	134.1	31.9	9.7	−0.8	163.6	…	…	4.4
XVCCS[4]	22.9	79.1	20.9	9.3	10.6	106.7	…	…	11.1
XVCCSD{4}	23.5	30.1	6.7	5.1	6.7	38.1	31.6	34.6	6.4
XVCCSD(4)	24.1	31.6	6.6	5.4	7.5	38.9	32.2	35.4	6.8
XVCCSD[4]	22.9	24.2	6.0	4.2	6.0	31.0	35.5	33.0	5.6
XVCCSDT{4}	0.6	0.3	0.5	1.2	1.7	1.9	7.3	2.2	1.4
XVCCSDT(4)	0.6	−0.0	0.5	1.2	1.7	1.8	7.9	2.4	1.4
XVCCSDT[4]	0.8	0.1	0.5	1.2	1.9	1.3	6.1	1.4	1.4
XVCCSDTQ{4}	0.0	−0.2	0.0	0.1	0.1	0.5	2.4	0.9	0.1
XVCCSDTQ(4)	0.0	−0.3	0.0	0.1	0.1	0.4	2.4	0.9	0.1
XVCCSDTQ[4]	0.1	0.2	0.0	0.1	0.0	0.7	2.3	1.0	0.0
XVMP2(4)^b	1.3	−190.9	0.7	0.9	2.7	33.0	…	…	3.0
XVMP2[4]^b	−1.0	−15.9	0.1	1.3	1.6	−414.3	…	…	1.8
XVMP2(4)^c	1.3	−6.6	1.4	1.0	2.6	−0.2	6.6	29.3	3.0
XVMP2[4]^c	−1.0	−8.9	0.3	1.4	1.5	2.2	18.0	24.8	1.8

Method^a	E₀	ν₁	ν₂	ν₃	ν₄	ν₅	ν₃ν₆	ν₂ν₆	ν₆
VCI	5813.6	2837.2	1719.8	1501.5	1160.7	2873.3	2702.8	2989.2	1237.5
Harmonic	74.3	136.2	33.2	38.6	36.2	174.2	…	…	29.4
XVSCF(4)	73.0	185.0	41.0	10.8	−13.7	216.4	…	…	−3.3
XVSCF[4]	22.9	79.0	20.4	9.8	10.6	106.7	…	…	11.2
XVCCS{4}	24.0	134.8	31.9	9.4	−1.3	164.7	…	…	4.0
XVCCS(4)	25.5	134.1	31.9	9.7	−0.8	163.6	…	…	4.4
XVCCS[4]	22.9	79.1	20.9	9.3	10.6	106.7	…	…	11.1
XVCCSD{4}	23.5	30.1	6.7	5.1	6.7	38.1	31.6	34.6	6.4
XVCCSD(4)	24.1	31.6	6.6	5.4	7.5	38.9	32.2	35.4	6.8
XVCCSD[4]	22.9	24.2	6.0	4.2	6.0	31.0	35.5	33.0	5.6
XVCCSDT{4}	0.6	0.3	0.5	1.2	1.7	1.9	7.3	2.2	1.4
XVCCSDT(4)	0.6	−0.0	0.5	1.2	1.7	1.8	7.9	2.4	1.4
XVCCSDT[4]	0.8	0.1	0.5	1.2	1.9	1.3	6.1	1.4	1.4
XVCCSDTQ{4}	0.0	−0.2	0.0	0.1	0.1	0.5	2.4	0.9	0.1
XVCCSDTQ(4)	0.0	−0.3	0.0	0.1	0.1	0.4	2.4	0.9	0.1
XVCCSDTQ[4]	0.1	0.2	0.0	0.1	0.0	0.7	2.3	1.0	0.0
XVMP2(4)^b	1.3	−190.9	0.7	0.9	2.7	33.0	…	…	3.0
XVMP2[4]^b	−1.0	−15.9	0.1	1.3	1.6	−414.3	…	…	1.8
XVMP2(4)^c	1.3	−6.6	1.4	1.0	2.6	−0.2	6.6	29.3	3.0
XVMP2[4]^c	−1.0	−8.9	0.3	1.4	1.5	2.2	18.0	24.8	1.8

^a

XVCC for zero-point energy and EOM-XVCC for frequencies. Suffix “{4}” means the use of a harmonic reference, “(4)” means the use of a XVSCF(4) reference, and “[4]” means the use of a XVSCF[4] reference.

^b

In the frequency-independent approximation.⁵⁹

^c

Solutions of the Dyson equation with the frequency-dependent self-energy.⁵⁹

TABLE VIII.

The calculated zero-point energy (E₀) and frequencies (ν_i) of the fundamental transitions (in cm⁻¹) of the formaldehyde molecule. The errors from XVCC8 with the corresponding reference wave function are shown.

Method^a	E₀	ν₁	ν₂	ν₃	ν₄	ν₅	ν₃ν₆	ν₂ν₆	ν₆
XVCC5{4}	−0.02	−0.03	0.00	0.01	0.00	0.03	0.33	0.06	0.01
XVCC5(4)	−0.02	−0.06	0.00	0.01	0.00	0.01	0.32	0.05	0.01
XVCC5[4]	−0.02	0.02	0.00	0.01	−0.00	0.04	0.26	0.06	0.01
XVCC6{4}	−0.001	0.000	0.001	0.001	−0.001	0.005	0.067	0.011	0.001
XVCC6(4)	−0.001	−0.002	0.000	0.001	−0.001	0.003	0.062	0.010	0.001
XVCC6[4]	−0.000	0.002	0.000	0.000	−0.002	0.006	0.052	0.008	0.000
XVCC7{4}	−0.0003	0.0012	0.0002	0.0004	0.0003	0.0010	0.0131	0.0021	0.0005
XVCC7(4)	−0.0002	0.0014	0.0002	0.0004	0.0003	0.0010	0.0124	0.0022	0.0005
XVCC7[4]	−0.0004	−0.0011	0.0002	0.0006	0.0007	−0.0005	0.0100	0.0004	0.0007
XVCC8{4}	0	0	0	0	0	0	0	0	0
XVCC8(4)	0	0	0	0	0	0	0	0	0
XVCC8[4]	0	0	0	0	0	0	0	0	0

Method^a	E₀	ν₁	ν₂	ν₃	ν₄	ν₅	ν₃ν₆	ν₂ν₆	ν₆
XVCC5{4}	−0.02	−0.03	0.00	0.01	0.00	0.03	0.33	0.06	0.01
XVCC5(4)	−0.02	−0.06	0.00	0.01	0.00	0.01	0.32	0.05	0.01
XVCC5[4]	−0.02	0.02	0.00	0.01	−0.00	0.04	0.26	0.06	0.01
XVCC6{4}	−0.001	0.000	0.001	0.001	−0.001	0.005	0.067	0.011	0.001
XVCC6(4)	−0.001	−0.002	0.000	0.001	−0.001	0.003	0.062	0.010	0.001
XVCC6[4]	−0.000	0.002	0.000	0.000	−0.002	0.006	0.052	0.008	0.000
XVCC7{4}	−0.0003	0.0012	0.0002	0.0004	0.0003	0.0010	0.0131	0.0021	0.0005
XVCC7(4)	−0.0002	0.0014	0.0002	0.0004	0.0003	0.0010	0.0124	0.0022	0.0005
XVCC7[4]	−0.0004	−0.0011	0.0002	0.0006	0.0007	−0.0005	0.0100	0.0004	0.0007
XVCC8{4}	0	0	0	0	0	0	0	0	0
XVCC8(4)	0	0	0	0	0	0	0	0	0
XVCC8[4]	0	0	0	0	0	0	0	0	0

^a

XVCC for zero-point energy and EOM-XVCC for frequencies. Suffix “{4}” means the use of a harmonic reference, “(4)” means the use of a XVSCF(4) reference, and “[4]” means the use of a XVSCF[4] reference.

The ν₅ fundamental is known to undergo Fermi resonances with combinations ν₃ν₆ and ν₂ν₆. For this mode, XVMP2[4] in the frequency-independent approximation shows a sign of divergent perturbation series, with a large error of −414 cm⁻¹ (the VMP2 result is even worse).⁵⁹ XVMP2(4) in the same approximation, on the other hand, is reasonable (with an error of 33 cm⁻¹) but breaks down unexpectedly for nonresonant ν₁ with an error of −191 cm⁻¹.⁵⁹ Therefore, the behavior of XVMP2 in the frequency-independent approximation is rather unpredictable. This method also lacks roots for Fermi-resonant counterparts.

Solving the Dyson equation with the frequency-dependent XVMP2 self-energy rectifies the aforementioned problems of the frequency-independent approximation.⁵⁹ With this approach, the frequency of ν₅ is reproduced within 3 cm⁻¹ of VCI and the roots corresponding to Fermi-resonant counterparts, ν₃ν₆ and ν₂ν₆, can be located within 30 cm⁻¹. The calculated results of ν₁ also regain stability with an error of less than 10 cm⁻¹. The price one pays is that a numerically cumbersome pole search needs to be performed along each mode.

Without any pole search, EOM-XVCC yields results that display stable convergence and have generally high accuracy for these strongly correlated modes. EOM-XVCCSDT seems to systematically outperform XVMP2 with the frequency-dependent self-energy for all listed modes including the Fermi resonances. EOM-XVCCSDTQ yields near-exact results with the largest error being merely 2.4 cm⁻¹. Remarkably, EOM-XVCCSD, though lacking an ${\hat{R}}_{3}$ operator to fully describe cubic anharmonicity, locates all three Fermi-resonant transitions (ν₅, ν₃ν₆, and ν₂ν₆) with uniform, medium accuracy. XVSCF and EOM-XVCCS, as expected, do not have roots corresponding to ν₃ν₆ and ν₂ν₆ and give rather poor results for ν₁ and ν₅.

The combinations and overtones accessible by XVMP2 with the frequency-dependent self-energy are limited to those that are in resonance with a fundamental, which may be an advantage in some cases, but a disadvantage in others. EOM-XVCCm does not have such a limitation and has roots for all low-lying states involving up to m modes and is accurate for those involving up to (m − 1) modes. This is analogous to the electronic case: EOM-CCSD has roots for up to two-electron excited states and is accurate for one-electron excitation energies; for accurate results for two-electron excitation energies, one needs to use EOM-CCSDT.⁸⁹ Table IX attests to this general trend for two- and three-mode combinations or overtones. The errors in EOM-XVCCSD for two-mode excitations and those in EOM-XVCCSDT for three-mode excitations are on the same order of magnitude as those in EOM-XVCCS for nonresonant fundamentals (ν₂, ν₃, ν₄, and ν₆). EOM-XVCCSDTQ for three-mode excitations, EOM- XVCCSDT for two-mode excitations, and EOM-XVCCSD for nonresonant fundamentals have comparable, high accuracy.

TABLE IX.

The calculated frequencies of selected overtones and combinations (in cm⁻¹) of the formaldehyde molecule. Except for VCI, the errors from VCI are shown.

Method	(000200)	(000101)	(000002)	(001100)	(010100)	(002000)	(000300)	(000201)
VCI	2318.1	2399.6	2472.5	2660.6	2871.2	3000.8	3464.2	3555.2
EOM-XVCCSD[4]	35.7	28.3	31.0	77.7	38.3	23.2	…	…
EOM-XVCCSDT[4]	7.6	6.0	5.3	5.7	4.5	4.7	62.5	49.6
EOM-XVCCSDTQ[4]	2.1	2.0	2.1	2.0	1.8	2.0	14.2	10.9
EOM-XVCC5[4]	0.1	0.1	0.2	0.1	0.1	0.2	3.3	3.0
EOM-XVCC6[4]	0.0	0.0	0.0	0.0	0.0	0.0	0.5	0.3
EOM-XVCC7[4]	0.0	0.0	0.0	0.0	0.0	0.0	0.1	0.1
EOM-XVCC8[4]	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0

Method	(000200)	(000101)	(000002)	(001100)	(010100)	(002000)	(000300)	(000201)
VCI	2318.1	2399.6	2472.5	2660.6	2871.2	3000.8	3464.2	3555.2
EOM-XVCCSD[4]	35.7	28.3	31.0	77.7	38.3	23.2	…	…
EOM-XVCCSDT[4]	7.6	6.0	5.3	5.7	4.5	4.7	62.5	49.6
EOM-XVCCSDTQ[4]	2.1	2.0	2.1	2.0	1.8	2.0	14.2	10.9
EOM-XVCC5[4]	0.1	0.1	0.2	0.1	0.1	0.2	3.3	3.0
EOM-XVCC6[4]	0.0	0.0	0.0	0.0	0.0	0.0	0.5	0.3
EOM-XVCC7[4]	0.0	0.0	0.0	0.0	0.0	0.0	0.1	0.1
EOM-XVCC8[4]	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0

For transition frequencies to strongly correlated states involving at least two modes (Fermi resonances, combinations, and overtones), the following is the order of accuracy: XVCCSD < XVMP2 with pole search < XVCCSDT < XVCCSDTQ < XVCC5 < XVCC6 < XVCC7 ≤ XVCC8. The methods not mentioned in the list either do not have the corresponding roots or tend to display signs of nonconvergence.

V. CONCLUSION

We have presented an accurate, stable, efficient, and widely applicable vibrational CC (XVCC) theory for both zero-point energies and transition frequencies, the latter through the EOM formalism, of a molecule with an anharmonic PES. It relies on a Hartree product of harmonic-oscillator modals as a reference wave function and a truncated Taylor-series PES. They lead to streamlined sets of equations defining XVCC and EOM-XVCC, which mirror exactly those of electronic CC and EOM-CC theories. XVCC and EOM-XVCC are diagrammatically size consistent, free of modal basis sets, and thus an inherently infinite-basis theory.

We have established the rules for enumerating Feynman– Goldstone diagrams for XVCC and EOM-XVCC theories and the ones for translating them into algebraic equations. This method of derivation is not only expedient but also insightful because the diagrammatic topology can more immediately determine size consistency of a theory or the lack thereof as well as the relationships among various diagrammatic theories. The rules are justified by the underlying normal-ordered second-quantized method of derivation, which gives the same equations as the diagrammatic derivation in each case.

XVCC and EOM-XVCC theories are general in the sense that their formalisms are unchanged for any single-reference Hartree product of harmonic-oscillator modals, including the wave functions of XVSCF(n), XVSCF[n], SCP, or the harmonic approximation. Reflecting the centrality of XVSCF (which includes SCP) as the size-consistent single-Hartree-product method in many-body vibrational theory, some diagrammatic contributions in XVCC and EOM-XVCC vanish in this reference, simplifying their equations further. Low-order XVCC and EOM-XVCC results are also the same or very close to the XVSCF results for the same reason.

The XVCCm and EOM-XVCCm methods up to m = 8 have been formulated and implemented by computer algebra of the aforementioned diagram logic as well as of algorithmic optimizations such as strength reduction, intermediate reuse, and factorization. The numerical tests of these methods show that the results are smoothly and rapidly convergent towards the exact values. This is in contrast to XVMP2 in the frequency-independent approximation to the self-energy, which displays signs of divergent perturbation series not only for resonant modes but also for nonresonant modes. While this problem is rectified by solving the Dyson equation with the frequency-dependent self-energy in XVMP2, EOM-XVCC may be said to achieve the same with a more well-studied diagonalization algorithm instead of a pole search. XVCCm and EOM-XVCCm with m = n, where n is the truncation rank of the Taylor-series PES, seem to be ones with the best cost-accuracy balance and hence the greatest practical utility.

The VCC theory of Prasad and coworkers^27,32 can be identified as XVCC and EOM-XVCC with the XVSCF[n] reference. Our work is, therefore, built upon the previous studies of Prasad and coworkers,^27,32 but it generalizes the theories to a variety of reference wave functions, establishes the normal-ordered, diagrammatic, and automated derivations, and implements higher-order methods up to connected octuple excitations using computerized code optimization and synthesis. Christiansen’s VCC theory,³⁶ on the other hand, relies on VSCF and modals that are numerically defined and differs from XVCC; unlike XVCC, VCC likely retains many terms that are not size consistent and is expected to be fundamentally less efficient for large molecules and solids despite the exponential parameterization of its wave function. In VCC, the rank of disconnected excitations is also limited by the number of basis functions, which is not the case with basis-set-free XVCC.

VSCF is exact for any one-dimensional potential in the infinite basis-set limit, whereas XVSCF captures only a small portion of its anharmonic effect. For solids, VSCF reduces to XVSCF, the former giving the same results as the latter, but at a much greater computational cost. This suggests that Christiansen’s VCC may also perform well for small molecules with particularly strong anharmonicity, while reducing to more streamlined XVCC or the VCC theory of Prasad and coworkers in the thermodynamic limit. Our numerical results show that, for Morse-like anharmonicity, overtones and combinations, and anharmonic resonances, XVCC can nevertheless rectify the deficiency of XVSCF for small molecules. However, XVSCF and XVCC are expected to not work well for a double-well potential, for instance.

Furthermore, rectilinear coordinates such as normal coordinates, also underlying XVCC, may become awkward for large-amplitude vibrations. The methods that are suitable for such strong correlation and handle different coordinates flexibly include the multiconfiguration time-dependent Hartree (MCTDH) method,⁹⁰ which is variational, but may lack size consistency.⁹¹ XVMP2 and XVCC are designed primarily for weak correlation in the ground state and are nonvariational and diagrammatically size consistent.⁹¹ Excited states are usually strongly correlated but can be handled by either class of methods.⁵⁹ Therefore, it may be said that XVCC is not well suited for strong correlation in the ground state, but it is an excellent—arguably the best—size-consistent method for weak correlation in the ground state and for strong and weak correlation in excited states.

Acknowledgments

This work is supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Grant No. DE-FG02-11ER16211. J.A.F. is supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1144245. The authors thank Matthew R. Hermes for useful discussions and supplying the QFF’s and VCI results used in this work.

APPENDIX A: RIGOROUS DERIVATION OF THE CONNECTED FORM OF THE XVCC EQUATIONS

As before, we demand the XVCC wave function, $e^{\hat{T}} | Φ_{0} 〉$ ⁠, satisfy the vibrational Schrödinger equation,

{\hat{H}}_{N} e^{\hat{T}} | Φ_{0} 〉 = Δ E_{0}^{(m)} e^{\hat{T}} | Φ_{0} 〉,

(A1)

but only in the space of Hartree products reachable by $\hat{T}$ from |Φ₀〉. Prior to projecting this equation onto this space, we first act $e^{- \hat{T}}$ on both sides of the equation,

e^{- \hat{T}} {\hat{H}}_{N} e^{\hat{T}} | Φ_{0} 〉 = Δ E_{0}^{(m)} | Φ_{0} 〉 .

(A2)

The projection of this equation then yields the energy and amplitude equations, which read

〈 Φ_{0} | e^{- \hat{T}} {\hat{H}}_{N} e^{\hat{T}} | Φ_{0} 〉 = Δ E_{0}^{(m)},

(A3)

〈 Φ_{ν_{1} ν_{2} \dots ν_{n}} | e^{- \hat{T}} {\hat{H}}_{N} e^{\hat{T}} | Φ_{0} 〉 = 0,

(A4)

for all ν₁, …, ν_n and each n (1 ≤ n ≤ m). It may be understood that the projection space 〈Φ₀|∪{〈Φ_{ν₁ν₂⋯ν_n}|} used in the derivation in the main text spans the same space $〈 Φ_{0} | e^{- \hat{T}} \cup {〈 Φ_{ν_{1} ν_{2} \dots ν_{n}} | e^{- \hat{T}}}$ as above and, therefore, the two derivations lead to one and the same theory.

The similarity-transformed Hamiltonian, $e^{- \hat{T}} {\hat{H}}_{N} e^{\hat{T}}$ ⁠, can be rewritten using the Baker–Campbell–Hausdorff expansion^1,2 as

e^{- \hat{T}} {\hat{H}}_{N} e^{\hat{T}} = {\hat{H}}_{N} + [{\hat{H}}_{N}, \hat{T}] + \frac{1}{2!} [[{\hat{H}}_{N}, \hat{T}], \hat{T}] + \frac{1}{3!} [[[{\hat{H}}_{N}, \hat{T}], \hat{T}], \hat{T}] + \dots,

(A5)

which terminates with the n-fold commutator, where n is the truncation rank of the Taylor-series PES. The commutators immediately guarantee the diagrammatically connected structure of each surviving term. Therefore, we may write

e^{- \hat{T}} {\hat{H}}_{N} e^{\hat{T}} = {({\hat{H}}_{N} e^{\hat{T}})}_{C},

(A6)

underscoring the equivalence of Eqs. (A3) and (A4) with Eqs. (29) and (40), respectively.

APPENDIX B: XVCCSDT

Here, we document the algebraic and diagrammatic forms of the XVCCSDT energy and amplitude equations in a QFF. The diagrams in the energy equation that are not already listed in the XVCCSD energy equation are shown in Fig. 16, and the corresponding algebraic expressions are given in Table X. As the numerical results illustrate (see above), the inclusion of connected triple excitations is critical for highly accurate zero-point energies. This is attributed to the importance of cubic force constants and diagram E₉. The T₁ and T₂ amplitude equations are given by

0 = T_{1}^{(2)} + \sum_{i = 13}^{18} S_{i} \equiv T_{1}^{(3)},

(B1)

0 = T_{2}^{(2)} + \sum_{i = 13}^{24} D_{i} \equiv T_{2}^{(3)},

(B2)

where $T_{1}^{(2)}$ and $T_{2}^{(2)}$ are defined by Eqs. (51) and (53), respectively. The additional diagrams in the T₁ and T₂ amplitude equations are shown in Figs. 17 and 18. The algebraic interpretation of the diagrams are given in Tables XI and XII. The T₃ amplitude equation is written as

0 = \sum_{i = 1}^{21} T_{i} \equiv T_{3}^{(3)} .

(B3)

The diagrammatic form of this equation is shown in Fig. 19 with the corresponding algebraic interpretation in Table XIII.

TABLE X.

The energy equation of XVCCSDT in a QFF. Einstein summation convention is adopted.

Δ E_{0}^{(3)} = Δ E_{0}^{(2)} + \frac{1}{3!} W_{i j k} τ_{i j k} + \frac{1}{3!} W_{i j k l} τ_{i} τ_{j k l}

FIG. 17.

View large Download slide

Additional diagrams in the T₁ amplitude equation of XVCCSDT in a QFF.

FIG. 18.

View large Download slide

Additional diagrams in the T₂ amplitude equation of XVCCSDT in a QFF.

TABLE XI.

The T₁ amplitude equation of XVCCSDT in a QFF. Einstein summation convention is adopted.

T_{1}^{(3)} \equiv T_{1}^{(2)} + \frac{1}{2!} W_{i j}^{-} τ_{i j ν} + \frac{1}{3!} W_{ν i j k} τ_{i j k} + \frac{1}{2!} W_{i j k} τ_{i} τ_{j k ν}

+ \frac{1}{2!^{2}} W_{i j k l} τ_{i j} τ_{k l ν} + \frac{1}{3!} W_{i j k l} τ_{i ν} τ_{j k l} + \frac{1}{2!^{2}} W_{i j k l} τ_{i} τ_{j} τ_{k l ν}

= 0

TABLE XII.

The T₂ amplitude equation of XVCCSDT in a QFF. Einstein summation convention is adopted.

T_{2}^{(3)}

≡

T_{2}^{(2)} + \hat{P} (ν μ) (\frac{1}{2!} W_{i} τ_{i ν μ} + \frac{1}{2!} W_{ν i j} τ_{i j μ} + \frac{1}{2!} W_{i j}^{-} τ_{i} τ_{j ν μ}

+ \frac{1}{2!} W_{ν i j k} τ_{i} τ_{j k μ} + \frac{1}{2!} W_{i j k} τ_{i ν} τ_{j k μ} + \frac{1}{2!^{2}} W_{i j k} τ_{i j} τ_{k ν μ}

+ \frac{1}{2!^{2}} W_{i j k} τ_{i} τ_{j} τ_{k ν μ} + \frac{1}{2! 3!} W_{i j k l} τ_{i j k} τ_{l ν μ} + \frac{1}{2!^{3}} W_{i j k l} τ_{i j ν} τ_{k l μ}

+ \frac{1}{2!^{2}} W_{i j k l} τ_{i} τ_{j k} τ_{l ν μ} + \frac{1}{2!} W_{i j k l} τ_{i} τ_{j ν} τ_{k l μ} + \frac{1}{2! 3!} W_{i j k l} τ_{i} τ_{j} τ_{k} τ_{l ν μ})

=

0

FIG. 19.

View large Download slide

Diagrams in the T₃ amplitude equation of XVCCSDT in a QFF.

TABLE XIII.

The T₃ amplitude equation of XVCCSDT in a QFF. Einstein summation convention is adopted.

T_{3}^{(3)}

≡

\hat{P} (ν μ λ) (\frac{1}{3!} W_{ν μ λ} + \frac{1}{3!} W_{ν μ λ i} τ_{i} + \frac{1}{2!} W_{ν μ i} τ_{i λ} + \frac{1}{2!} W_{ν i}^{+} τ_{i μ λ}

+ \frac{1}{2!^{2}} W_{ν μ i j} τ_{i j λ} + \frac{1}{2!} W_{ν μ i j} τ_{i} τ_{j λ} + \frac{1}{2!} W_{ν i j} τ_{i μ} τ_{j λ} + \frac{1}{2!} W_{ν i j} τ_{i} τ_{j μ λ}

+ \frac{1}{2!} W_{i j}^{-} τ_{i ν} τ_{j μ λ} + \frac{1}{2!} W_{ν i j k} τ_{i μ} τ_{j k λ} + \frac{1}{2!^{2}} W_{ν i j k} τ_{i j} τ_{k μ λ}

+ \frac{1}{2!^{2}} W_{i j k} τ_{i j ν} τ_{k μ λ} + \frac{1}{2!} W_{ν i j k} τ_{i} τ_{j μ} τ_{k λ} + \frac{1}{3!} W_{i j k} τ_{i ν} τ_{j μ} τ_{k λ}

+ \frac{1}{2!^{2}} W_{ν i j k} τ_{i} τ_{j} τ_{k μ λ} + \frac{1}{2!} W_{i j k} τ_{i} τ_{j ν} τ_{k μ λ} + \frac{1}{2!^{2}} W_{i j k l} τ_{i} τ_{j k ν} τ_{l μ λ}

+ \frac{1}{2!^{2}} W_{i j k l} τ_{i ν} τ_{j μ} τ_{k l λ} + \frac{1}{2!^{2}} W_{i j k l} τ_{i j} τ_{k ν} τ_{l μ λ} + \frac{1}{3!} W_{i j k l} τ_{i} τ_{j ν} τ_{k μ} τ_{l λ}

+ \frac{1}{2!^{2}} W_{i j k l} τ_{i} τ_{j} τ_{k ν} τ_{l μ λ})

=

0

APPENDIX C: EOM-XVCCSD

The R₁ and R₂ amplitude equations of EOM-XVCCSD in a QFF are written as

R_{1}^{(2)} \equiv R_{1}^{(1)} + \sum_{i = 4}^{18} S_{i}^{'} = ω_{k}^{(2)} r_{ν},

(C1)

R_{2}^{(2)} \equiv \sum_{i = 1}^{17} D_{i}^{'} = ω_{k}^{(2)} r_{μ ν},

(C2)

where $R_{1}^{(1)}$ is defined by Eq. (70). The additional diagrams in the R₁ amplitude equation are depicted in Fig. 20. Their corresponding algebraic interpretation is given in Table XIV. The diagrammatic form of the R₂ amplitude equation is drawn in Fig. 21, while the algebraic form is given in Table XV. EOM-XVCCSD provides the first sizable correction to the vibrational frequencies. This is attributed to the inclusion of diagrams such as $S_{5}^{'}$ ⁠, which account for the important cubic force constants that are not multiplied by near-zero T₂ amplitudes. In particular, using a XVSCF[4] reference, EOM-XVCCS yields frequencies that are equivalent to the eigenvalues of the $W_{i j}^{+}$ matrix, completely neglecting cubic force constants. This is rectified by EOM-XVCCSD giving results with accuracy comparable to XVMP2.

TABLE XIV.

The R₁ amplitude equation of EOM-XVCCSD in a QFF. Einstein summation convention is adopted.

R_{1}^{(2)}

≡

R_{1}^{(1)} + W_{i} r_{i ν} + \frac{1}{2!} W_{ν i j} r_{i j} + W_{i j}^{-} τ_{j ν} r_{i} + W_{i j}^{-} τ_{j} r_{i ν}

+ \frac{1}{2!} W_{ν i j k} τ_{j k} r_{i} + \frac{1}{2!} W_{ν i j k} τ_{k} r_{i j} + \frac{1}{2!} W_{i j k} τ_{j k} r_{i ν} + \frac{1}{2!} W_{i j k} τ_{k ν} r_{i j}

+ W_{i j k} τ_{j} τ_{k ν} r_{i} + \frac{1}{2!} W_{i j k} τ_{j} τ_{k} r_{i ν} + \frac{1}{2!} W_{i j k l} τ_{j k} τ_{l ν} r_{i}

+ \frac{1}{2!} W_{i j k l} τ_{k} τ_{l ν} r_{i j} + \frac{1}{2!} W_{i j k l} τ_{j} τ_{k l} r_{i ν} + \frac{1}{2!} W_{i j k l} τ_{j} τ_{k} τ_{l ν} r_{i}

+ \frac{1}{3!} W_{i j k l} τ_{j} τ_{k} τ_{l} r_{i ν}

=

ω_{k}^{(2)} r_{ν}

FIG. 21.

View large Download slide

Diagrams in the R₂ amplitude equation of EOM-XVCCSD in a QFF.

TABLE XV.

The R₂ amplitude equation of EOM-XVCCSD in a QFF. Einstein summation convention is adopted.

R_{2}^{(2)}

≡

\hat{P} (ν μ) (\frac{1}{2!} W_{ν μ i} r_{i} + W_{ν i}^{+} r_{i μ} + \frac{1}{2!^{2}} W_{ν μ i j} r_{i j} + \frac{1}{2!} W_{ν μ i j} τ_{j} r_{i}

+ W_{ν i j} τ_{j μ} r_{i} + W_{ν i j} τ_{j} r_{i μ} + W_{i j}^{-} τ_{j μ} r_{i ν} + \frac{1}{2!} W_{ν i j k} τ_{j k} r_{i μ}

+ \frac{1}{2!} W_{ν i j k} τ_{k μ} r_{i j} + W_{ν i j k} τ_{j} τ_{k μ} r_{i} + \frac{1}{2!} W_{i j k} τ_{j ν} τ_{k μ} r_{i}

+ \frac{1}{2!} W_{ν i j k} τ_{j} τ_{k} r_{i μ} + W_{i j k} τ_{j} τ_{k μ} r_{i ν} + \frac{1}{2!} W_{i j k l} τ_{j k} τ_{l μ} r_{i ν}

+ \frac{1}{2!^{2}} W_{i j k l} τ_{k ν} τ_{l μ} r_{i j} + \frac{1}{2!} W_{i j k l} τ_{j} τ_{k ν} τ_{l μ} r_{i} + \frac{1}{2!} W_{i j k l} τ_{j} τ_{k} τ_{l μ} r_{i ν})

=

ω_{k}^{(2)} r_{ν μ}

APPENDIX D: COST SCALING OF XVCCm AND EOM-XVCCm IN AN nTH-ORDER PES

Theorem.

The peak operation cost of XVCCm or EOM-XVCCm in an nth-order PES can be expressed as O(N^C), where N is the number of vibrational degrees of freedom and C is an integer. The lowest value of C (C_min) is m + ⌊n/2⌋ for m ≥ n and is not greater than n + ⌊(m − 1)/2⌋ for m < n.

Remark.

The proof is given below separately for m ≥ n and for m < n. In both cases, we first show that the aforementioned values of C are an upper bound; for any tensor contraction appearing in XVCCm or EOM-XVCCm, there is at least one contraction order whose cost function rank does not exceed the respective value. For m ≥ n, we additionally show that the aforementioned value of C is also a lower bound; there is always one tensor contraction whose minimal cost function rank is exactly this value. Together they prove that the expression of C given above is minimal for m ≥ n.

In what follows, we call an R amplitude of EOM-XVCC a T amplitude because they need not be distinguished in this context. We also only need to consider sequential contraction orders in which the force-constant tensor is always contracted first without losing rigor of the proof. Let T^(k) be the kth T amplitude to be contracted, having $T_{ext}^{(k)}$ external lines and $T_{int}^{(k)}$ internal lines. Let I^(k) be the kth intermediate. The symbol I^(k) also denotes its size, i.e., the total number of lines. I⁽⁰⁾ is the force-constant tensor, i.e., I⁽⁰⁾ = n.

In the contraction of I^(k−1) with T^(k), the number of summation indexes is $T_{int}^{(k)}$ ⁠, because two T amplitudes never share common internal lines and, therefore, all internal lines of T^(k) must connect with I^(k−1). The size of the product, i.e., the kth intermediate, is $I^{(k - 1)} + T_{ext}^{(k)} - T_{int}^{(k)}$ and the cost function rank of the contraction is $I^{(k - 1)} + T_{ext}^{(k)}$ ⁠.

The theorem applies not only to XVCCm and EOM-XVCCm but also to other CC and EOM-CC methods for electrons and nucleons.

Lemma.

Let h be the total number of T amplitudes in a tensor contraction. Given intermediate I⁽ⁱ⁾ (0 ≤ i ≤ h), there is at least one contraction order of remaining (h − i) T amplitudes in which none of the subsequent intermediates has size exceeding max(I⁽ⁱ⁾, I^(h)).

Proof of the lemma.

The size of the kth intermediate is given by

I^{(k)} = I^{(i)} + \sum_{l = i + 1}^{k} (T_{ext}^{(l)} - T_{int}^{(l)}) .

(D1)

Let T^(k) amplitudes be ordered as follows: T^(k)’s with negative values of the difference in the external and internal lines (⁠

T_{ext}^{(k)} - T_{int}^{(k)}

⁠) are contracted first, next by T^(k)’s with zero difference, and finally by T^(k)’s with positive values of the difference. Then, I^(k) as a function of k for this contraction order is a convex function and thus cannot exceed the greater of their values at the end points, I^(k) ≤ max(I⁽ⁱ⁾, I^(h)).

Proof of the theorem for m ≥ n.

First, we show that C = m + ⌊n/2⌋ is an upper bound. We consider the following two cases: (1) there is at least one T^(k) with $T_{int}^{(k)} \geq ⌈ n / 2 ⌉$ ⁠; (2) otherwise.

In case 1, the T amplitude with $T_{int}^{(1)} \geq ⌈ n / 2 ⌉$ is contracted first with the force-constant tensor. The cost function rank of this contraction is given by

C = n + T_{ext}^{(1)} \leq n + m - T_{int}^{(1)} \leq n + m - ⌈ n / 2 ⌉ = m + ⌊ n / 2 ⌋,

(D2)

where we use $T_{ext}^{(1)} + T_{int}^{(1)} \leq m$ in the first inequality. The size of the resulting intermediate is

I^{(1)} = n + T_{ext}^{(1)} - T_{int}^{(1)} \leq m + ⌊ n / 2 ⌋ - ⌈ n / 2 ⌉ \leq m,

(D3)

where Eq. (D2) is used in the first inequality.

Next, we consider case 2 and the second and subsequent contractions in case 1 together. In the former, $T_{int}^{(k)} \leq ⌊ n / 2 ⌋$ for all k. In the latter, $T_{int}^{(k)} \leq ⌊ n / 2 ⌋$ for all k ≥ 2. Here, we contract T amplitudes in one of the orders of the lemma. As per the lemma, this guarantees I^(k) ≤ m in both cases [see Eq. (D3) in case 1]. The sizes of the kth and (k − 1)th intermediates are related to each other by

I^{(k)} = I^{(k - 1)} + T_{ext}^{(k)} - T_{int}^{(k)} \leq m .

(D4)

Therefore, the cost function rank of the kth contraction is

C = I^{(k - 1)} + T_{ext}^{(k)} \leq m + T_{int}^{(k)} \leq m + ⌊ n / 2 ⌋ .

(D5)

Equations (D2) and (D5) prove C_min ≤ m + ⌊n/2⌋.

Second, we show that the same value of C is a lower bound. There is always a binary product of the force-constant tensor, I⁽⁰⁾, and a T_m amplitude with ⌈n/2⌉ internal lines. Note that ⌈n/2⌉ cannot exceed n or m because m ≥ n ≥ 1. Furthermore, the number of external lines is m + n − 2⌈n/2⌉ ≤ m. Hence, this is a valid tensor contraction. The cost function rank of this contraction is

C = I^{(0)} + T_{ext}^{(0)} = n + m - ⌈ n / 2 ⌉ = m + ⌊ n / 2 ⌋,

(D6)

which proves C_min ≥ m + ⌊n/2⌋.

Together, we obtain C_min = m + ⌊n/2⌋ for m ≥ n.

Remark.

This is an example of the optimization problem in computer science known as the matrix chain multiplication.^84,92,93 In this case (m ≥ n), we have shown that there is an analytical solution that minimizes the exponent of the leading-order cost function and its value.

Proof of the theorem for m < n.

We contract T amplitudes in one of the orders of the lemma. As per the lemma, this ensures I^(k) ≤ n. We consider the following three cases in each contraction: (1) $T_{ext}^{(k)} < T_{int}^{(k)}$ ⁠, (2) $T_{ext}^{(k)} > T_{int}^{(k)}$ ⁠, and (3) $T_{ext}^{(k)} = T_{int}^{(k)}$ ⁠.

In case 1, since $T_{ext}^{(k)} + T_{int}^{(k)} \leq m$ ⁠, we have $T_{ext}^{(k)} \leq ⌊ (m - 1) / 2 ⌋$ ⁠. The cost function rank of the kth contraction is

C = I^{(k - 1)} + T_{ext}^{(k)} \leq n + ⌊ (m - 1) / 2 ⌋ .

(D7)

In case 2, similarly, $T_{int}^{(k)} \leq ⌊ (m - 1) / 2 ⌋$ holds. The size of the kth intermediate is given by

I^{(k)} = I^{(k - 1)} + T_{ext}^{(k)} - T_{int}^{(k)} \leq n,

(D8)

which implies

C = I^{(k - 1)} + T_{ext}^{(k)} \leq n + T_{int}^{(k)} \leq n + ⌊ (m - 1) / 2 ⌋ .

(D9)

In case 3, m must be even and $T_{ext}^{(k)} = T_{int}^{(k)} = m / 2$ ⁠. The intermediate resulting from the contraction of this T amplitude has the size,

I^{(k)} = I^{(k - 1)} + T_{ext}^{(k)} - T_{int}^{(k)} \leq m < n .

(D10)

The first inequality follows from the fact that I^(l) is a convex function with its minimum occurring at l = k, where $T_{ext}^{(k)} - T_{int}^{(k)} = 0$ ⁠. The cost function rank of this step is

C = I^{(k - 1)} + T_{ext}^{(k)} \leq n + T_{int}^{(k)} - 1 \leq n + m / 2 - 1 .

(D11)

Equations (D7), (D9), and (D11) prove C_min ≤ n + ⌊(m − 1)/2⌋ for m < n.

Remark.

This upper-bound expression for m < n happens to agree with the optimal cost function ranks for all cases of m < n in Table IV. However, we have identified many other combinations of m < n whose optimal ranks are lower.

REFERENCES

1.

R. J.

Bartlett

and

M.

Musiał

,

Rev. Mod. Phys.

79

,

291

(

2007

).

https://doi.org/10.1103/RevModPhys.79.291

Google Scholar

Crossref

2.

I.

Shavitt

and

R. J.

Bartlett

,

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

(

Cambridge University Press

,

Cambridge

,

2009

).

Google Scholar

Crossref

3.

F.

Coester

,

Nucl. Phys.

7

,

421

(

1958

).

https://doi.org/10.1016/0029-5582(58)90280-3

Google Scholar

Crossref

4.

F.

Coester

and

H.

Kümmel

,

Nucl. Phys.

9

,

225

(

1958

).

https://doi.org/10.1016/0029-5582(58)90394-8

Google Scholar

Crossref

5.

J.

Čížek

,

J. Chem. Phys.

45

,

4256

(

1966

).

https://doi.org/10.1063/1.1727484

Google Scholar

Crossref

6.

J.

Paldus

,

J.

Čížek

, and

I.

Shavitt

,

Phys. Rev. A

5

,

50

(

1972

).

https://doi.org/10.1103/PhysRevA.5.50

Google Scholar

Crossref

7.

J. A.

Pople

,

R.

Krishnan

,

H. B.

Schlegel

, and

J. S.

Binkley

,

Int. J. Quantum Chem.

14

,

545

(

1978

).

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

Google Scholar

Crossref

8.

R. J.

Bartlett

and

G. D.

Purvis

,

Int. J. Quantum Chem.

14

,

561

(

1978

).

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

Google Scholar

Crossref

9.

G. D.

Purvis

III and

R. J.

Bartlett

,

J. Chem. Phys.

76

,

1910

(

1982

).

https://doi.org/10.1063/1.443164

Google Scholar

Crossref

10.

S.

Hirata

,

Theor. Chem. Acc.

129

,

727

(

2011

).

https://doi.org/10.1007/s00214-011-0954-4

Google Scholar

Crossref

11.

S.

Hirata

,

M.

Keçeli

,

Y.

Ohnishi

,

O.

Sode

, and

K.

Yagi

,

Annu. Rev. Phys. Chem.

63

,

131

(

2012

).

https://doi.org/10.1146/annurev-physchem-032511-143718

Google Scholar

Crossref

PubMed

12.

R. F.

Bishop

,

Theor. Chim. Acta

80

,

95

(

1991

).

https://doi.org/10.1007/BF01119617

Google Scholar

Crossref

13.

W.

Förner

,

Int. J. Quantum Chem.

43

,

221

(

1992

).

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

Google Scholar

Crossref

14.

M.

Yu

,

S.

Kalvoda

, and

M.

Dolg

,

Chem. Phys.

224

,

121

(

1997

).

https://doi.org/10.1016/S0301-0104(97)00257-7

Google Scholar

Crossref

15.

P.

Reinhardt

,

Theor. Chem. Acc.

104

,

426

(

2000

).

https://doi.org/10.1007/s002140000169

Google Scholar

Crossref

16.

S.

Hirata

,

I.

Grabowski

,

M.

Tobita

, and

R. J.

Bartlett

,

Chem. Phys. Lett.

345

,

475

(

2001

).

https://doi.org/10.1016/S0009-2614(01)00897-1

Google Scholar

Crossref

17.

S.

Hirata

,

R.

Podeszwa

,

M.

Tobita

, and

R. J.

Bartlett

,

J. Chem. Phys.

120

,

2581

(

2004

).

https://doi.org/10.1063/1.1637577

Google Scholar

Crossref

PubMed

18.

J. J.

Shepherd

,

T. M.

Henderson

, and

G. E.

Scuseria

,

J. Chem. Phys.

140

,

124102

(

2014

).

https://doi.org/10.1063/1.4867783

Google Scholar

Crossref

PubMed

19.

D. J.

Dean

,

J. R.

Gour

,

G.

Hagen

,

M.

Hjorth-Jensen

,

K.

Kowalski

,

T.

Papenbrock

,

P.

Piecuch

, and

M.

Wloch

,

Nucl. Phys. A

752

,

299C

(

2005

).

https://doi.org/10.1016/j.nuclphysa.2005.02.041

Google Scholar

Crossref

20.

P.

Piecuch

,

M.

Wloch

,

J. R.

Gour

,

D. J.

Dean

,

M.

Hjorth-Jensen

, and

T.

Papenbrock

,

AIP Conf. Proc.

777

,

28

(

2005

).

https://doi.org/10.1063/1.1996870

Google Scholar

Crossref

21.

O.

Christiansen

,

Phys. Chem. Chem. Phys.

9

,

2942

(

2007

).

https://doi.org/10.1039/b618764a

Google Scholar

Crossref

PubMed

22.

M. R.

Hermes

and

S.

Hirata

,

Int. Rev. Phys. Chem.

34

,

71

(

2015

).

https://doi.org/10.1080/0144235X.2014.1001220

Google Scholar

Crossref

23.

U. B.

Kaulfuss

and

M.

Altenbokum

,

Phys. Rev. D

33

,

3658

(

1986

).

https://doi.org/10.1103/PhysRevD.33.3658

Google Scholar

Crossref

24.

R. F.

Bishop

and

M. F.

Flynn

,

Phys. Rev. A

38

,

2211

(

1988

).

https://doi.org/10.1103/PhysRevA.38.2211

Google Scholar

Crossref

25.

R. F.

Bishop

,

M. C.

Bosca

, and

M. F.

Flynn

,

Phys. Rev. A

40

,

3484

(

1989

).

https://doi.org/10.1103/PhysRevA.40.3484

Google Scholar

Crossref

26.

M. D.

Prasad

,

J. Chem. Phys.

88

,

7005

(

1988

).

https://doi.org/10.1063/1.454399

Google Scholar

Crossref

27.

V.

Nagalakshmi

,

V.

Lakshminarayana

,

G.

Sumithra

, and

M. D.

Prasad

,

Chem. Phys. Lett.

217

,

279

(

1994

).

https://doi.org/10.1016/0009-2614(93)E1380-Y

Google Scholar

Crossref

28.

G. M.

Sastry

and

M. D.

Prasad

,

Chem. Phys. Lett.

228

,

213

(

1994

).

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

Google Scholar

Crossref

29.

G.

Latha

and

M. D.

Prasad

,

Chem. Phys. Lett.

241

,

215

(

1995

).

https://doi.org/10.1016/0009-2614(95)00631-D

Google Scholar

Crossref

30.

G. S.

Latha

and

M. D.

Prasad

,

J. Chem. Phys.

105

,

2972

(

1996

).

https://doi.org/10.1063/1.472170

Google Scholar

Crossref

31.

M. D.

Prasad

,

Int. J. Mol. Sci.

3

,

447

(

2002

).

https://doi.org/10.3390/i3050447

Google Scholar

Crossref

32.

S.

Banik

,

S.

Pal

, and

M. D.

Prasad

,

J. Chem. Phys.

129

,

134111

(

2008

).

https://doi.org/10.1063/1.2982502

Google Scholar

Crossref

PubMed

33.

S.

Banik

,

S.

Pal

, and

M. D.

Prasad

,

J. Chem. Theory Comput.

6

,

3198

(

2010

).

https://doi.org/10.1021/ct1003669

Google Scholar

Crossref

PubMed

34.

S.

Banik

,

S.

Pal

, and

M. D.

Prasad

,

J. Chem. Phys.

137

,

114108

(

2012

).

https://doi.org/10.1063/1.4753422

Google Scholar

Crossref

PubMed

35.

M. R.

Hermes

,

M.

Keçeli

, and

S.

Hirata

,

J. Chem. Phys.

136

,

234109

(

2012

).

https://doi.org/10.1063/1.4729602

Google Scholar

Crossref

PubMed

36.

O.

Christiansen

,

J. Chem. Phys.

120

,

2149

(

2004

).

https://doi.org/10.1063/1.1637579

Google Scholar

Crossref

PubMed

37.

O.

Christiansen

,

J. Chem. Phys.

122

,

194105

(

2005

).

https://doi.org/10.1063/1.1899156

Google Scholar

Crossref

PubMed

38.

P.

Seidler

and

O.

Christiansen

,

J. Chem. Phys.

126

,

204101

(

2007

).

https://doi.org/10.1063/1.2734970

Google Scholar

Crossref

PubMed

39.

P.

Seidler

,

M. B.

Hansen

, and

O.

Christiansen

,

J. Chem. Phys.

128

,

154113

(

2008

).

https://doi.org/10.1063/1.2907860

Google Scholar

Crossref

PubMed

40.

P.

Seidler

,

E.

Matito

, and

O.

Christiansen

,

J. Chem. Phys.

131

,

034115

(

2009

).

https://doi.org/10.1063/1.3158946

Google Scholar

Crossref

PubMed

41.

P.

Seidler

and

O.

Christiansen

,

J. Chem. Phys.

131

,

234109

(

2009

).

https://doi.org/10.1063/1.3272796

Google Scholar

Crossref

PubMed

42.

P.

Seidler

,

M.

Sparta

, and

O.

Christiansen

,

J. Chem. Phys.

134

,

054119

(

2011

).

https://doi.org/10.1063/1.3536499

Google Scholar

Crossref

PubMed

43.

K. M.

Christoffel

and

J. M.

Bowman

,

Chem. Phys. Lett.

85

,

220

(

1982

).

https://doi.org/10.1016/0009-2614(82)80335-7

Google Scholar

Crossref

44.

L. S.

Norris

,

M. A.

Ratner

,

A. E.

Roitberg

, and

R. B.

Gerber

,

J. Chem. Phys.

105

,

11261

(

1996

).

https://doi.org/10.1063/1.472922

Google Scholar

Crossref

45.

J. M.

Bowman

,

J. Chem. Phys.

68

,

608

(

1978

).

https://doi.org/10.1063/1.435782

Google Scholar

Crossref

46.

M. A.

Ratner

and

R. B.

Gerber

,

J. Phys. Chem.

90

,

20

(

1986

).

https://doi.org/10.1021/j100273a008

Google Scholar

Crossref

47.

J. M.

Bowman

,

Acc. Chem. Res.

19

,

202

(

1986

).

https://doi.org/10.1021/ar00127a002

Google Scholar

Crossref

48.

S.

Hirata

,

M.

Keçeli

, and

K.

Yagi

,

J. Chem. Phys.

133

,

034109

(

2010

).

https://doi.org/10.1063/1.3462237

Google Scholar

Crossref

PubMed

49.

M.

Keçeli

and

S.

Hirata

,

J. Chem. Phys.

135

,

134108

(

2011

).

https://doi.org/10.1063/1.3644895

Google Scholar

Crossref

PubMed

50.

M. R.

Hermes

and

S.

Hirata

,

J. Phys. Chem. A

117

,

7179

(

2013

).

https://doi.org/10.1021/jp4008834

Google Scholar

Crossref

PubMed

51.

N.

Makri

,

J. Phys. Chem. B

103

,

2823

(

1999

).

https://doi.org/10.1021/jp9847540

Google Scholar

Crossref

52.

S.

Hirata

and

M. R.

Hermes

,

J. Chem. Phys.

141

,

184111

(

2014

).

https://doi.org/10.1063/1.4901061

Google Scholar

Crossref

PubMed

53.

D. J.

Hooton

,

Philos. Mag.

3

,

49

(

1958

).

https://doi.org/10.1080/14786435808243224

Google Scholar

Crossref

54.

T. R.

Koehler

,

Phys. Rev. Lett.

17

,

89

(

1966

).

https://doi.org/10.1103/PhysRevLett.17.89

Google Scholar

Crossref

55.

P.

Choquard

,

The Anharmonic Crystal

(

Benjamin

,

New York

,

1967

).

Google Scholar

56.

N. S.

Gillis

,

N. R.

Werthamer

, and

T. R.

Koehler

,

Phys. Rev.

165

,

951

(

1968

).

https://doi.org/10.1103/PhysRev.165.951

Google Scholar

Crossref

57.

T.

Luty

,

A.

van der Avoird

,

R. M.

Berns

, and

T.

Wasiutynski

,

J. Chem. Phys.

75

,

1451

(

1981

).

https://doi.org/10.1063/1.442152

Google Scholar

Crossref

58.

K.

Yagi

,

K.

Hirao

,

T.

Taketsugu

,

M. W.

Schmidt

, and

M. S.

Gordon

,

J. Chem. Phys.

121

,

1383

(

2004

).

https://doi.org/10.1063/1.1764501

Google Scholar

Crossref

PubMed

59.

M. R.

Hermes

and

S.

Hirata

,

J. Chem. Phys.

139

,

034111

(

2013

).

https://doi.org/10.1063/1.4813123

Google Scholar

Crossref

PubMed

60.

K.

Yagi

,

M.

Keçeli

, and

S.

Hirata

,

J. Chem. Phys.

137

,

204118

(

2012

).

https://doi.org/10.1063/1.4767776

Google Scholar

Crossref

PubMed

61.

B.

Thomsen

,

K.

Yagi

, and

O.

Christiansen

,

J. Chem. Phys.

140

,

034111

(

2014

).

https://doi.org/10.1063/1.4870775

Google Scholar

Crossref

PubMed

62.

X.

Cheng

and

R. P.

Steele

,

J. Chem. Phys.

141

,

104105

(

2014

).

https://doi.org/10.1063/1.4894507

Google Scholar

Crossref

PubMed

63.

G. C.

Wick

,

Phys. Rev.

80

,

268

(

1950

).

https://doi.org/10.1103/PhysRev.80.268

Google Scholar

Crossref

64.

M. R.

Hermes

and

S.

Hirata

,

J. Chem. Phys.

141

,

084105

(

2014

).

https://doi.org/10.1063/1.4892614

Google Scholar

Crossref

PubMed

65.

J.

Goldstone

,

Proc. R. Soc. A

239

,

267

(

1957

).

https://doi.org/10.1098/rspa.1957.0037

Google Scholar

Crossref

66.

See supplementary material at http://dx.doi.org/10.1063/1.4931472 for the higher-order XVCC and EOM-XVCC equations.

67.

K.

Emrich

,

Nucl. Phys. A

351

,

379

(

1981

).

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

Google Scholar

Crossref

68.

K.

Emrich

,

Nucl. Phys. A

351

,

397

(

1981

).

https://doi.org/10.1016/0375-9474(81)90180-9

Google Scholar

Crossref

69.

H.

Sekino

and

R. J.

Bartlett

,

Int. J. Quantum Chem.

S18

,

255

(

1984

).

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

Google Scholar

Crossref

70.

J.

Geertsen

,

M.

Rittby

, and

R. J.

Bartlett

,

Chem. Phys. Lett.

164

,

57

(

1989

).

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

Google Scholar

Crossref

71.

D. C.

Comeau

and

R. J.

Bartlett

,

Chem. Phys. Lett.

207

,

414

(

1993

).

https://doi.org/10.1016/0009-2614(93)89023-B

Google Scholar

Crossref

72.

J. F.

Stanton

and

R. J.

Bartlett

,

J. Chem. Phys.

98

,

7029

(

1993

).

https://doi.org/10.1063/1.464746

Google Scholar

Crossref

73.

H. J.

Monkhorst

,

Int. J. Quantum Chem.

12

(

S11

),

421

(

1977

).

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

Google Scholar

Crossref

74.

S.

Ghosh

,

D.

Mukherjee

, and

S.

Bhattacharyya

,

Mol. Phys.

43

,

173

(

1981

).

https://doi.org/10.1080/00268978100101261

Google Scholar

Crossref

75.

M.

Takahashi

and

J.

Paldus

,

J. Chem. Phys.

85

,

1486

(

1986

).

https://doi.org/10.1063/1.451241

Google Scholar

Crossref

76.

H.

Koch

and

P.

Jørgensen

,

J. Chem. Phys.

93

,

3333

(

1990

).

https://doi.org/10.1063/1.458814

Google Scholar

Crossref

77.

H.

Koch

,

H. J. A.

Jensen

,

P.

Jørgensen

, and

T.

Helgaker

,

J. Chem. Phys.

93

,

3345

(

1990

).

https://doi.org/10.1063/1.458815

Google Scholar

Crossref

78.

R. J.

Rico

and

M.

Head-Gordon

,

Chem. Phys. Lett.

213

,

224

(

1993

).

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

Google Scholar

Crossref

79.

D.

Knoll

and

D.

Keyes

,

J. Comput. Phys.

193

,

357

(

2004

).

https://doi.org/10.1016/j.jcp.2003.08.010

Google Scholar

Crossref

80.

C.

Lanczos

,

J. Res. Natl. Bur. Stand.

45

,

255

(

1950

).

https://doi.org/10.6028/jres.045.026

Google Scholar

Crossref

81.

E. R.

Davidson

,

J. Comput. Phys.

17

,

87

(

1975

).

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

Google Scholar

Crossref

82.

K.

Hirao

and

H.

Nakatsuji

,

J. Comput. Phys.

45

,

246

(

1982

).

https://doi.org/10.1016/0021-9991(82)90119-X

Google Scholar

Crossref

83.

D.

Sorensen

,

R.

Lehoucq

,

C.

Yang

, and

K.

Maschhoff

,

arpack

(

Rice University

,

1996-2008

).

Google Scholar

84.

S.

Hirata

,

J. Phys. Chem. A

107

,

9887

(

2003

).

https://doi.org/10.1021/jp034596z

Google Scholar

Crossref

85.

T.

Shiozaki

,

M.

Kamiya

,

S.

Hirata

, and

E. F.

Valeev

,

Phys. Chem. Chem. Phys.

10

,

3358

(

2008

).

https://doi.org/10.1039/b803704n

Google Scholar

Crossref

PubMed

86.

M.

Valiev

,

E. J.

Bylaska

,

N.

Govind

,

K.

Kowalski

,

T. P.

Straatsma

,

H. J. J.

Van Dam

,

D.

Wang

,

J.

Nieplocha

,

E.

Apra

,

T. L.

Windus

, and

W. A.

de Jong

,

Comput. Phys. Commun.

181

,

1477

(

2010

).

https://doi.org/10.1016/j.cpc.2010.04.018

Google Scholar

Crossref

87.

K.

Yagi

,

sindo

(

The University of Tokyo

,

2006

).

Google Scholar

88.

M.

Keçeli

,

mavi

(

University of Illinois at Urbana-Champaign

,

2011

).

Google Scholar

89.

S.

Hirata

,

M.

Nooijen

, and

R. J.

Bartlett

,

Chem. Phys. Lett.

326

,

255

(

2000

).

https://doi.org/10.1016/S0009-2614(00)00772-7

Google Scholar

Crossref

90.

M. H.

Beck

,

A.

Jäckle

,

G. A.

Worth

, and

H. D.

Meyer

,

Phys. Rep.

324

,

1

(

2000

).

https://doi.org/10.1016/S0370-1573(99)00047-2

Google Scholar

Crossref

91.

S.

Hirata

and

I.

Grabowski

,

Theor. Chem. Acc.

133

,

1440

(

2014

).

https://doi.org/10.1007/s00214-013-1440-y

Google Scholar

Crossref

92.

A. A.

Auer

,

G.

Baumgartner

,

D. E.

Bernholdt

,

A.

Bibireata

,

V.

Choppella

,

D.

Cociorva

,

X. Y.

Gao

,

R.

Harrison

,

S.

Krishnamoorthy

,

S.

Krishnan

,

C. C.

Lam

,

Q. D.

Lu

,

M.

Nooijen

,

R.

Pitzer

,

J.

Ramanujam

,

P.

Sadayappan

, and

A.

Sibiryakov

,

Mol. Phys.

104

,

211

(

2006

).

https://doi.org/10.1080/00268970500275780

Google Scholar

Crossref

93.

A.

Hartono

,

Q. D.

Lu

,

T.

Henretty

,

S.

Krishnamoorthy

,

H. J.

Zhang

,

G.

Baumgartner

,

D. E.

Bernholdt

,

M.

Nooijen

,

R.

Pitzer

,

J.

Ramanujam

, and

P.

Sadayappan

,

J. Phys. Chem. A

113

,

12715

(

2009

).

https://doi.org/10.1021/jp9051215

Google Scholar

Crossref

PubMed

2015

AIP Publishing LLC

Higher-order diagrammatic vibrational coupled-cluster theory

I. INTRODUCTION

II. THEORY

A. XVCC formalism

B. Algebraic derivation of XVCC

C. Diagrammatic derivation of XVCC

D. EOM-XVCC formalism

E. Algebraic derivation of EOM-XVCC

F. Diagrammatic derivation of EOM-XVCC

III. AUTOMATED COMPUTER IMPLEMENTATION

A. Strength reduction

B. Intermediate reuse

C. Factorization

D. Observed computational cost

IV. NUMERICAL TESTS

A. Water

B. Formaldehyde

V. CONCLUSION

Acknowledgments

APPENDIX A: RIGOROUS DERIVATION OF THE CONNECTED FORM OF THE XVCC EQUATIONS

APPENDIX B: XVCCSDT

APPENDIX C: EOM-XVCCSD

APPENDIX D: COST SCALING OF XVCCm AND EOM-XVCCm IN AN nTH-ORDER PES

REFERENCES

Supplementary Material

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

(1)	Draw an ${\hat{H}}_{N}$ vertex having up to n lines.
(2)	Draw 0 to $n {\hat{T}}_{k}$ vertices (1 ≤ k ≤ m) beneath the ${\hat{H}}_{N}$ vertex.
(3)	Connect the lines of the ${\hat{T}}_{k}$ vertices with the lines of the ${\hat{H}}_{N}$
	vertex without altering (upgoing or downgoing) line directions,
	leaving l external (i.e., dangling) lines to form a connected
	diagram, which is either closed (l = 0) or open (l > 0).
(4)	Repeat steps (1)–(3) to enumerate all topologically distinct
	diagrams. If permutation of $\hat{T}$ vertices allows two diagrams
	to overlap, they are deemed topologically equivalent.

(1)	Draw an ${\hat{H}}_{N}$ vertex having up to n lines.
(2)	Draw a ${\hat{R}}_{h}$ vertex (1 ≤ h ≤ m) and zero to $(n - 1) {\hat{T}}_{k}$ vertices
	(1 ≤ k ≤ m) beneath the ${\hat{H}}_{N}$ vertex.
(3)	Connect the lines of the ${\hat{R}}_{h}$ and ${\hat{T}}_{k}$ vertices with the lines of
	the ${\hat{H}}_{N}$ vertex without altering (upgoing or downgoing) line
	directions, leaving l external (i.e., dangling) lines to form
	a connected diagram.
(4)	Repeat steps (1)–(3) to enumerate all topologically distinct
	diagrams. If permutation of $\hat{T}$ vertices allows two
	diagrams to overlap, they are deemed topologically equivalent.

$T_{2}^{(3)}$	≡	$T_{2}^{(2)} + \hat{P} (ν μ) (\frac{1}{2!} W_{i} τ_{i ν μ} + \frac{1}{2!} W_{ν i j} τ_{i j μ} + \frac{1}{2!} W_{i j}^{-} τ_{i} τ_{j ν μ}$
		$+ \frac{1}{2!} W_{ν i j k} τ_{i} τ_{j k μ} + \frac{1}{2!} W_{i j k} τ_{i ν} τ_{j k μ} + \frac{1}{2!^{2}} W_{i j k} τ_{i j} τ_{k ν μ}$
		$+ \frac{1}{2!^{2}} W_{i j k} τ_{i} τ_{j} τ_{k ν μ} + \frac{1}{2! 3!} W_{i j k l} τ_{i j k} τ_{l ν μ} + \frac{1}{2!^{3}} W_{i j k l} τ_{i j ν} τ_{k l μ}$
		$+ \frac{1}{2!^{2}} W_{i j k l} τ_{i} τ_{j k} τ_{l ν μ} + \frac{1}{2!} W_{i j k l} τ_{i} τ_{j ν} τ_{k l μ} + \frac{1}{2! 3!} W_{i j k l} τ_{i} τ_{j} τ_{k} τ_{l ν μ})$
	=	0

$T_{3}^{(3)}$	≡	$\hat{P} (ν μ λ) (\frac{1}{3!} W_{ν μ λ} + \frac{1}{3!} W_{ν μ λ i} τ_{i} + \frac{1}{2!} W_{ν μ i} τ_{i λ} + \frac{1}{2!} W_{ν i}^{+} τ_{i μ λ}$
		$+ \frac{1}{2!^{2}} W_{ν μ i j} τ_{i j λ} + \frac{1}{2!} W_{ν μ i j} τ_{i} τ_{j λ} + \frac{1}{2!} W_{ν i j} τ_{i μ} τ_{j λ} + \frac{1}{2!} W_{ν i j} τ_{i} τ_{j μ λ}$
		$+ \frac{1}{2!} W_{i j}^{-} τ_{i ν} τ_{j μ λ} + \frac{1}{2!} W_{ν i j k} τ_{i μ} τ_{j k λ} + \frac{1}{2!^{2}} W_{ν i j k} τ_{i j} τ_{k μ λ}$
		$+ \frac{1}{2!^{2}} W_{i j k} τ_{i j ν} τ_{k μ λ} + \frac{1}{2!} W_{ν i j k} τ_{i} τ_{j μ} τ_{k λ} + \frac{1}{3!} W_{i j k} τ_{i ν} τ_{j μ} τ_{k λ}$
		$+ \frac{1}{2!^{2}} W_{ν i j k} τ_{i} τ_{j} τ_{k μ λ} + \frac{1}{2!} W_{i j k} τ_{i} τ_{j ν} τ_{k μ λ} + \frac{1}{2!^{2}} W_{i j k l} τ_{i} τ_{j k ν} τ_{l μ λ}$
		$+ \frac{1}{2!^{2}} W_{i j k l} τ_{i ν} τ_{j μ} τ_{k l λ} + \frac{1}{2!^{2}} W_{i j k l} τ_{i j} τ_{k ν} τ_{l μ λ} + \frac{1}{3!} W_{i j k l} τ_{i} τ_{j ν} τ_{k μ} τ_{l λ}$
		$+ \frac{1}{2!^{2}} W_{i j k l} τ_{i} τ_{j} τ_{k ν} τ_{l μ λ})$
	=	0

$R_{1}^{(2)}$	≡	$R_{1}^{(1)} + W_{i} r_{i ν} + \frac{1}{2!} W_{ν i j} r_{i j} + W_{i j}^{-} τ_{j ν} r_{i} + W_{i j}^{-} τ_{j} r_{i ν}$
		$+ \frac{1}{2!} W_{ν i j k} τ_{j k} r_{i} + \frac{1}{2!} W_{ν i j k} τ_{k} r_{i j} + \frac{1}{2!} W_{i j k} τ_{j k} r_{i ν} + \frac{1}{2!} W_{i j k} τ_{k ν} r_{i j}$
		$+ W_{i j k} τ_{j} τ_{k ν} r_{i} + \frac{1}{2!} W_{i j k} τ_{j} τ_{k} r_{i ν} + \frac{1}{2!} W_{i j k l} τ_{j k} τ_{l ν} r_{i}$
		$+ \frac{1}{2!} W_{i j k l} τ_{k} τ_{l ν} r_{i j} + \frac{1}{2!} W_{i j k l} τ_{j} τ_{k l} r_{i ν} + \frac{1}{2!} W_{i j k l} τ_{j} τ_{k} τ_{l ν} r_{i}$
		$+ \frac{1}{3!} W_{i j k l} τ_{j} τ_{k} τ_{l} r_{i ν}$
	=	$ω_{k}^{(2)} r_{ν}$

$R_{2}^{(2)}$	≡	$\hat{P} (ν μ) (\frac{1}{2!} W_{ν μ i} r_{i} + W_{ν i}^{+} r_{i μ} + \frac{1}{2!^{2}} W_{ν μ i j} r_{i j} + \frac{1}{2!} W_{ν μ i j} τ_{j} r_{i}$
		$+ W_{ν i j} τ_{j μ} r_{i} + W_{ν i j} τ_{j} r_{i μ} + W_{i j}^{-} τ_{j μ} r_{i ν} + \frac{1}{2!} W_{ν i j k} τ_{j k} r_{i μ}$
		$+ \frac{1}{2!} W_{ν i j k} τ_{k μ} r_{i j} + W_{ν i j k} τ_{j} τ_{k μ} r_{i} + \frac{1}{2!} W_{i j k} τ_{j ν} τ_{k μ} r_{i}$
		$+ \frac{1}{2!} W_{ν i j k} τ_{j} τ_{k} r_{i μ} + W_{i j k} τ_{j} τ_{k μ} r_{i ν} + \frac{1}{2!} W_{i j k l} τ_{j k} τ_{l μ} r_{i ν}$
		$+ \frac{1}{2!^{2}} W_{i j k l} τ_{k ν} τ_{l μ} r_{i j} + \frac{1}{2!} W_{i j k l} τ_{j} τ_{k ν} τ_{l μ} r_{i} + \frac{1}{2!} W_{i j k l} τ_{j} τ_{k} τ_{l μ} r_{i ν})$
	=	$ω_{k}^{(2)} r_{ν μ}$

Higher-order diagrammatic vibrational coupled-cluster theory

I. INTRODUCTION

II. THEORY

A. XVCC formalism

B. Algebraic derivation of XVCC

C. Diagrammatic derivation of XVCC

D. EOM-XVCC formalism

E. Algebraic derivation of EOM-XVCC

F. Diagrammatic derivation of EOM-XVCC

III. AUTOMATED COMPUTER IMPLEMENTATION

A. Strength reduction

B. Intermediate reuse

C. Factorization

D. Observed computational cost

IV. NUMERICAL TESTS

A. Water

B. Formaldehyde

V. CONCLUSION

Acknowledgments

APPENDIX A: RIGOROUS DERIVATION OF THE CONNECTED FORM OF THE XVCC EQUATIONS

APPENDIX B: XVCCSDT

APPENDIX C: EOM-XVCCSD

APPENDIX D: COST SCALING OF XVCCm AND EOM-XVCCm IN AN nTH-ORDER PES

REFERENCES

Supplementary Material

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only