In this paper, we report on the development of a parallel implementation of the coupled-cluster (CC) Green function formulation (GFCC) employing single and double excitations in the cluster operator (GFCCSD). A key aspect of this work is the determination of the frequency dependent self-energy, Σ(*ω*). The detailed description of the underlying algorithm is provided, including approximations used that preserve the pole structure of the full GFCCSD method, thereby reducing the computational costs while maintaining an accurate character of methodology. Furthermore, for systems with strong local correlation, our formulation reveals a diagonally dominate block structure where as the non-local correlation increases, the block size increases proportionally. To demonstrate the accuracy of our approach, several examples including calculations of ionization potentials for benchmark systems are presented and compared against experiment.

## I. INTRODUCTION

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,^{43} 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 _{x}c*, 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)^{42} 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.

## II. COUPLED CLUSTER GREEN FUNCTION APPROACH

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 ($\u3008\Psi 0(N)|$) and ket ($|\Psi 0(N)\u3009$) ground-state wavefunctions of a N-electron system,

The cluster (*T*) operator, energy *E*_{0}, 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*and Λ

_{n}_{n}operators can be given by the following expressions:

where $ta1\u2026ani1\u2026in$ and $\lambda i1\u2026ina1\u2026an$ are the amplitudes determining *T* and Λ operators. The indices *i*, *j*, *k*, … (*i*_{1}, *i*_{2}, …) and *a*, *b*, *c*, … (*a*_{1}, *a*_{2}, …) correspond to occupied and unoccupied spinorbitals, respectively. The *a _{p}* ($ap\u2020$) operator is the annihilation (creation) operator for electron in the

*p*th state. By employing the bi-variational approach, the Green function can be expressed as

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

where the similarity transformed operators $a\u0304p$, $a\u0304p\u2020$, and $H\u0304$ are given by the equations,

Using the Campbell-Hausdorff formula

one can write an explicit form for the similarity transformed creation and annihilation operators $a\u0304p$ an $a\u0304p\u2020$ as

Using the spectral resolution form of the $H\u0304$ 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 *all**N* − 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*and

_{p}*Y*defined on

_{q}*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\chi (N\xb11)$ and $|\Psi \chi (N\xb11)\u3009$ 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\u0304$ in terms of the right $|R\chi (N\xb11)\u3009$ and left $\u3008L\chi (N\xb11)|$ eigenvectors of $H\u0304$ 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 $\u3008\Psi 0(N)|$, $|\Psi 0(N)\u3009$, $\u3008\Psi \mu (N\u22121)|$, $|\Psi \mu (N\u22121)\u3009$, $\u3008\Psi \nu (N+1)|$, and $|\Psi \nu (N+1)\u3009$ 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\mu (N\u22121)$ and $R\mu (N\u22121)$ (operators $L\nu (N+1)$ and $R\nu (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.

## III. GFCCSD

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*) in the basic operators, which are required to determine the Green function matrix. For example, the

_{e}*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.

^{43}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$, $\lambda ia$, $\lambda ijab$, *x ^{i}*(

*ω*)

_{p}, $xaij(\omega )p$,

*y*(

_{a}*ω*)

_{q}, and $yabi(\omega )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*operators used in standard IP/EA-EOMCCSD formulations, respectively. The process of calculating the CCSD Green function is composed of four distinct computational steps:

_{K}### A. Solving CCSD equations for the *T* and Λ operators

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 $|\Phi ia\u3009$ and $|\Phi ijab\u3009$ 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.

### B. Solving for the intermediate ionization-potential-type operators *X*_{p}(*ω*) in the CCSD parameterization

_{p}

In analogy to the standard IP-EOMCCSD parameterization (see Eq. (34)), the *X _{p}*(

*ω*) operator is represented as

where the *x ^{i}*(

*ω*) and $xaij(\omega )$ amplitudes are determined by solving the linear set of equations in the

*N*− 1 particle space

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

The configurations |Φ_{i}〉 and $|\Phi ija\u3009$ are defined by the action of the string of creation/annihilation operators onto the reference function |Φ〉, i.e.,

In the CCSD approximation, the $a\u0304p$ 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\u0304Xp(\omega )|\Phi \u3009$ 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.

### C. Solving for the intermediate electron-affinity-type operators *Y*_{q}(*ω*) in the CCSD parameterization

_{q}

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(\omega )$ 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 $|\Phi iab\u3009$ are defined as

In the CCSD approximation, the $a\u0304q\u2020$ operator in Eq. (55) is given by the expression,

or equivalently by

### D. Forming elements of the CCSD Green function

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

*p*th row of the CCSD Green function matrix, and the

*Y*(

_{q}*ω*)-dependent terms contribute to the

*q*th 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*’s and

_{p}*Y*’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.

_{q}^{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*

_{1},

*T*

_{2}, Λ

_{1}, and Λ

_{2}operators (the so-called world group of Refs. 72–75).

The present algorithm allows one to include analytically *all poles* corresponding to the 1*h*/2*h* − 1*p* (IP) and 1*p*/1*h* − 2*p* (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\u0304$ in 1*h*/2*h* − 1*p* (IP) and 1*p*/1*h* − 2*p* (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\xb11)+Q2(N\xb11))H\u0304(Q1(N\xb11)+Q2(N\xb11))$ blocks of the CCSD similarity transformed Hamiltonian $H\u0304=e\u2212T1\u2212T2HeT1+T2$ are utilized in Eqs. (46) and (55). This procedure can be extended to the case when $H\u0304$ 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\mu (N\xb11)$ energies). In this case, $(Q1(N\xb11)+Q2(N\xb11)+Q3(N\xb11))H\u0304(Q1(N\xb11)+Q2(N\xb11)+Q3(N\xb11))$ blocks of $H\u0304$ 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\u0304$ 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\mu N\xb11\u3009$ and $\u3008L\mu (N\xb11)|$ and the corresponding energies $E\mu (N\xb11)$. This is an important feature of the GFCCSD formalism that in a natural way provides a path for going beyond 1

*h*-2

*p*and 2

*h*-1

*p*configurations and to include excitations of the 2

*h*-3

*p*and 3

*h*-2

*p*(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**_{0}(*ω*) 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**_{0}. 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\u0304p=1\u2212np$) and

*ϵ*’s are the occupancy of the HF orbitals and HF orbital energies, respectively.

_{p}## IV. IMPLEMENTATION DETAILS

The present parallel implementation of the GFCCSD approach has been developed using Tensor Contraction Engine functionalities^{76,77} available in NWChem.^{78} 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\u0304$ operator onto *X _{p}*(

*ω*)|Φ〉 (IP-EOMCCSD) and

*Y*(

_{q}*ω*)|Φ〉 (EA-EOMCCSD) vectors. We will symbolically designate these terms as

**A**

^{IP}

**x**

_{p}(

*ω*) and

**A**

^{EA}

**y**

_{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(\omega )$ can be obtained from the equation

which can be solved using standard DIIS procedure.^{79} Once $xpRe(\omega )$ is known, the $xpIm(\omega )$ can be obtained from the linear equation

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

and for the $yqIm(\omega )$ 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*^{6} (per iteration). For given *η* and *ω*, the numerical cost of solving of all *N* IP/EA linear equations (75)–(78) is 2*N* × *N*^{5} + 2*N* × *N*^{5} (per iteration), and the corresponding cost of forming the Green function matrix elements is proportional to 2*N*^{2} × *N*^{4}. 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.^{45} The schematic representation of the CCSD Green function algorithm is shown in Fig. 1.

## V. RESULTS AND DISCUSSION

For illustrating some of the key features of the CCSD Green function, we have chosen some well-known benchmark systems including Ne, HF, H_{2}O, and N_{2}. 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_{2}O), and Ref. 80 (N_{2}). 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_{2} 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^{98} and Sadlej PVTZ basis sets (denoted hereafter by POL1).^{99} 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),^{105} two-particle hole Tamm-Dancoff approximation (2p-h TDA),^{4} the partial third-order quasiparticle method P3,^{106} and the third-order algebraic diagrammatic construction method (ADC(3)).^{3,4}

