Coupled cluster Green function: Model involving single and double excitations

Many-body Green function (GF) techniques are powerful methods for addressing problems where strong correlation and excited states dictate key properties of chemical and material systems. These methods were introduced to chemical, condensed matter, and materials sciences’ communities via electronic structure methods where non-spectral GFs permit hierarchical embedding of physics and chemistry to address beyond-groundstate phenomena. Many-body GF approaches are^1–24 ideal methods because they provide a direct way of calculating electronic properties, excited states, electron transport, and equally important a direct procedure for calculating response and correlation functions, which allow for direct comparison with experiment. Of particular importance to chemistry, condensed matter physics and materials science is the accurate description of excited states that is necessary for explaining the observed properties in a number of technologically important areas including light harvesting materials,^25–31 semiconducting heterostructures,^32–34 colossal magnetoresistance, meta-materials, and novel superconductivity.^35–37 

A crucial outcome of our formulation is the determination of a first principles frequency dependent self-energy Σ(ω) that does not contain any adjustable parameters, and equally important, it contains the complete excitation spectrum corresponding to excitations of (N − 1) particle and (N + 1) particle systems. The development of a frequency dependent self-energy based on a coupled cluster (CC)^38–42 Green function that uses single and double excitations (GFCCSD) would provide the needed information for embedding single-particle formulations of the Green function. Our formulation provides an energy dependent GFCCSD matrix that can be constructed for any particular energy window. This approach should be viewed as a stepping-stone towards more accurate formulations along with several high performance algorithmic developments. For example, the efficiency of pole-shift algorithms in linear-response CC theory has been generalized to the Green function formulation (GFCC) approach,⁴³ which provides a simple, direct, and numerically efficient way of introducing a class of higher-rank correlation effects into the approximate GFCC approaches (such as the GFCCSD formulation) in a systematic and controllable way.

There are several advantages in using the GFCC to generate the self-energy.^43,44 First, the self-energy is obtained from first principles where it implicitly contains the level of theory of the reference state and the associated potentials (i.e., Hartree, Hartree-Fock exchange, or if density functional theory (DFT), v_xc, etc.), which provides a direct path for treating double counting terms. In addition, this approach naturally contains electron-electron, electron-hole, and hole-hole interactions and preserves the spin properties of the reference function, thereby, eliminating the need for any further approximations for determining Σ(ω). Second, the CCSD methods correspond to an infinite resummation of several classes of many-body perturbation theory (MBPT) contributions that also contain dispersive interactions (e.g., van der Waals). Third, the computation of the GFCC is a strong scaling method where for a fixed problem size, N, the numerical mathematics can be continually parallelized producing a drastic reduction in the time-to-solution that is crucial for rapid throughput and for allowing the simulation of large systems (i.e., surfaces and interfaces) with large basis sets.^45–47 

A rapid evolution of computational platforms and efficient parallel tools offers the possibility of constructing Green functions methodologies based on the CC parameterization of the ground-state wavefunction,^38–41 which is currently one of the most accurate approaches for encapsulating electronic correlation. In contrast to other wavefunction-based approaches, it combines several important features that play a critical role in high-level characterization of electron correlation effects in many-electron systems. The CC equations, amplitudes, and energies are described by connected diagrams that provide correct scaling of the CC energies as a function of the number of correlated particles. Furthermore, the high-order excitations can be represented by products of cluster operators of lower rank. This feature, in particular, is consistent with MBPT analysis. Equally important, the exact energy/wavefunction limit exists, which corresponds to the inclusion of all possible excitations in the cluster operator. The inclusion of higher-and-higher collective excitations in the cluster operator provided the foundation for the hierarchical structure of CC approximations starting from the electron-pair (or CCD, CC with doubles) approximation and CC with singles and doubles (CCSD)⁴² to CC with singles, doubles, and triples (CCSDT)^48–50 approaches and higher-rank formulations.