In the IP-EOMCC/GFCCSD calculations, the geometries for uracil, thymine, cytosine, adenine, and guanine were optimized employing the POL1 basis set and B3LYP level of theory. In Table I we collected IPs calculated with various *ab initio* methods including GW, CASPT2, CCSD, CCSD(T), and IP-EOMCCSD/GFCCSD methodologies.

Method . | Uracil . | Thymine . | Cytosine . | Adenine . | Guanine . |
---|---|---|---|---|---|

G_{0}W_{0}(PBE)^{a} | 8.99 | 8.63 | 8.18 | 7.99 | 7.64 |

G_{0}W_{0}(LDA)^{a} | 9.03 | 8.64 | 8.21 | 7.90 | 7.49 |

GW(LDA)^{a} | 9.47 | 9.05 | 8.73 | 8.22 | 7.81 |

CASPT2/CASSCF/ANO-L 431/21^{b} | 9.22 | 8.87 | 8.54 | 8.18 | 7.91 |

CASPT2(IPEA=0.25)/ANO-L 431/21^{b} | 9.42 | 9.07 | 8.73 | 8.37 | 8.09 |

CCSD//CCSD/aug-cc-pVDZ^{b} | 9.40 | 9.01 | 8.72 | 8.39 | 8.07 |

CCSD(T)//CCSD/aug-cc-pVDZ^{b} | 9.43 | 9.04 | 8.76 | 8.40 | 8.09 |

IP-EOMCCSD//cc-pVTZ^{c} | … | 9.13 | 8.78 | 8.37 | 8.15 |

IP-EOMCCSD/cc-pVDZ | 9.16 | 8.78 | 8.34 | 7.97 | 7.65 |

IP-EOMCCSD/aug-cc-pVDZ | 9.43 | 9.02 | 8.66 | 8.26 | 7.97 |

IP-EOMCCSD/aug-cc-pVTZ | 9.63 | 9.22 | 8.86 | 8.45 | 8.15 |

IP-EOMCCSD/POL1 | 9.42 | 9.01 | 8.65 | 8.26 | 7.96 |

Expt. range | [9.4,9.6] | [9.0,9.2] | [8.8,9.0] | [8.3,8.5] | [8.0,8.3] |

Method . | Uracil . | Thymine . | Cytosine . | Adenine . | Guanine . |
---|---|---|---|---|---|

G_{0}W_{0}(PBE)^{a} | 8.99 | 8.63 | 8.18 | 7.99 | 7.64 |

G_{0}W_{0}(LDA)^{a} | 9.03 | 8.64 | 8.21 | 7.90 | 7.49 |

GW(LDA)^{a} | 9.47 | 9.05 | 8.73 | 8.22 | 7.81 |

CASPT2/CASSCF/ANO-L 431/21^{b} | 9.22 | 8.87 | 8.54 | 8.18 | 7.91 |

CASPT2(IPEA=0.25)/ANO-L 431/21^{b} | 9.42 | 9.07 | 8.73 | 8.37 | 8.09 |

CCSD//CCSD/aug-cc-pVDZ^{b} | 9.40 | 9.01 | 8.72 | 8.39 | 8.07 |

CCSD(T)//CCSD/aug-cc-pVDZ^{b} | 9.43 | 9.04 | 8.76 | 8.40 | 8.09 |

IP-EOMCCSD//cc-pVTZ^{c} | … | 9.13 | 8.78 | 8.37 | 8.15 |

IP-EOMCCSD/cc-pVDZ | 9.16 | 8.78 | 8.34 | 7.97 | 7.65 |

IP-EOMCCSD/aug-cc-pVDZ | 9.43 | 9.02 | 8.66 | 8.26 | 7.97 |

IP-EOMCCSD/aug-cc-pVTZ | 9.63 | 9.22 | 8.86 | 8.45 | 8.15 |

IP-EOMCCSD/POL1 | 9.42 | 9.01 | 8.65 | 8.26 | 7.96 |

Expt. range | [9.4,9.6] | [9.0,9.2] | [8.8,9.0] | [8.3,8.5] | [8.0,8.3] |

The G_{0}W_{0} results presented here use different numerical formulations. Marzari, *et al.*^{84} used an optimal plane-wave basis set for representing operators such as the polarizability, and they replace the sum over states approach for calculating the irreducible dynamic polarizability and self-energy^{107} by a Lanczos chain algorithm.^{108} Whereas, Faber, *et al.*^{109} use a localized basis set containing a Gaussian auxiliary basis combined with the radial part of the numerical basis, which is fit with up to 5 contracted Gaussian for calculating the Coulomb matrix elements and the matrix elements associated with the auxiliary Gaussian basis.^{109} In addition, it should be pointed out that the GW(local-density approximation (LDA)) is an G_{0}W_{0} calculation but with self-consistent eigenvalues as opposed to the standard one-shot approach.^{109} This was done to overcome the well-known problem of overscreening contained in the LDA eigenvalues that is known to lead to a significant underestimation of the band gaps, and it tends to produce rather small ionization energies.

The ionization potentials of nucleobases can be adequately described by the CASPT2, CC, IP-EOMCCSD/GFCCSD, and GW(LDA) methods. From the IP plot (see Figure 2) one can see that all methods display a nearly linear trend when going from uracil to adenine with G_{0}W_{0}(LDA) and G_{0}W_{0}(PBE), both underestimating the IPs by around 0.5 eV. From Figure 2, one can see that CCSD(T) IP values all fall within the experimental range except cytosine, where a very small difference can be observed. Both GW(LDA) and IP-EOMCCSD/GFCCSD are very close to falling within the experimental range. Interestingly, IP-EOMCCSD calculation using an aug-cc-pVTZ basis produces slightly better results than the IP-EOMCCSD/GFCCSD calculation that employs the POL1 basis set. The G_{0}W_{0}(LDA) and G_{0}W_{0}(PBE) yield very similar results for the IP; however, they display a significant difference in going from cytosine to adenine. Going from cytosine to adenine, both G_{0}W_{0}(PBE) and GW(LDA) change slope with G_{0}W_{0}(PBE) having a positive slope while GW(LDA) has a negative one. The remaining higher order quantum chemistry methods and G_{0}W_{0}(LDA) slopes remain virtually unchanged. It is also worth mentioning that the GW(LDA) approach provides better IP values as compared to experiment, but from cytosine to adenine, the values all fall below the experimental range.