Our main idea is to combine the effectiveness of CC approach in encapsulating correlation effects within the Green function framework where the self-energy can naturally be obtained. Although the GFCCs were introduced to the field of quantum chemistry by Nooijen and Snijders two decades ago^51–53 (see also Ref. 54), there has been no instance of utilizing CC Green function and/or CC self-energies in electronic structure studies (as opposed to closely related electron-affinities (EAs)/ionization potentials (IPs) equation-of-motion CC (EOMCC) formulations^{45–47,55–61}). Green function techniques are natural embedding type methods where the many-body self-energy obtained from a GFCCSD calculation can potentially be embedded in a lower-order density functional type methods.

The paper is organized as follows. In Sec. II, we give a brief description of the GFCC formalism. In Sec. III, we focus on the details of the GFCCSD formalism including intermediate steps in calculating GFCCSD and the corresponding self-energy. In Sec. IV, we discuss implementation details of the GFCCSD formalism. In Sec. V, we discuss numerical results for small but well-known benchmark systems including neon atom, water molecule, and nitrogen dimer. We also discuss the accuracies of the GFCCSD poles for DNA/RNA base, 1,4-naphthalenedione, benzoquinone, and anthracene. In particular, we compare the quality of the CCSD self-energy calculations with GW based approaches and experimentally measured ionization potentials.

In this section, we will give the basic tenets of the GFCC formalism introduced by Nooijen and Snijders^51–53 The GFCC formalism hinges upon the CC bi-variational formalism^62–64 utilizing different parametrizations of the bra (⁠$〈Ψ0(N)|$⁠) and ket (⁠$|Ψ0(N)〉$⁠) ground-state wavefunctions of a N-electron system,

The cluster (T) operator, energy E₀, and de-excitation (Λ) operator are determined from the standard CC equations that are solved in the following order:

where Q is the projection operator onto the subspace spanned by Slater determinants generated by the T operator when acting on the reference function |Φ〉. In the exact formulation, the T and Λ operators are represented as sums of their many-body components (T_n and Λ_n),

where N_e stands for the total number of correlated electrons. The T_n and Λ_n operators can be given by the following expressions:

where $ta1…ani1…in$ and $λi1…ina1…an$ are the amplitudes determining T and Λ operators. The indices i, j, k, … (i₁, i₂, …) and a, b, c, … (a₁, a₂, …) correspond to occupied and unoccupied spinorbitals, respectively. The a_p (⁠$ap†$⁠) operator is the annihilation (creation) operator for electron in the pth state. By employing the bi-variational approach, the Green function can be expressed as

By introducing the resolution of identity e^−Te^T = 1, the Green function can be rewritten as

where the similarity transformed operators $āp$⁠, $āp†$⁠, and $H̄$ are given by the equations,

Using the Campbell-Hausdorff formula

one can write an explicit form for the similarity transformed creation and annihilation operators $āp$ an $āp†$ as

Using the spectral resolution form of the $H̄$ operator in N − 1 and N + 1 particle Hilbert spaces, one can show that the spectral resolution of the corresponding CC Green function involves sums over allN − 1 and N + 1 particle electronic states. The same holds for approximate GFCC approaches, where sums run over all electronic states obtained in corresponding (defined by the excitation rank of the T and Λ operators) subspaces of the N − 1 and N + 1 particle Hilbert spaces. To evaluate the G_pq(ω) matrix elements in a numerically efficient way, two sets of intermediate operators X_p and Y_q defined on N − 1 and N + 1 particle Hilbert spaces are introduced as follows:

This leads to the following compact expression for the CC Green function matrix elements:

where the first term on the right-hand-side of the above equation is the IP term, and the second term is the EA term.

We will now establish a direct connection between Eq. (20) and the following form for the spectral resolution:

where $Eχ(N±1)$ and $|Ψχ(N±1)〉$ denote energies and states of N ± 1 electron systems, and M_N+1 and M_N−1 refer to the total number of electronic states in the N + 1 and N − 1 electron spaces. The analogous form of the representation can be derived by employing spectral representations of the similarity transformed Hamiltonian $H̄$ in terms of the right $|Rχ(N±1)〉$ and left $〈Lχ(N±1)|$ eigenvectors of $H̄$ in N ± 1 Hilbert spaces (defined by the projection operators Q^(N±1)),

Introducing these representations into the formal solutions of Eqs. (18) and (19) naturally leads to the sum-over-states (SOS) form of the CC Green function,

which leads to the SOS form of the CC Green function (20), given by the following form:

It should be stressed that by solving the set of linear equations (18) and (19), the X_p(ω) and Y_q(ω) operators account for all poles. The above expression becomes identical with the spectral form of the Green function (21), where $〈Ψ0(N)|$⁠, $|Ψ0(N)〉$⁠, $〈Ψμ(N−1)|$⁠, $|Ψμ(N−1)〉$⁠, $〈Ψν(N+1)|$⁠, and $|Ψν(N+1)〉$ are replaced by the corresponding wavefunctions for N, N − 1, and N + 1 particle systems in the CC, IP-EOMCC, and EA-EOMCC parameterization, respectively, i.e.,

where the operators $Lμ(N−1)$ and $Rμ(N−1)$ (operators $Lν(N+1)$ and $Rν(N+1)$ are defined in an analogous way) are defined as

Since there are no clear criteria for how many states are needed in the SOS expansion, the SOS expansion is rather of limited use from the point of view of practical applications. Our approach has two clear advantages: (1) it allows to include higher-order correlation effects in a straightforward way, and (2) it does not require the diagonalization of large similarity transformed Hamiltonian matrix to obtain excited-state eigenvalues and eigenvectors necessary for calculating the SOS expansion. Similar problems can be encountered in closely related linear-response CC theory. Moreover, the choice of states to be included in the SOS expansion may strongly dependent on the ω-values of interest. By introducing intermediate X_p(ω) and Y_q(ω) operators, one can by-pass the above mentioned problems and have a formulation of the CC Green function method where all poles are included analytically, and no further adjustments are needed for various ω values.

In practical applications, the approximate forms of the T, Λ, X_p(ω), and Y_q(ω) operators need to be considered. A hierarchy of approximations based on inclusion of increasing excitation level in the cluster operator T has evolved into a standard procedure for constructing hierarchical approaches for inclusion of various classes of correlation effects and to control the accuracy of the resulting approximations. We will follow this procedure and define approximate GFCC methods based on the expansions truncated at certain excitation level n (n ≪ N_e) in the basic operators, which are required to determine the Green function matrix. For example, the n = 2 case corresponds to the GFCCSD formalism, while n = 3 produces the higher order GFCC with singles, doubles, and triples (GFCCSDT) method, etc. Other possibilities for inclusion of higher order-effects can be considered by using perturbation theory, for instance, one can use perturbative estimates of three body components of T, Λ, X_p(ω), Y_q(ω) operators based on converged singly and doubly excited amplitudes. In our recent paper, we have introduced the pole-shift procedure to include selected correlation effects due to triples in the GFCCSD formalism.⁴³ As in the single-reference CC calculations for the ground-state energies, the cost of the approximation grows dramatically with the excitation rank included in the corresponding operators.

In this paper, we will discuss the GFCCSD approximation and its numerical implementation. In particular, we will focus our attention on the GFCCSD implementation for closed-shell systems where the reference function |Φ〉 corresponds to the restricted Hartree-Fock (RHF) Slater determinant. In the GFCCSD method, we invoke the following form of the T, Λ, X_p(ω), and Y_q(ω) operators:

where the tensors $tai$⁠, $tabij$⁠, $λia$⁠, $λijab$⁠, xⁱ(ω)_p, $xaij(ω)p$⁠, y_a(ω)_q, and $yabi(ω)q$ are the amplitudes determining the T, Λ, X_p(ω), and Y_q(ω) operators. The algebraic forms of X_p(ω) and Y_q(ω) operators are analogous to the form of the R_K operators used in standard IP/EA-EOMCCSD formulations, respectively. The process of calculating the CCSD Green function is composed of four distinct computational steps:

The CCSD T and Λ operators are determined by solving the left and right CC Equations (3) and (5) using standard formulations,