In summary, the GFCCSD/EOMCCSD ionization potentials, which correspond to the poles of the CCSD Green function, are in a very good agreement with experimentally observed data. From the above discussion, the GFCCSD/EOMCCSD formulations can avoid large extent problems with balancing individual correlation effects for the ground states of neutral and ionized systems.

The accuracy of the approximate GF is contingent upon the accuracies of higher-order ionization potentials and electron affinities. While this subject has recently been studied by us in the context of higher-order IPs and EAs for acenes Ref. 46, for the purpose of this paper we decided to compare the GFCCSD IP results with ones obtained with popular Green function/electron-propagator methods D3, P3, OVGF, NR2, 2p-h TDA, ADC(2), G_{0}W_{0} + SOSEX@PBE for 1,4-naphthalenedione (Table II), benzoquinone (Table III), and anthracene (Table IV).

Method . | ^{2}B_{2}
. | ^{2}B_{1}
. | ^{2}A_{2}
. |
---|---|---|---|

D3^{a} | 10.90 | 9.89 | 9.85 |

OVGF^{a} | 10.09 | 9.84 | 9.79 |

NR2^{a} | 9.80 | 10.06 | 9.95 |

2p-h TDA^{a} | 8.77 | 9.46 | 9.38 |

ADC(3)^{a} | 10.37 | 9.93 | 9.90 |

G_{0}W_{0} + SOSEX@PBE^{a} | 9.54 | 9.75 | 9.73 |

IP-EOMCCSD (aug-cc-pVDZ) | 9.71 | 9.78 | 9.72 |

IP-EOMCCSD (aug-cc-pVTZ) | 9.93 | 9.92 | 9.86 |

Expt.^{b} | 9.60 | 9.90 | 9.77 |

Method . | ^{2}B_{2}
. | ^{2}B_{1}
. | ^{2}A_{2}
. |
---|---|---|---|

D3^{a} | 10.90 | 9.89 | 9.85 |

OVGF^{a} | 10.09 | 9.84 | 9.79 |

NR2^{a} | 9.80 | 10.06 | 9.95 |

2p-h TDA^{a} | 8.77 | 9.46 | 9.38 |

ADC(3)^{a} | 10.37 | 9.93 | 9.90 |

G_{0}W_{0} + SOSEX@PBE^{a} | 9.54 | 9.75 | 9.73 |

IP-EOMCCSD (aug-cc-pVDZ) | 9.71 | 9.78 | 9.72 |

IP-EOMCCSD (aug-cc-pVTZ) | 9.93 | 9.92 | 9.86 |

Expt.^{b} | 9.60 | 9.90 | 9.77 |

Method . | ^{2}B_{2u}
. | ^{2}B_{3g}
. | ^{2}B_{3u}
. | ^{2}B_{1g}
. |
---|---|---|---|---|

D3^{a} | 11.77 | 11.13 | 11.24 | 11.07 |

OVGF^{a} | 10.90 | 10.43 | 11.00 | 11.07 |

NR2^{a} | 10.35 | 10.13 | 11.17 | 11.34 |

2p-h TDA^{a} | 9.88 | 9.17 | 10.66 | 10.86 |

ADC(3)^{a} | 11.02 | 10.67 | 11.12 | 11.20 |

G_{0}W_{0} + SOSEX@PBE^{a} | 10.30 | 9.95 | 10.76 | 10.89 |

IP-EOMCCSD (aug-cc-pVDZ) | 10.32 | 10.04 | 11.06 | 11.14 |

IP-EOMCCSD (aug-cc-pVTZ) | 10.56 | 10.27 | 11.21 | 11.29 |

Expt.^{b} | 10.41 | 10.11 | 11.06 | 11.23 |

Method . | ^{2}B_{2u}
. | ^{2}B_{3g}
. | ^{2}B_{3u}
. | ^{2}B_{1g}
. |
---|---|---|---|---|

D3^{a} | 11.77 | 11.13 | 11.24 | 11.07 |

OVGF^{a} | 10.90 | 10.43 | 11.00 | 11.07 |

NR2^{a} | 10.35 | 10.13 | 11.17 | 11.34 |

2p-h TDA^{a} | 9.88 | 9.17 | 10.66 | 10.86 |

ADC(3)^{a} | 11.02 | 10.67 | 11.12 | 11.20 |

G_{0}W_{0} + SOSEX@PBE^{a} | 10.30 | 9.95 | 10.76 | 10.89 |

IP-EOMCCSD (aug-cc-pVDZ) | 10.32 | 10.04 | 11.06 | 11.14 |

IP-EOMCCSD (aug-cc-pVTZ) | 10.56 | 10.27 | 11.21 | 11.29 |

Expt.^{b} | 10.41 | 10.11 | 11.06 | 11.23 |

Method . | B_{2g}
. | B_{3g}
. | A
. _{u} |
---|---|---|---|

OVGF^{a} | 7.17 | 8.27 | 9.07 |

P3^{a} | 7.31 | 8.46 | 9.09 |

IP-EOMCCSD (aug-cc-pVDZ)^{b} | 7.28 | 8.41 | 9.23 |

IP-EOMCCSD (aug-cc-pVTZ)^{b} | 7.42 | 8.55 | 9.36 |

Expt.^{c} | 7.47 | 8.57 | 9.23 |

In Table II, we report on the tree lowest ionization potentials for the ^{2}*B*_{2}, ^{2}*B*_{1}, and ^{2}*A*_{2} symmetries of the 1-4 naphthalenedione. The results are shown in Table II. While the IP-EOMCCSD aug-cc-pVTZ results for the ^{2}*B*_{1} and ^{2}*A*_{2} symmetries lie very close to the experimental values, for the ^{2}*B*_{2} state, the IP-EOMCCSD/aug-cc-pVTZ ionization potential is characterized by ≃0.3 eV error. Similar behavior for the ^{2}*B*_{2} ionization potential is disclosed by the D3, OVGF, NR2, and ADC(3) approaches. For this particular state, the G_{0}W_{0}+SOSEX@PBE methodology produces almost perfect agreement with the experiment.