where $E0(N)$ corresponds to the ground-state CCSD energy $E0CCSD$⁠, and the projection operators $Q1(N)$ and $Q2(N)$ are defined as

The singly and doubly excited configurations $|Φia〉$ and $|Φijab〉$ are given by the formulas,

Since the algebraic forms of the above equations have been amply discussed in the literature,^43,65 we will not further discuss the structure of these equations.

In analogy to the standard IP-EOMCCSD parameterization (see Eq. (34)), the X_p(ω) operator is represented as

where the xⁱ(ω) and $xaij(ω)$ amplitudes are determined by solving the linear set of equations in the N − 1 particle space

where the projection operators $Q1(N−1)$ and $Q2(N−1)$ are defined as follows:

The configurations |Φ_i〉 and $|Φija〉$ are defined by the action of the string of creation/annihilation operators onto the reference function |Φ〉, i.e.,

In the CCSD approximation, the $āp$ operator in Eq. (46) is given by the expression

which can be written in a more explicit form as

where O and V designate sets of occupied and unoccupied spinorbitals, respectively. The algebraic form of the $H̄Xp(ω)|Φ〉$ term in Eq. (46) is identical with the energy-independent terms corresponding to the standard IP-EOMCCSD equations, which have been extensively discussed in Ref. 55.

Following the same strategy as for the IP contribution to GFCCSD, the Y_q(ω) operator is represented in a way analogous to the standard EA-EOMCCSD parameterization (see Eq. (35)),

where the y_a(ω) and $yabi(ω)$ amplitudes are determined by solving the linear set of equations in the N + 1 particle space

with the projection operators $Q1(N+1)$ and $Q2(N+1)$ cast in the following form:

and the configurations |Φ^a〉 and $|Φiab〉$ are defined as

In the CCSD approximation, the $āq†$ operator in Eq. (55) is given by the expression,

or equivalently by

Again, the algebraic form of the $H̄Yq(ω)|Φ〉$ term in Eq. (55) is identical with the energy-independent terms corresponding to the standard EA-EOMCCSD equations (for more details, see Ref. 56).

Once the X_p(ω) and Y_q(ω) operators are determined, the corresponding contributions to the CCSD Green function matrix can be calculated from the equation,

An algebraic form for the above matrix elements can be derived by applying Wick’s theorem, which results in the following expressions for the EA and IP contributions:

and

where the X_p(ω)-dependent terms contribute to the pth row of the CCSD Green function matrix, and the Y_q(ω)-dependent terms contribute to the qth row of the CCSD Green function matrix. The equations for X_p(ω)’s and Y_q(ω)’s are decoupled and they can be calculated independently, thereby providing a highly parallel method for calculating the CCSD Green function. For example, one can envision the algorithm based on coarse grain parallelism where all X_p’s and Y_q’s are calculated independently by using a processor group strategy, which has been successfully tested in the context of multi-reference coupled cluster methods representing a conglomerate of loosely coupled set of nonlinear equations.^66–71 The workflow chart for calculating CCSD Green function is schematically shown in Figure 1. In this algorithm, each process associated with calculating X_p(ω) (p = 1, …, N) or Y_q(ω) (q = 1, …, N) is delegated to a separate processor group, which has access to basic tensors corresponding to 1- and 2-electron integrals defining the Hamiltonian matrix along with T₁, T₂, Λ₁, and Λ₂ operators (the so-called world group of Refs. 72–75).

The present algorithm allows one to include analytically all poles corresponding to the 1h/2h − 1p (IP) and 1p/1h − 2p (EA) subspaces of the N − 1 and N + 1 particle spaces and thus provides a crucial advantage over SOS type expansion based on inclusion of a small number of electronic states obtained from Eq. (26) via diagonalization of $H̄$ in 1h/2h − 1p (IP) and 1p/1h − 2p (EA) subspaces. There is also a practical aspect of our formalism associated with the possibility of obtaining Green function for ω values falling into an arbitrary energy window. In an SOS-type expansion, two problems are of particular concern: (1) what N − 1 and N + 1 particle states should be selected for the sum-over-state expansion? and (2) how to numerically obtain these states? Based on our experience with various types of EOMCC methodologies, obtaining states with excitation energies falling into a certain energy window (not necessarily the valence one) is very problematic for small molecules and due to high density of excited states unfeasible for larger systems. Our procedure replaces this problem by the set of 2 × N linear equations, which are much easier to solve than the EOMCC eigenvalue problems without a priori knowledge of excited states that needs to be included in the SOS expansion. Given the recent progress in the development of parallel implementations of the IP/EA-EOMCCSD implementations, the full form of the GFCCSD can be constructed for systems described by 1000-1500 orbitals (or even more).

The CCSD Green function formulation requires singly and doubly excited components of the X_p(ω) and Y_q(ω) operators, which are determined from Eqs. (46) and (55), respectively. In that case, the rank of excitation in the X_p(ω) and Y_q(ω) operators corresponds to the excitation ranks of T and Λ operators, which is reflected in the fact that only $(Q1(N±1)+Q2(N±1))H̄(Q1(N±1)+Q2(N±1))$ blocks of the CCSD similarity transformed Hamiltonian $H̄=e−T1−T2HeT1+T2$ are utilized in Eqs. (46) and (55). This procedure can be extended to the case when $H̄$ is represented in a larger space, for example, including single, double, and triple excitations (in fact, diagonalizing CCSD similarity transformed Hamiltonian in full N ± 1 particle spaces leads to the exact $Eμ(N±1)$ energies). In this case, $(Q1(N±1)+Q2(N±1)+Q3(N±1))H̄(Q1(N±1)+Q2(N±1)+Q3(N±1))$ blocks of $H̄$ are used to calculate X_p,i(ω) and Y_q,i(ω) (i = 1, 2, 3) operators. Subsequently, only X_p,i(ω) and Y_q,i(ω) (i = 1, 2) are utilized in constructing GFCCSD matrix (63). However, even the limited subset of amplitudes (singles and doubles) carries the information about the poles obtained through the diagonalization of $H̄$ in larger (than singles and doubles) subspaces of N ± 1 particle spaces. We refer to this approach as the GFCCSD with inner triples excitations (or GFCCSD-iT for short). In terms of complete SOS expansion, this procedure leads to an improved (beyond CCSD approximation) description of the N ± 1 electron wavefunctions $|RμN±1〉$ and $〈Lμ(N±1)|$ and the corresponding energies $Eμ(N±1)$⁠. This is an important feature of the GFCCSD formalism that in a natural way provides a path for going beyond 1h-2p and 2h-1p configurations and to include excitations of the 2h-3p and 3h-2p (or in principle any nh-(n + 1) p and (n + 1) h-np excitations) types to describe higher-order shake-up configurations in a non-perturbative way. The performance of the GFCCSD-iT methods will be discussed in the separate paper.

A key quantity of considerable interest is the self-energy (Σ(ω)) that describes many-body effects felt by an electron due to the interactions with the surrounding medium including all particles comprising the system. The self-energy Σ(ω) can be obtained from the definition of the Green function, whose representation is obtained via a Dyson type expansion,

G₀(ω) is the reference Green functions that is usually chosen to be the solution to the free-particle Helmholtz equation, and G(ω) is the CCSD Green function matrix. The Hartree-Fock Green function will be used for the G₀. This choice provides a better reference for assessing many body affects since the Hartree-Fock Green function contains the single electron interactions except correlation. The Hartree-Fock Green function is written as

where the n_p (⁠$n̄p=1−np$⁠) and ϵ_p’s are the occupancy of the HF orbitals and HF orbital energies, respectively.

The present parallel implementation of the GFCCSD approach has been developed using Tensor Contraction Engine functionalities^76,77 available in NWChem.⁷⁸ In this paper we discuss the GFCCSD implementation for the closed-shell systems where the reference wavefunction |Φ〉 is represented by the RHF Slater determinant. We also assume that the non-relativistic form of the Hamiltonian is used, which for closed-shell systems means the α − α and β − β diagonal blocks of the CCSD Green function matrix, G_αα and G_ββ, respectively, are equal, i.e.,