For the benzoquinone molecule, the IP-EOMCCSD the aug-cc-pVTZ results are characterized by 0.15 (^{2}*B*_{2u}), 0.16 (^{2}*B*_{3g}), 0.15 (^{2}*B*_{3u}), 0.06 eV of error with respect to the experimental values. It is interesting to notice that IP-EOMCCSD performs for this case better than the ADC(3) formalism. Moreover, the IP-EOMCCSD formalism is competitive to the G_{0}W_{0}+SOSEX@PBE approach. The efficiency of the GFCCSD formalism is clearly seen especially for the lowest ^{2}*B*_{2u} state, where D3, OVGF, and ADC(3) methodologies consistently overestimate the experimental value of 10.41 eV.

Our last example of anthracene highlights the ability of the GFCCSD formalism to accurately locate higher-order ionization potentials (see Table IV). The IP-EOMCCSD aug-cc-pVTZ errors with respect to the experiment are 0.05, 0.02, and 0.13 eV for the states of *B*_{2g}, *B*_{3g}, and *A _{u}* symmetry, respectively. For comparison, the corresponding OVGF and P3 errors are 0.3, 0.3, 0.16 eV and 0.16, 0.09, 0.14 eV, respectively. One should expect further improvements in the location of poles by the inclusion of triple excitations (for example, as offered by the GFCCSD-iT approach).

Several criteria can be established to determine the strength of the correlation effects in the system from the GFCCSD matrix. Understanding how the electron correlation effects are built into the GFCCSD matrix would play a crucial role in understanding how to derive approximate formulations of the GFCCSD matrix to reduce computational cost and thereby provide a path for extending the GFCCSD approach to larger molecular systems. Two possible ways to characterize the impact of the correlation effects on the form of the CCSD Green function matrix is either by calculating the root-mean square of the off-diagonal elements or via differences between diagonal elements of the GFCCSD and Hartree-Fock Green function (*G*_{0,ii})

As expected (see Table V), the maximum values of the off-diagonal elements in the GFCCSD matrix grow proportionally with the strength of the ground-state correlation effects. For example, the maximum off-diagonal element of GFCCSD matrix for N_{2} is 0.104, while the maximum for Ne atom is only 0.005. Similar trends are found in the Δ_{off−D} values in Table V.

System . | $maxGi\u2260j$ . | $\Delta off\u2212D=\u2211i\u2260jGij2$ . |
---|---|---|

Ne | 0.005 | 0.010 |

HF | 0.021 | 0.055 |

H_{2}O | 0.036 | 0.090 |

N_{2} | 0.104 | 0.177 |

System . | $maxGi\u2260j$ . | $\Delta off\u2212D=\u2211i\u2260jGij2$ . |
---|---|---|

Ne | 0.005 | 0.010 |

HF | 0.021 | 0.055 |

H_{2}O | 0.036 | 0.090 |

N_{2} | 0.104 | 0.177 |

In Tables VI and VII we show the differences between diagonal elements CCSD (**G**(*ω* = 0)) and the Hartree-Fock (**G**_{0}(*ω* = 0)) Green function matrices (Δ_{ii} values), Hartree-Fock orbital energies (*ϵ _{HF}*), and the eigenvalues

*σ*of

_{CCSD}**F**+

**Σ**(

*ω*= 0) matrix (where

**F**is the Fock matrix) and differences between corresponding

*ϵ*’s and

_{HF}*σ*’s for the Ne and HF systems. Again, the Δ

_{CCSD}_{ii}values grow with the correlation level of the ground-state wavefunction. For example, the largest Δ

_{ii}for Ne is 0.08 and while for N

_{2}, it is 0.35. Similar behavior is found for

*σ*

_{CCSD}−

*ϵ*

_{HF}values, especially for the virtual orbitals. It is also interesting to notice that for

*ω*= 0, the largest values of Δ

_{ii}are “localized” in the valence orbital region, which can be seen in Fig. 3.

Ne: . | Ne: . | Ne: . | Ne: . | HF: . | HF: . | HF: . | HF: . |
---|---|---|---|---|---|---|---|

Δ_{ii} = G − _{ii}G_{0;ii}
. | ϵ_{HF}(occ.)
. | σ_{CCSD}
. | ϵ_{HF} − σ_{CCSD}
. | Δ_{ii} = G − _{ii}G_{0;ii}
. | ϵ_{HF}(occ.)
. | σ_{CCSD}
. | ϵ_{HF} − σ_{CCSD}
. |

0.00 | −32.56(2) | −32.48 | −0.09 | 0.00 | −26.11(2) | −26.03 | −0.09 |

0.01 | −1.87(2) | −1.83 | −0.04 | 0.02 | −1.55(2) | −1.52 | −0.04 |

0.08 | −0.79(2) | −0.75 | −0.04 | 0.03 | −0.71(2) | −0.69 | −0.01 |

0.08 | −0.79(2) | −0.75 | −0.04 | 0.16 | −0.60(2) | −0.55 | −0.05 |

0.08 | −0.79(2) | −0.75 | −0.04 | 0.16 | −0.60(2) | −0.55 | −0.05 |

−0.00 | 2.69(0) | 2.67 | 0.02 | −0.03 | 0.27(0) | 0.27 | 0.00 |

−0.00 | 2.69(0) | 2.67 | 0.02 | −0.01 | 1.23(0) | 1.21 | 0.02 |

−0.00 | 2.69(0) | 2.67 | 0.02 | −0.01 | 2.23(0) | 2.21 | 0.03 |

−0.00 | 4.08(0) | 4.06 | 0.03 | −0.01 | 2.23(0) | 2.21 | 0.03 |

−0.01 | 2.42(0) | 2.37 | 0.05 | ||||

−0.00 | 3.69(0) | 3.65 | 0.04 |

Ne: . | Ne: . | Ne: . | Ne: . | HF: . | HF: . | HF: . | HF: . |
---|---|---|---|---|---|---|---|

Δ_{ii} = G − _{ii}G_{0;ii}
. | ϵ_{HF}(occ.)
. | σ_{CCSD}
. | ϵ_{HF} − σ_{CCSD}
. | Δ_{ii} = G − _{ii}G_{0;ii}
. | ϵ_{HF}(occ.)
. | σ_{CCSD}
. | ϵ_{HF} − σ_{CCSD}
. |

0.00 | −32.56(2) | −32.48 | −0.09 | 0.00 | −26.11(2) | −26.03 | −0.09 |

0.01 | −1.87(2) | −1.83 | −0.04 | 0.02 | −1.55(2) | −1.52 | −0.04 |

0.08 | −0.79(2) | −0.75 | −0.04 | 0.03 | −0.71(2) | −0.69 | −0.01 |

0.08 | −0.79(2) | −0.75 | −0.04 | 0.16 | −0.60(2) | −0.55 | −0.05 |

0.08 | −0.79(2) | −0.75 | −0.04 | 0.16 | −0.60(2) | −0.55 | −0.05 |

−0.00 | 2.69(0) | 2.67 | 0.02 | −0.03 | 0.27(0) | 0.27 | 0.00 |

−0.00 | 2.69(0) | 2.67 | 0.02 | −0.01 | 1.23(0) | 1.21 | 0.02 |