Our algorithm is based on existing CCSD implementations for determining CCSD T and Λ operators. Steps B and C of Sec. III were coded using recent EA/IP-EOMCCSD implementations in NWChem^45,46 that allow for the calculation of the most expensive terms of Eqs. (46) and (55) characterized by the action of $H̄$ operator onto X_p(ω)|Φ〉 (IP-EOMCCSD) and Y_q(ω)|Φ〉 (EA-EOMCCSD) vectors. We will symbolically designate these terms as A^IPx_p(ω) and A^EAy_q(ω) and rewrite Eqs. (46) and (55) in matrix form

where $bpIP$ and $bqEA$ correspond to free terms in Eqs. (46) and (55). The general solutions to the above equations are complex. For this reason, the x_p(ω) and y_q(ω) vectors are represented as sums of their real and imaginary parts,

The equation for the real parts of the x_p(ω) and y_q(ω) vectors can be easily decoupled from their corresponding imaginary counterparts. For example, for the IP problem (Eq. (71)), the $xpRe(ω)$ can be obtained from the equation

which can be solved using standard DIIS procedure.⁷⁹ Once $xpRe(ω)$ is known, the $xpIm(ω)$ can be obtained from the linear equation

A similar procedure can be applied for solving the EA problem (Eq. (72)) for the $yqRe(ω)$

and for the $yqIm(ω)$ vectors,

An important advantage of our approach is the fact that the calculations for the T and Λ operators need to be done only once,which results in the numerical cost of 2 × N⁶ (per iteration). For given η and ω, the numerical cost of solving of all N IP/EA linear equations (75)–(78) is 2N × N⁵ + 2N × N⁵ (per iteration), and the corresponding cost of forming the Green function matrix elements is proportional to 2N² × N⁴. However, one should not focus on the number of floating point operations since our approach is a strong scaling method that uses process groups to take advantage of the inherent independencies in the theory. This means as N grows, the numerical algorithm further decomposes the equations into processor groups, which means the real-time to solution goes from hours to minutes that allows for the simulation of larger systems with large basis sets.⁴⁵ The schematic representation of the CCSD Green function algorithm is shown in Fig. 1.

For illustrating some of the key features of the CCSD Green function, we have chosen some well-known benchmark systems including Ne, HF, H₂O, and N₂. A 3-21G and cc-pVDZ basis sets were used for all systems, and the equilibrium geometries were taken from Ref. 80 (HF), Ref. 81 (H₂O), and Ref. 80 (N₂). In all calculations, all electrons were correlated, and the imaginary part of η was set equal to zero.

The above benchmark systems are also aligned according to the increasing importance of the correlations effects. While for the Ne atom, the correlation effects are small, for the N₂ system, they play an important role, which is reflected by a significant increase in singly and doubly excited cluster amplitudes. To illustrate the accuracy of the ionization potentials of the GFCCSD (which correspond to the eigenvalues of the IP-EOMCCSD matrix), we have chosen several systems (uracil, thymine, cytosine, adenine, guanine) whose IPs have been intensively studied both theoretically and experimentally.^82–97 In these simulations, we employed aug-cc-pVTZ basis sets⁹⁸ and Sadlej PVTZ basis sets (denoted hereafter by POL1).⁹⁹ Finally, to illustrate the quality of higher order poles, we chose several benchmark systems (1,4-naphthalenedione, benzoquinone, and anthracene), which have recently been studied with various Green function/electron-propagator methods. In particular, we focus our attention on well-established diagonal third-order self-energy approximation (D3) (for the review of diagonal Green function approximations, see Ref. 100), outer valence Green function (OVGF),^2,4,104 the non-diagonal, renormalized, second-order approach (NR2),¹⁰⁵ two-particle hole Tamm-Dancoff approximation (2p-h TDA),⁴ the partial third-order quasiparticle method P3,¹⁰⁶ and the third-order algebraic diagrammatic construction method (ADC(3)).^3,4