−0.00 | 2.69(0) | 2.67 | 0.02 | −0.01 | 2.23(0) | 2.21 | 0.03 |

−0.00 | 4.08(0) | 4.06 | 0.03 | −0.01 | 2.23(0) | 2.21 | 0.03 |

−0.01 | 2.42(0) | 2.37 | 0.05 | ||||

−0.00 | 3.69(0) | 3.65 | 0.04 |

H_{2}O:
. | H_{2}O:
. | H_{2}O:
. | H_{2}O:
. | N_{2}:
. | N_{2}:
. | N_{2}:
. | N_{2}:
. |
---|---|---|---|---|---|---|---|

Δ_{ii} = G − _{ii}G_{0;ii}
. | ϵ_{HF}(occ.)
. | σ
. _{CCSD} | ϵ_{HF} − σ_{CCSD}
. | Δ_{ii} = G − _{ii}G_{0;ii}
. | ϵ_{HF}(occ.)
. | σ
. _{CCSD} | ϵ_{HF} − σ_{CCSD}
. |

0.00 | −20.43(2) | −20.34 | −0.09 | 0.00 | −15.61(2) | −15.52 | −0.10 |

0.02 | −1.33(2) | −1.29 | −0.03 | 0.00 | −15.61(2) | −15.52 | −0.10 |

0.00 | −0.69(2) | −0.68 | −0.00 | 0.02 | −1.51(2) | −1.46 | −0.05 |

0.13 | −0.54(2) | −0.50 | −0.04 | 0.12 | −0.76(2) | −0.70 | −0.06 |

0.26 | −0.48(2) | −0.43 | −0.05 | 0.15 | −0.61(2) | −0.64 | 0.03 |

−0.14 | 0.26(0) | 0.25 | 0.01 | −0.07 | −0.61(2) | −0.64 | 0.03 |

−0.08 | 0.36(0) | 0.35 | 0.01 | −0.08 | −0.61(2) | −0.56 | −0.05 |

−0.01 | 1.19(0) | 1.17 | 0.02 | 0.35 | 0.18(0) | 0.19 | −0.01 |

−0.02 | 1.31(0) | 1.28 | 0.03 | 0.34 | 0.18(0) | 0.19 | −0.01 |

−0.01 | 1.78(0) | 1.75 | 0.03 | −0.02 | 0.75(0) | 0.74 | 0.01 |

−0.01 | 1.87(0) | 1.82 | 0.04 | −0.01 | 1.21(0) | 1.18 | 0.02 |

−0.02 | 2.02(0) | 1.96 | 0.06 | −0.03 | 1.28(0) | 1.23 | 0.05 |

−0.00 | 3.11(0) | 3.07 | 0.04 | −0.03 | 1.28(0) | 1.23 | 0.05 |

−0.03 | 1.42(0) | 1.37 | 0.05 | ||||

−0.03 | 1.42(0) | 1.37 | 0.05 | ||||

−0.02 | 1.47(0) | 1.44 | 0.03 | ||||

−0.01 | 1.85(0) | 1.82 | 0.03 | ||||

−0.00 | 2.56(0) | 2.53 | 0.03 |

H_{2}O:
. | H_{2}O:
. | H_{2}O:
. | H_{2}O:
. | N_{2}:
. | N_{2}:
. | N_{2}:
. | N_{2}:
. |
---|---|---|---|---|---|---|---|

Δ_{ii} = G − _{ii}G_{0;ii}
. | ϵ_{HF}(occ.)
. | σ
. _{CCSD} | ϵ_{HF} − σ_{CCSD}
. | Δ_{ii} = G − _{ii}G_{0;ii}
. | ϵ_{HF}(occ.)
. | σ
. _{CCSD} | ϵ_{HF} − σ_{CCSD}
. |

0.00 | −20.43(2) | −20.34 | −0.09 | 0.00 | −15.61(2) | −15.52 | −0.10 |

0.02 | −1.33(2) | −1.29 | −0.03 | 0.00 | −15.61(2) | −15.52 | −0.10 |

0.00 | −0.69(2) | −0.68 | −0.00 | 0.02 | −1.51(2) | −1.46 | −0.05 |

0.13 | −0.54(2) | −0.50 | −0.04 | 0.12 | −0.76(2) | −0.70 | −0.06 |

0.26 | −0.48(2) | −0.43 | −0.05 | 0.15 | −0.61(2) | −0.64 | 0.03 |

−0.14 | 0.26(0) | 0.25 | 0.01 | −0.07 | −0.61(2) | −0.64 | 0.03 |

−0.08 | 0.36(0) | 0.35 | 0.01 | −0.08 | −0.61(2) | −0.56 | −0.05 |

−0.01 | 1.19(0) | 1.17 | 0.02 | 0.35 | 0.18(0) | 0.19 | −0.01 |

−0.02 | 1.31(0) | 1.28 | 0.03 | 0.34 | 0.18(0) | 0.19 | −0.01 |

−0.01 | 1.78(0) | 1.75 | 0.03 | −0.02 | 0.75(0) | 0.74 | 0.01 |

−0.01 | 1.87(0) | 1.82 | 0.04 | −0.01 | 1.21(0) | 1.18 | 0.02 |

−0.02 | 2.02(0) | 1.96 | 0.06 | −0.03 | 1.28(0) | 1.23 | 0.05 |

−0.00 | 3.11(0) | 3.07 | 0.04 | −0.03 | 1.28(0) | 1.23 | 0.05 |

−0.03 | 1.42(0) | 1.37 | 0.05 | ||||

−0.03 | 1.42(0) | 1.37 | 0.05 | ||||

−0.02 | 1.47(0) | 1.44 | 0.03 | ||||

−0.01 | 1.85(0) | 1.82 | 0.03 | ||||

−0.00 | 2.56(0) | 2.53 | 0.03 |

The lack of strong correlation effects in Ne atoms results in small differences between diagonal elements. The largest discrepancies occur for matrix elements corresponding to the valence 2*p* orbitals. All differences corresponding to virtual orbitals are smaller than 0.005, which is a consequence of the fact that the Hartree-Fock determinant provides a good approximation to the exact wavefunction. For HF and H_{2}O systems, one can observe a gradual increase in Δ_{ii} values especially for the valence orbitals (or active orbitals) involved in forming bonds. For example, for the HF molecule, the active space is composed of HOMO and LUMO orbitals, or equivalently, a (1,1) model space, where in general (*n _{oa}*,

*n*) stands for the active space composed of

_{va}*n*(

_{oa}*n*) highest lying occupied (lowest lying unoccupied) orbitals. In the case of the H

_{ov}_{2}O molecule, the largest differences correspond to molecular orbitals contributing to the (2,2) active space. A similar pattern is observed for the N

_{2}molecule.

The development of cost-reduced versions of Green function approaches plays an important role in enabling these formalisms for large molecular systems. Several approximate strategies based on the inclusion of diagonal terms or on the resummation of several classes of diagrams in the perturbative expansion of the self-energy have been extensively discussed in the literature including expanding self-energies through low order (second/third) and diagonal and non-diagonal approximations. These strategies have extensively been discussed in Refs. 2, 4, and 100–106 by analyzing accuracies of approximate approaches such as mentioned earlier OVGF, 2p-h TDA, ADC(3), D3, P3, and NR2 methods. Our approach is however slightly different: using well-defined hierarchy of the correlation effects in the CC wavefunction ansatz, we approximate the form of the self-energy through imposing the approximations on the form of the corresponding Green function.

Since the largest discrepancies between diagonal elements of the HF and CCSD Green function matrices for all benchmark systems considered above occur for the valence orbitals (or active spinorbitals), it is reasonable to represent the CCSD Green function matrix, at least for the frequencies *ω* in the valence/active region, by its approximate block form (referred to as B-GFCCSD) composed of diagonal elements of the HF Green function for inactive orbitals and active-orbitals sub-block of the CCSD Green function matrix, i.e.,

where

and $\delta p(A)=1\u2212\delta p(I)$. Although the above approximation does not eliminate the expensive CCSD and Λ-CCSD steps, it does significantly reduce the number of intermediate steps associated with calculating *X _{p}* and

*Y*operators required to form B-GFCCSD (note: only those

_{q}*X*and

_{p}*Y*operators defined by the active

_{q}*p*and

*q*indices need to be calculated). Additionally, the B-GFCCSD approximation does not affect the location of poles of the corresponding Green function (at least for the poles corresponding to electronic states of

*N*+ 1 and

*N*− 1 particle systems of the symmetry of active orbitals)—the pole structure of the B-GFCCSD will be exactly the same as the pole structure of the full GFCCSD matrix. Similar ideas, although in the context of solving self-consistent Dyson equation for second-order representation of the self-energy, were explored in the Basis Generated by Lanczos method (BAGEL) approach of Refs. 115 and 116. In contrast, our B-GFCCSD formulation does not introduce a truncation of the 2

*h*− 1

*p*and 1

*h*− 2

*p*spaces in order to select the most important configurations as well as to reduce the numerical footprint of the underlying numerical procedure. The B-GFCCSD approach can also be applied in the context of the GFCCSD-iT extension.

Because of the block-diagonal form of the **G**^{B−CCSD}(*ω*) matrix, the corresponding self-energy matrix (Σ^{B−CCSD}(*ω*)) only has non-zero elements for active-active block (i.e., non-zero matrix elements $\Sigma pqB\u2212CCSD(\omega )$ corresponding to active *p* and *q* spinorbital indices). This fact can be generalized to multiple sub-blocks and the coupling elements between these sub-blocks. In this case, a proper permutation of spinorbital basis transforms the multi-block structure into single block form (see Fig. 4). Although the same *T* and Λ operators are employed to construct Σ^{CCSD}(*ω*) and Σ^{B−CCSD}(*ω*) matrices, the active-active block of Σ^{CCSD}(*ω*) will be different from the Σ^{B−CCSD}(*ω*) one. Since the exact *X _{p}*(

*ω*) and

*Y*(

_{q}*ω*) operators are used in constructing the block approximation, the location of the poles in the B-CCSD approximation is the same as in the full GFCCSD approach (at least for the spatial symmetries of active orbitals).

To illustrate the performance of the B-CCSD approximation, we investigate the H_{2}O and N_{2} benchmark systems using the 3-21G basis set. For H_{2}O, we chose HOMO-1, HOMO, LUMO, LUMO+1 orbitals as active ones ((2,2) model space), whereas the active orbitals for N_{2} are defined by HOMO-2, HOMO-1, HOMO, LUMO, LUMO+1, LUMO+2 Hartree-Fock orbitals ((3,3) model space). In both cases, the active valence orbitals are energetically well separated from remaining (inactive) occupied and unoccupied orbitals. The eigenvalues of **F** + **Σ**^{CCSD}(*ω*) (*σ*_{i,CCSD}) and **F** + **Σ**^{B−CCSD}(*ω*) (*σ*_{i,B−CCSD}) matrices are shown in Tables VIII and IX. It is interesting to observe that for both systems and for all values of *ω* parameter considered here, the differences between *σ*_{i,CCSD} and *σ*_{i,B−CCSD} corresponding to active orbitals are negligibly small. For these two benchmark systems, this means the decoupling of various energy regimes can be achieved by a proper choice of the active orbitals defining the block approximation. Our block approximation approach can be viewed as an embedding method where the blocks of the CCSD Green function are embedded into the Hartree-Fock Green function. More advanced embedding schemes may utilize many-body perturbation theory to provide a “correlated background” for the GFCCSD block embedding. The related ideas will be discussed in a future paper. The utilization of larger cc-pVDZ basis set^{98} (see Tables X and XI) leads to similar conclusions as obtained for 3-21G basis set.

MO . | ϵ_{HF}
. | σ_{CC}(0)
. | σ_{B−CC}(0)
. | σ_{CC}(−0.1)
. | σ_{B−CC}(−0.1)
. | σ_{CC}(−0.15)
. | σ_{B−CC}(−0.15)
. | σ_{CC}(−0.2)
. | σ_{B−CC}(−0.2)
. |
---|---|---|---|---|---|---|---|---|---|

1 | −20.43 | −20.34 | −20.43 | −20.24 | −20.43 | −20.19 | −20.43 | −20.16 | −20.43 |

2 | −1.33 | −1.29 | −1.33 | −1.19 | −1.33 | −1.14 | −1.33 | −1.09 | −1.33 |

3 | −0.69 | −0.68 | −0.69 | −0.58 | −0.69 | −0.53 | −0.69 | −0.48 | −0.69 |

4 | −0.54 | −0.50 | −0.50 | −0.40 | −0.40 | −0.35 | −0.35 | −0.30 | −0.30 |

5 | −0.48 | −0.43 | −0.43 | −0.32 | −0.32 | −0.27 | −0.27 | −0.22 | −0.22 |

6 | 0.26 | 0.25 | 0.25 | 0.15 | 0.15 | 0.10 | 0.10 | 0.05 | 0.05 |

7 | 0.36 | 0.35 | 0.35 | 0.25 | 0.25 | 0.20 | 0.20 | 0.15 | 0.15 |

8 | 1.19 | 1.17 | 1.19 | 1.07 | 1.19 | 1.02 | 1.19 | 0.97 | 1.19 |

9 | 1.31 | 1.28 | 1.31 | 1.18 | 1.31 | 1.13 | 1.31 | 1.08 | 1.31 |

10 | 1.78 | 1.75 | 1.78 | 1.65 | 1.78 | 1.60 | 1.78 | 1.55 | 1.78 |

11 | 1.87 | 1.82 | 1.87 | 1.72 | 1.87 | 1.67 | 1.87 | 1.62 | 1.87 |

12 | 2.02 | 1.96 | 2.02 | 1.85 | 2.02 | 1.80 | 2.02 | 1.75 | 2.02 |

13 | 3.11 | 3.07 | 3.11 | 2.97 | 3.11 | 2.92 | 3.11 | 2.87 | 3.11 |

MO . | ϵ_{HF}
. | σ_{CC}(0)
. | σ_{B−CC}(0)
. | σ_{CC}(−0.1)
. | σ_{B−CC}(−0.1)
. | σ_{CC}(−0.15)
. | σ_{B−CC}(−0.15)
. | σ_{CC}(−0.2)
. | σ_{B−CC}(−0.2)
. |
---|---|---|---|---|---|---|---|---|---|

1 | −20.43 | −20.34 | −20.43 | −20.24 | −20.43 | −20.19 | −20.43 | −20.16 | −20.43 |

2 | −1.33 | −1.29 | −1.33 | −1.19 | −1.33 | −1.14 | −1.33 | −1.09 | −1.33 |

3 | −0.69 | −0.68 | −0.69 | −0.58 | −0.69 | −0.53 | −0.69 | −0.48 | −0.69 |

4 | −0.54 | −0.50 | −0.50 | −0.40 | −0.40 | −0.35 | −0.35 | −0.30 | −0.30 |

5 | −0.48 | −0.43 | −0.43 | −0.32 | −0.32 | −0.27 | −0.27 | −0.22 | −0.22 |

6 | 0.26 | 0.25 | 0.25 | 0.15 | 0.15 | 0.10 | 0.10 | 0.05 | 0.05 |

7 | 0.36 | 0.35 | 0.35 | 0.25 | 0.25 | 0.20 | 0.20 | 0.15 | 0.15 |

8 | 1.19 | 1.17 | 1.19 | 1.07 | 1.19 | 1.02 | 1.19 | 0.97 | 1.19 |

9 | 1.31 | 1.28 | 1.31 | 1.18 | 1.31 | 1.13 | 1.31 | 1.08 | 1.31 |

10 | 1.78 | 1.75 | 1.78 | 1.65 | 1.78 | 1.60 | 1.78 | 1.55 | 1.78 |

11 | 1.87 | 1.82 | 1.87 | 1.72 | 1.87 | 1.67 | 1.87 | 1.62 | 1.87 |

12 | 2.02 | 1.96 | 2.02 | 1.85 | 2.02 | 1.80 | 2.02 | 1.75 | 2.02 |

13 | 3.11 | 3.07 | 3.11 | 2.97 | 3.11 | 2.92 | 3.11 | 2.87 | 3.11 |

MO . | ϵ_{HF}
. | σ_{CC}(0)
. | σ_{B−CC}(0)
. | σ_{CC}(−0.1)
. | σ_{B−CC}(−0.1)
. | σ_{CC}(−0.15)
. | σ_{B−CC}(−0.15)
. |
---|---|---|---|---|---|---|---|

1 | −15.61 | −15.52 | −15.61 | −15.42 | −15.61 | −15.37 | −15.61 |

2 | −15.61 | −15.52 | −15.61 | −15.42 | −15.61 | −15.37 | −15.61 |

3 | −1.51 | −1.46 | −1.51 | −1.36 | −1.51 | −1.31 | −1.51 |

4 | −0.76 | −0.70 | −0.76 | −0.59 | −0.76 | −0.54 | −0.76 |

5 | −0.61 | −0.64 | −0.64 | −0.54 | −0.54 | −0.49 | −0.49 |

6 | −0.61 | −0.64 | −0.64 | −0.54 | −0.54 | −0.49 | −0.49 |

7 | −0.61 | −0.56 | −0.56 | −0.46 | −0.46 | −0.41 | −0.41 |

8 | 0.18 | 0.19 | 0.19 | 0.09 | 0.09 | 0.04 | 0.04 |

9 | 0.18 | 0.19 | 0.19 | 0.09 | 0.09 | 0.04 | 0.04 |

10 | 0.75 | 0.74 | 0.74 | 0.64 | 0.64 | 0.59 | 0.59 |

11 | 1.21 | 1.18 | 1.21 | 1.08 | 1.21 | 1.03 | 1.21 |

12 | 1.28 | 1.23 | 1.28 | 1.13 | 1.28 | 1.07 | 1.28 |

13 | 1.28 | 1.23 | 1.28 | 1.13 | 1.28 | 1.07 | 1.28 |

14 | 1.42 | 1.37 | 1.42 | 1.27 | 1.42 | 1.22 | 1.42 |

15 | 1.42 | 1.37 | 1.42 | 1.27 | 1.42 | 1.22 | 1.42 |

16 | 1.47 | 1.44 | 1.47 | 1.34 | 1.47 | 1.28 | 1.47 |

17 | 1.85 | 1.82 | 1.85 | 1.72 | 1.85 | 1.67 | 1.85 |

18 | 2.56 | 2.53 | 2.56 | 2.43 | 2.56 | 2.38 | 2.56 |

MO . | ϵ_{HF}
. | σ_{CC}(0)
. | σ_{B−CC}(0)
. | σ_{CC}(−0.1)
. | σ_{B−CC}(−0.1)
. | σ_{CC}(−0.15)
. | σ_{B−CC}(−0.15)
. |
---|---|---|---|---|---|---|---|

1 | −15.61 | −15.52 | −15.61 | −15.42 | −15.61 | −15.37 | −15.61 |

2 | −15.61 | −15.52 | −15.61 | −15.42 | −15.61 | −15.37 | −15.61 |

3 | −1.51 | −1.46 | −1.51 | −1.36 | −1.51 | −1.31 | −1.51 |

4 | −0.76 | −0.70 | −0.76 | −0.59 | −0.76 | −0.54 | −0.76 |

5 | −0.61 | −0.64 | −0.64 | −0.54 | −0.54 | −0.49 | −0.49 |

6 | −0.61 | −0.64 | −0.64 | −0.54 | −0.54 | −0.49 | −0.49 |

7 | −0.61 | −0.56 | −0.56 | −0.46 | −0.46 | −0.41 | −0.41 |

8 | 0.18 | 0.19 | 0.19 | 0.09 | 0.09 | 0.04 | 0.04 |

9 | 0.18 | 0.19 | 0.19 | 0.09 | 0.09 | 0.04 | 0.04 |

10 | 0.75 | 0.74 | 0.74 | 0.64 | 0.64 | 0.59 | 0.59 |

11 | 1.21 | 1.18 | 1.21 | 1.08 | 1.21 | 1.03 | 1.21 |

12 | 1.28 | 1.23 | 1.28 | 1.13 | 1.28 | 1.07 | 1.28 |

13 | 1.28 | 1.23 | 1.28 | 1.13 | 1.28 | 1.07 | 1.28 |

14 | 1.42 | 1.37 | 1.42 | 1.27 | 1.42 | 1.22 | 1.42 |

15 | 1.42 | 1.37 | 1.42 | 1.27 | 1.42 | 1.22 | 1.42 |

16 | 1.47 | 1.44 | 1.47 | 1.34 | 1.47 | 1.28 | 1.47 |

17 | 1.85 | 1.82 | 1.85 | 1.72 | 1.85 | 1.67 | 1.85 |

18 | 2.56 | 2.53 | 2.56 | 2.43 | 2.56 | 2.38 | 2.56 |

MO . | ϵ_{HF}
. | σ_{CC}(0)
. | σ_{B−CC}(0)
. | σ_{CC}(−0.1)
. | σ_{B−CC}(−0.1)
. | σ_{CC}(−0.15)
. | σ_{B−CC}(−0.15)
. |
---|---|---|---|---|---|---|---|

4 | −0.57 | −0.54 | −0.54 | −0.44 | −0.44 | −0.39 | −0.39 |

5 | −0.49 | −0.46 | −0.46 | −0.35 | −0.35 | −0.30 | −0.30 |

6 | 0.19 | 0.17 | 0.17 | 0.07 | 0.07 | 0.02 | 0.02 |

7 | 0.26 | 0.24 | 0.24 | 0.14 | 0.14 | 0.09 | 0.09 |

MO . | ϵ_{HF}
. | σ_{CC}(0)
. | σ_{B−CC}(0)
. | σ_{CC}(−0.1)
. | σ_{B−CC}(−0.1)
. | σ_{CC}(−0.15)
. | σ_{B−CC}(−0.15)
. |
---|---|---|---|---|---|---|---|

4 | −0.57 | −0.54 | −0.54 | −0.44 | −0.44 | −0.39 | −0.39 |

5 | −0.49 | −0.46 | −0.46 | −0.35 | −0.35 | −0.30 | −0.30 |

6 | 0.19 | 0.17 | 0.17 | 0.07 | 0.07 | 0.02 | 0.02 |

7 | 0.26 | 0.24 | 0.24 | 0.14 | 0.14 | 0.09 | 0.09 |

MO . | ϵ_{HF}
. | σ_{CC}(0)
. | σ_{B−CC}(0)
. | σ_{CC}(−0.1)
. | σ_{B−CC}(−0.1)
. | σ_{CC}(−0.15)
. | σ_{B−CC}(−0.15)
. |
---|---|---|---|---|---|---|---|

5 | −0.63 | −0.64 | −0.64 | −0.54 | −0.54 | −0.48 | −0.48 |

6 | −0.61 | −0.64 | −0.64 | −0.54 | −0.54 | −0.48 | −0.48 |

7 | −0.61 | −0.59 | −0.59 | −0.48 | −0.48 | −0.43 | −0.43 |

8 | 0.18 | 0.16 | 0.16 | 0.05 | 0.05 | 0.00 | 0.00 |

9 | 0.18 | 0.16 | 0.16 | 0.05 | 0.05 | 0.00 | 0.00 |

10 | 0.59 | 0.57 | 0.57 | 0.47 | 0.47 | 0.42 | 0.42 |

MO . | ϵ_{HF}
. | σ_{CC}(0)
. | σ_{B−CC}(0)
. | σ_{CC}(−0.1)
. | σ_{B−CC}(−0.1)
. | σ_{CC}(−0.15)
. | σ_{B−CC}(−0.15)
. |
---|---|---|---|---|---|---|---|

5 | −0.63 | −0.64 | −0.64 | −0.54 | −0.54 | −0.48 | −0.48 |

6 | −0.61 | −0.64 | −0.64 | −0.54 | −0.54 | −0.48 | −0.48 |

7 | −0.61 | −0.59 | −0.59 | −0.48 | −0.48 | −0.43 | −0.43 |

8 | 0.18 | 0.16 | 0.16 | 0.05 | 0.05 | 0.00 | 0.00 |

9 | 0.18 | 0.16 | 0.16 | 0.05 | 0.05 | 0.00 | 0.00 |

10 | 0.59 | 0.57 | 0.57 | 0.47 | 0.47 | 0.42 | 0.42 |

## VI. CONCLUSIONS

In this paper, we discussed key aspects of the CCSD Green function formulation including equations for calculating the frequency-dependent self-energy. We demonstrated that the main computational steps are defined by solving the CCSD and Λ-CCSD equations along with the *N* IP- and EA-type linear equations for the intermediate *X _{p}*(

*ω*) and

*Y*(

_{p}*ω*) operators. We demonstrated that our approach analytically includes all poles corresponding to the eigenvalues of the CCSD similarity transformed Hamiltonian in the 1

*h*/2

*h*− 1

*p*(IP) and 1

*p*/1

*h*− 2

*p*(EA) subspaces. This is an important advantage over the limited sum-over-states representation of the Green function. Consequently, our approach can be universally applied to calculate

*ω*-dependent CCSD self-energies for an arbitrary energy window. This procedure can be easily extended to GFCC formulations involving higher-rank excitations (for example, GFCCSDT). For the DNA/RNA bases, we obtain highly accurate IPs as compared to experiment. On the examples of 1,4-naphthtalenedion, benzoquinone, and anthracene, we showed that the location of the GFCCSD poles corresponding to higher-order ionization potentials is in a good agreement with the experimental data. The accuracies of the GFCCSD IP poles are competitive to the estimates obtained with popular 2p-h TDA, OVGA, D3, P3, NR2, and ADC(3) approximations.

We also introduced two new theoretical approximations: (1) the extension of the GFCCSD method enabling one to achieve arbitrary level of accuracy in the description of *N* ± 1 particle electronic states (see, for example, the GFCCSD-iT formalism) while preserving the CCSD description of the *N*-particle ground-state function and (2) block approximation of the GFCCSD matrix. We demonstrated that the novelty of the block approximation lies in the fact that the B-GFCCSD representation is not a perturbative approach for approximating the self-energy but rather is a direct formulation of the self-energy, and as such, it naturally preserves the pole structure of the full GFCCSD formalism (at least for the symmetries of active orbitals). This feature of the B-GFCCSD method is also valid in the case of the GFCCSD-iT type approaches. Also, since the diagonal GFCCSD approximation requires all *X _{p}*(

*ω*) and

*Y*(

_{q}*ω*) intermediates, it is as expensive as the full GFCCSD formulation. Due to recent progress in the development of parallel CC implementations that has resulted in significant reduction in time-to-solution, we believe this will enable widespread utilization of non-perturbative GFCC formulations.

## Acknowledgments

This work has been supported by the Extreme Scale Computing Initiative (K.K.), a Laboratory Directed Research and Development Program at Pacific Northwest National Laboratory. Computational support from Environmental Molecular Sciences Laboratory located at the Pacific Northwest National Laboratory that is operated for the US Department of Energy by the Battelle Memorial Institute under Contract No. DE-AC06.76RLO-1830. The computational work conducted by W.A.S. is supported by the U.S. Department of Energy under EPSCoR Grant No. DE-SC0012432 with additional support from the Louisiana Board of Regents.