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.

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 are1–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), vxc, 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 ago51–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) formulations45–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 Snijders51–53 The GFCC formalism hinges upon the CC bi-variational formalism62–64 utilizing different parametrizations of the bra (Ψ0(N)|) and ket (|Ψ0(N)) ground-state wavefunctions of a N-electron system,

Ψ0(N)|=Φ|(1+Λ)eT,
(1)
|Ψ0(N)=eT|Φ.
(2)

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

QeTHeT|Φ=0,
(3)
E0(N)=Φ|eTHeT|Φ,
(4)
Φ|(1+Λ)eTHeTQ=E0Φ|(1+Λ)Q,
(5)

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 (Tn and Λn),

T=n=1NeTn,
(6)
Λ=n=1NeΛn,
(7)

where Ne stands for the total number of correlated electrons. The Tn and Λn operators can be given by the following expressions:

Tn=1(n!)2i1,,in;a1,,anta1ani1inaa1aanainai1,
(8)
Λn=1(n!)2i1,,in;a1,,anλi1ina1anai1ainaanaa1,
(9)

where ta1ani1in and λi1ina1an are the amplitudes determining T and Λ operators. The indices i, j, k, … (i1, i2, …) and a, b, c, … (a1, a2, …) correspond to occupied and unoccupied spinorbitals, respectively. The ap (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

Gpq(ω)=Φ|(1+Λ)eTaq(ω+(HE0)iη)1apeT|Φ+Φ|(1+Λ)eTap(ω(HE0)+iη)1aqeT|Φ.
(10)

By introducing the resolution of identity eTeT = 1, the Green function can be rewritten as

Gpq(ω)=Φ|(1+Λ)āq(ω+(H̄E0)iη)1āp|Φ+Φ|(1+Λ)āp(ω(H̄E0)+iη)1āq|Φ,
(11)

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

āp=eTapeT,
(12)
āp=eTapeT,
(13)
H̄=eTHeT.
(14)

Using the Campbell-Hausdorff formula

eBAeB=A+[A,B]+12[[A,B],B]+,
(15)

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

āp=ap+[ap,T],
(16)
āp=ap+[ap,T].
(17)

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 Gpq(ω) matrix elements in a numerically efficient way, two sets of intermediate operators Xp and Yq defined on N − 1 and N + 1 particle Hilbert spaces are introduced as follows:

(ω+(H̄E0)iη)Xp(ω)|Φ=āp|Φ,
(18)
(ω(H̄E0)+iη)Yq(ω)|Φ=āq|Φ.
(19)

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

Gpq(ω)=Φ|(1+Λ)āqXp(ω)|Φ+Φ|(1+Λ)āpYq(ω)|Φ,
(20)

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:

Gpq(ω)=μ=1MN+1Ψ0(N)|ap|Ψμ(N+1)Ψμ(N+1)|aq|Ψ0(N)ω(Eμ(N+1)E0(N))+iη+ν=1MN1Ψ0(N)|aq|Ψν(N1)Ψν(N1)|ap|Ψ0(N)ω(E0(N)Eν(N1))iη,
(21)

where Eχ(N±1) and |Ψχ(N±1) denote energies and states of N ± 1 electron systems, and MN+1 and MN−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)),

Q(N+1)H̄Q(N+1)=μ=1MN+1|Rμ(N+1)Eμ(N+1)Lμ(N+1)|,
(22)
Q(N1)H̄Q(N1)=ν=1MN1|Rν(N1)Eν(N1)Rν(N1)|.
(23)

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,

Xp(ω)|Φ=(ω+(H̄E0)iη)1āp|Φ,
(24)
Yq(ω)|Φ=(ω(H̄E0)+iη)1āq|Φ,
(25)

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

Gpq(ω)=μ=1MN+1Φ|(1+Λ)āp|Rμ(N+1)Lμ(N+1)|āq|Φω(Eμ(N+1)E0(N))+iη+ν=1MN1Φ|(1+Λ)āq|Rν(N1)Lν(N1)|āp|Φω(E0(N)Eν(N1))iη.
(26)

It should be stressed that by solving the set of linear equations (18) and (19), the Xp(ω) and Yq(ω) operators account for all poles. The above expression becomes identical with the spectral form of the Green function (21), where Ψ0(N)|, |Ψ0(N), Ψμ(N1)|, |Ψμ(N1), Ψν(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.,

Ψ0(N)|=Φ|(1+Λ)eT,|Ψ0(N)=eT|Φ,
(27)
Ψμ(N1)|=Φ|Lμ(N1)eT,|Ψμ(N1)=Rμ(N1)eT|Φ,
(28)
Ψν(N+1)|=Φ|Lν(N+1)eT,|Ψν(N+1)=Rν(N+1)eT|Φ,
(29)

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

Φ|Lμ(N1)eT=Lμ(N1)|,
(30)
Rμ(N1)eT|Φ=|Rμ(N1).
(31)

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 Xp(ω) and Yq(ω) 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, Λ, Xp(ω), and Yq(ω) 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 (nNe) 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, Λ, Xp(ω), Yq(ω) 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, Λ, Xp(ω), and Yq(ω) operators:

Ti,ataiaaai+i<j,a<btabijaaabajai,
(32)
Λi,aλiaaiaa+i<j,a<bλijabaiajabaa,
(33)
Xp(ω)ixi(ω)pai+i<j,axaij(ω)paaajai,
(34)
Yq(ω)aya(ω)qaa+i,a<byabi(ω)qaaabai,
(35)

where the tensors tai, tabij, λia, λijab, xi(ω)p, xaij(ω)p, ya(ω)q, and yabi(ω)q are the amplitudes determining the T, Λ, Xp(ω), and Yq(ω) operators. The algebraic forms of Xp(ω) and Yq(ω) operators are analogous to the form of the RK 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,

T=T1+T2,
(36)
Λ=Λ1+Λ2,
(37)
(Q1(N)+Q2(N))eT1T2HeT1+T2|Φ=0,
(38)
E0(N)=Φ|eT1T2HeT1+T2|Φ,
(39)
Φ|(1+Λ1+Λ2)eT1T2HeT1+T2(Q1(N)+Q2(N))=E0(N)Φ|(1+Λ1+Λ2)(Q1(N)+Q2(N)),
(40)

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

Q1(N)=i,a|ΦiaΦia|,
(41)
Q2(N)=i<j;a<b|ΦijabΦijab|.
(42)

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

|Φia=XaXi|Φ,
(43)
|Φijab=XaXbXjXi|Φ.
(44)

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 Xp(ω) operator is represented as

Xp(ω)=Xp,1(ω)+Xp,2(ω)=ixi(ω)pai+i<j,axaij(ω)paaajai,
(45)

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

(Q1(N1)+Q2(N1))(ω+(H̄E0(N))iη)Xp(ω)|Φ=(Q1(N1)+Q2(N1))āp|Φ,
(46)

where the projection operators Q1(N1) and Q2(N1) are defined as follows:

Q1(N1)=i|ΦiΦi|,
(47)
Q2(N1)=i<j;a|ΦijaΦija|.
(48)

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

|Φi=ai|Φ,
(49)
|Φija=aa+ajai|Φ.
(50)

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

āp=ap+[ap,T1]+[ap,T2],
(51)

which can be written in a more explicit form as

āp={(52)apforpO,(53)ap+itpiai+i<j;atpaijaaajaiforpV,

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 Yq(ω) operator is represented in a way analogous to the standard EA-EOMCCSD parameterization (see Eq. (35)),

Yq(ω)=Yq,1(ω)+Yq,2(ω)=aya(ω)qaa+i;a<byabi(ω)paaabai,
(54)

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

(Q1(N+1)+Q2(N+1))(ω(H̄E0(N))+iη)Yq(ω)|Φ=(Q1(N+1)+Q2(N+1))āq|Φ,
(55)

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

Q1(N+1)=a|ΦaΦa|,
(56)
Q2(N+1)=i;a<b|ΦiabΦiab|,
(57)

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

|Φa=aa|Φ,
(58)
|Φiab=aaabai|Φ.
(59)

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

āq=aq+[aq,T1]+[aq,T2],
(60)

or equivalently by

āq={(61)aqataqaai;a<btabqiaaabaiforqO,(62)aqforqV.

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 Xp(ω) and Yq(ω) operators are determined, the corresponding contributions to the CCSD Green function matrix can be calculated from the equation,

Gpq(ω)=Φ|(1+Λ1+Λ2)āq(Xp,1(ω)+Xp,2(ω))|Φ+Φ|(1+Λ)āp(Yq,1(ω)+Yq,2(ω))|Φ.
(63)

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:

Φ|(1+Λ1+Λ2)āi(Xp,1+Xp,2)|Φ=xi(ω)pm,eλmeteixm(ω)pm,n,e<fλmneftefinxm(ω)p+m,eλmexeim(ω)pm<n,e,fλmnefteixfmn(ω)p,
(64)
Φ|(1+Λ1+Λ2)āa(Xp,1+Xp,2)|Φ=mλmaxm(ω)p+m<n,eλmnaexemn(ω)p,
(65)

and

Φ|(1+Λ1+Λ2)āi(Yq,1+Yq,2)|Φ=eλieye(ω)q+m,e<fλmiefyefm(ω)q,
(66)
Φ|(1+Λ1+Λ2)āa(Yq,1+Yq,2)|Φ=ya(ω)qm,eλmetamye(ω)q+m<n,e,fλmneftaemnyf(ω)q+m,eλmeyaem(ω)q+m,n,e<fλmneftanyefm(ω)q,
(67)

where the Xp(ω)-dependent terms contribute to the pth row of the CCSD Green function matrix, and the Yq(ω)-dependent terms contribute to the qth row of the CCSD Green function matrix. The equations for Xp(ω)’s and Yq(ω)’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 Xp’s and Yq’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 Xp(ω) (p = 1, …, N) or Yq(ω) (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 T1, T2, Λ1, and Λ2 operators (the so-called world group of Refs. 72–75).

FIG. 1.

Schematic representation of the CCSD Green function algorithm (see text for details).

FIG. 1.

Schematic representation of the CCSD Green function algorithm (see text for details).

Close modal

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 Xp(ω) and Yq(ω) operators, which are determined from Eqs. (46) and (55), respectively. In that case, the rank of excitation in the Xp(ω) and Yq(ω) 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̄=eT1T2HeT1+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 Xp,i(ω) and Yq,i(ω) (i = 1, 2, 3) operators. Subsequently, only Xp,i(ω) and Yq,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,

Σ(ω)=G01(ω)G1(ω).
(68)

G0(ω) 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 G0. 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

(G0)pq(ω)=δpqnpωϵpiη+δpqn̄pωϵp+iη,
(69)

where the np (n̄p=1np) 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 functionalities76,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.,

G=Gαα00Gββ,Gαα=Gββ.
(70)

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 NWChem45,46 that allow for the calculation of the most expensive terms of Eqs. (46) and (55) characterized by the action of H̄ operator onto Xp(ω)|Φ〉 (IP-EOMCCSD) and Yq(ω)|Φ〉 (EA-EOMCCSD) vectors. We will symbolically designate these terms as AIPxp(ω) and AEAyq(ω) and rewrite Eqs. (46) and (55) in matrix form

(ω+(AIPE0(N))iη)xp(ω)=bpIP,
(71)
(ω(AEAE0(N))+iη)yq(ω)=bqEA,
(72)

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 xp(ω) and yq(ω) vectors are represented as sums of their real and imaginary parts,

xp(ω)=xpRe(ω)+ixpIm(ω),
(73)
yq(ω)=yqRe(ω)+iyqIm(ω).
(74)

The equation for the real parts of the xp(ω) and yq(ω) 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

(AIP)2xpRe(ω)+2(ωE0(N))AIPxpRe(ω)+[(ωE0(N))2+η2]xpRe(ω)=(ωE0(N))bpIP+AIPbpIP,
(75)

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

(ω+(AIPE0(N)))xpIm(ω)=ηxpRe(ω).
(76)

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

(AEA)2yqRe(ω)2(ω+E0(N))AEAyqRe(ω)+[(ω+E0(N))2+η2]yqRe(ω)=(ω+E0(N))bqEAAEAbqEA,
(77)

and for the yqIm(ω) vectors,

(ω(AEAE0(N)))yqIm(ω)=ηyqRe(ω).
(78)

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 × N6 (per iteration). For given η and ω, the numerical cost of solving of all N IP/EA linear equations (75)(78) is 2N × N5 + 2N × N5 (per iteration), and the corresponding cost of forming the Green function matrix elements is proportional to 2N2 × N4. 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.

For illustrating some of the key features of the CCSD Green function, we have chosen some well-known benchmark systems including Ne, HF, H2O, and N2. 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 (H2O), and Ref. 80 (N2). 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 N2 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 sets98 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.

TABLE I.

The EOMCCSD ionization potentials for the uracil, thymine, cytosine, adenine, and guanine molecules calculated using various ab initio methods. For IP Expt. see Refs. 85–97.

MethodUracilThymineCytosineAdenineGuanine
G0W0(PBE)a 8.99 8.63 8.18 7.99 7.64 
G0W0(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/21b 9.22 8.87 8.54 8.18 7.91 
CASPT2(IPEA=0.25)/ANO-L 431/21b 9.42 9.07 8.73 8.37 8.09 
CCSD//CCSD/aug-cc-pVDZb 9.40 9.01 8.72 8.39 8.07 
CCSD(T)//CCSD/aug-cc-pVDZb 9.43 9.04 8.76 8.40 8.09 
IP-EOMCCSD//cc-pVTZc … 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] 
MethodUracilThymineCytosineAdenineGuanine
G0W0(PBE)a 8.99 8.63 8.18 7.99 7.64 
G0W0(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/21b 9.22 8.87 8.54 8.18 7.91 
CASPT2(IPEA=0.25)/ANO-L 431/21b 9.42 9.07 8.73 8.37 8.09 
CCSD//CCSD/aug-cc-pVDZb 9.40 9.01 8.72 8.39 8.07 
CCSD(T)//CCSD/aug-cc-pVDZb 9.43 9.04 8.76 8.40 8.09 
IP-EOMCCSD//cc-pVTZc … 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] 
a

Reference 84.

b

Reference 82.

c

Reference 83.

The G0W0 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-energy107 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 G0W0 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 G0W0(LDA) and G0W0(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 G0W0(LDA) and G0W0(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 G0W0(PBE) and GW(LDA) change slope with G0W0(PBE) having a positive slope while GW(LDA) has a negative one. The remaining higher order quantum chemistry methods and G0W0(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.

FIG. 2.

Comparison of ionization potentials obtained with various ab initio methods for uracil, thymine, cytosine, and adenine. Vertical dashed lines correspond to the experimental error bars.

FIG. 2.

Comparison of ionization potentials obtained with various ab initio methods for uracil, thymine, cytosine, and adenine. Vertical dashed lines correspond to the experimental error bars.

Close modal

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), G0W0 + SOSEX@PBE for 1,4-naphthalenedione (Table II), benzoquinone (Table III), and anthracene (Table IV).

TABLE II.

Comparison of higher order ionization potentials (in eV) obtained with the IP-EOMCCSD and other Green function/electron propagator methods for 1,4-naphthalenedione. In all IP-EOMCCSD calculations, core electrons were kept frozen.

Method2B22B12A2
D3a 10.90 9.89 9.85 
OVGFa 10.09 9.84 9.79 
NR2a 9.80 10.06 9.95 
2p-h TDAa 8.77 9.46 9.38 
ADC(3)a 10.37 9.93 9.90 
G0W0 + SOSEX@PBEa 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 
Method2B22B12A2
D3a 10.90 9.89 9.85 
OVGFa 10.09 9.84 9.79 
NR2a 9.80 10.06 9.95 
2p-h TDAa 8.77 9.46 9.38 
ADC(3)a 10.37 9.93 9.90 
G0W0 + SOSEX@PBEa 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 
a

Results taken from Ref. 110.

b

Experimental values of Ref. 113.

TABLE III.

Comparison of higher order ionization potentials (in eV) obtained with IP-EOMCCSD and other Green function/electron propagator methods for benzoquinone. In all IP-EOMCCSD calculations, core electrons were kept frozen.

Method2B2u2B3g2B3u2B1g
D3a 11.77 11.13 11.24 11.07 
OVGFa 10.90 10.43 11.00 11.07 
NR2a 10.35 10.13 11.17 11.34 
2p-h TDAa 9.88 9.17 10.66 10.86 
ADC(3)a 11.02 10.67 11.12 11.20 
G0W0 + SOSEX@PBEa 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 
Method2B2u2B3g2B3u2B1g
D3a 11.77 11.13 11.24 11.07 
OVGFa 10.90 10.43 11.00 11.07 
NR2a 10.35 10.13 11.17 11.34 
2p-h TDAa 9.88 9.17 10.66 10.86 
ADC(3)a 11.02 10.67 11.12 11.20 
G0W0 + SOSEX@PBEa 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 
a

Results taken from Ref. 110.

b

Experimental values of Ref. 114.

TABLE IV.

Comparison of higher order ionization potentials (in eV) obtained with IP-EOMCCSD and other Green function/electron propagator methods for anthracene. In all IP-EOMCCSD calculations, core electrons were kept frozen.

MethodB2gB3gAu
OVGFa 7.17 8.27 9.07 
P3a 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 
MethodB2gB3gAu
OVGFa 7.17 8.27 9.07 
P3a 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 
a

Results taken from Ref. 111.

b

Results taken from Ref. 46.

c

Experimental values of Ref. 112.

In Table II, we report on the tree lowest ionization potentials for the 2B2, 2B1, and 2A2 symmetries of the 1-4 naphthalenedione. The results are shown in Table II. While the IP-EOMCCSD aug-cc-pVTZ results for the 2B1 and 2A2 symmetries lie very close to the experimental values, for the 2B2 state, the IP-EOMCCSD/aug-cc-pVTZ ionization potential is characterized by ≃0.3 eV error. Similar behavior for the 2B2 ionization potential is disclosed by the D3, OVGF, NR2, and ADC(3) approaches. For this particular state, the G0W0+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 (2B2u), 0.16 (2B3g), 0.15 (2B3u), 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 G0W0+SOSEX@PBE approach. The efficiency of the GFCCSD formalism is clearly seen especially for the lowest 2B2u 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 B2g, B3g, and Au 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 (G0,ii)

Δoff−D=ij(Gij)2,
(79)
Δii=(GiiG0,ii).
(80)

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 N2 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.

TABLE V.

The largest off-diagonal element of the CCSD Green function matrix for the Ne, HF, H2O, and N2 systems. The second column shows the off-diagonality of the Green function matrix G(ω = 0; η = 0) defined as ijGij2. In all calculations, 3-21G basis set was used and all electrons were correlated.

SystemmaxGijΔoff−D=ijGij2
Ne 0.005 0.010 
HF 0.021 0.055 
H20.036 0.090 
N2 0.104 0.177 
SystemmaxGijΔoff−D=ijGij2
Ne 0.005 0.010 
HF 0.021 0.055 
H20.036 0.090 
N2 0.104 0.177 

In Tables VI and VII we show the differences between diagonal elements CCSD (G(ω = 0)) and the Hartree-Fock (G0(ω = 0)) Green function matrices (Δii values), Hartree-Fock orbital energies (ϵHF), and the eigenvalues σCCSD of F + Σ(ω = 0) matrix (where F is the Fock matrix) and differences between corresponding ϵHF’s and σCCSD’s for the Ne and HF systems. Again, the Δ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 N2, 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.

TABLE VI.

The differences between diagonal elements of CCSD (G(ω = 0)) and Hartree-Fock (G0(ω = 0)) Green function matrices, Hartree-Fock orbital energies (ϵHF), eigenvalues σCCSD of F + Σ(ω = 0) matrix (F and Σ(ω = 0) stand for the Fock matrix and CCSD self-energy calculated for ω = 0), and differences between corresponding ϵHF’s and σCCSD’s for the Ne and HF systems. All calculations were performed using 3-21G basis set, and all electrons were correlated in the calculations. Shown here in brackets are the occupation numbers of the orbitals.

Ne:Ne:Ne:Ne:HF:HF:HF:HF:
Δii = GiiG0;iiϵHF(occ.)σCCSDϵHFσCCSDΔii = GiiG0;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 = GiiG0;iiϵHF(occ.)σCCSDϵHFσCCSDΔii = GiiG0;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 
TABLE VII.

The differences between diagonal elements of CCSD (G(ω = 0)) and Hartree-Fock (G0(ω = 0)) Green function matrices, Hartree-Fock orbital energies (ϵHF), eigenvalues σCCSD of F + Σ(ω = 0) matrix (F and Σ(ω = 0) stand for the Fock matrix and CCSD self-energy calculated for ω = 0), and differences between corresponding ϵHF’s and σCCSD’s for the H2O and N2 systems. All calculations were performed using 3-21G basis set, and all electrons were correlated in the calculations. Shown here in brackets are the occupation numbers of the orbitals.

H2O:H2O:H2O:H2O:N2:N2:N2:N2:
Δii = GiiG0;iiϵHF(occ.)σCCSDϵHFσCCSDΔii = GiiG0;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 
H2O:H2O:H2O:H2O:N2:N2:N2:N2:
Δii = GiiG0;iiϵHF(occ.)σCCSDϵHFσCCSDΔii = GiiG0;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 
FIG. 3.

Graphical representation of the differences between diagonal elements of the CCSD Green function matrix (Gii) and their Hartree-Fock counterparts (G0,ii) for ω = 0 and for Ne, HF, H2O, and N2 benchmark systems in 3-21G basis set.117 All electrons were correlated in the GFCCSD calculations.

FIG. 3.

Graphical representation of the differences between diagonal elements of the CCSD Green function matrix (Gii) and their Hartree-Fock counterparts (G0,ii) for ω = 0 and for Ne, HF, H2O, and N2 benchmark systems in 3-21G basis set.117 All electrons were correlated in the GFCCSD calculations.

Close modal

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 2p 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 H2O 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 (noa, nva) stands for the active space composed of noa (nov) highest lying occupied (lowest lying unoccupied) orbitals. In the case of the H2O molecule, the largest differences correspond to molecular orbitals contributing to the (2,2) active space. A similar pattern is observed for the N2 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.,

GpqB−CCSD(ω)=δp(I)δpqnpωϵpiη+δp(I)δpqn̄pωϵp+iη+GpqCCSD(ω)δp(A)δq(A),
(81)

where

δp(I)={(82)1for inactivepspinorbital,(83)0for activepspinorbital,

and δp(A)=1δ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 Xp and Yq operators required to form B-GFCCSD (note: only those Xp and Yq operators defined by the active 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 2h − 1p and 1h − 2p 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 GB−CCSD(ω) matrix, the corresponding self-energy matrix (ΣB−CCSD(ω)) only has non-zero elements for active-active block (i.e., non-zero matrix elements ΣpqB−CCSD(ω) 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 Xp(ω) and Yq(ω) 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).

FIG. 4.

Graphical representation of the block approximation of the CCSD Green function (see text for details).

FIG. 4.

Graphical representation of the block approximation of the CCSD Green function (see text for details).

Close modal

To illustrate the performance of the B-CCSD approximation, we investigate the H2O and N2 benchmark systems using the 3-21G basis set. For H2O, we chose HOMO-1, HOMO, LUMO, LUMO+1 orbitals as active ones ((2,2) model space), whereas the active orbitals for N2 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 set98 (see Tables X and XI) leads to similar conclusions as obtained for 3-21G basis set.

TABLE VIII.

The ω-dependence (ω = 0.0, −0.1, −0.15, −0.2) of σCC(ω) and σB−CC(ω) eigenvalues of the F + ΣCCSD(ω) and F + ΣB−CCSD(ω) matrices, respectively, for the H2O system in 3-21G basis set (all electrons were correlated). The active orbitals used to define block approximation of the CCSD Green function correspond to the 4th, 5th, 6th, and 7th Hartree-Fock molecular orbitals. Here, ϵHF ’s denote Hartree-Fock orbital energies. The boldface values are the eigenvalues corresponding to active orbitals.

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)
−20.43 −20.34 −20.43 −20.24 −20.43 −20.19 −20.43 −20.16 −20.43 
−1.33 −1.29 −1.33 −1.19 −1.33 −1.14 −1.33 −1.09 −1.33 
−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 
1.19 1.17 1.19 1.07 1.19 1.02 1.19 0.97 1.19 
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)
−20.43 −20.34 −20.43 −20.24 −20.43 −20.19 −20.43 −20.16 −20.43 
−1.33 −1.29 −1.33 −1.19 −1.33 −1.14 −1.33 −1.09 −1.33 
−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 
1.19 1.17 1.19 1.07 1.19 1.02 1.19 0.97 1.19 
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 
TABLE IX.

The ω-dependence (ω = 0.0, − 0.1, − 0.15) of σCC(ω) and σB−CC(ω) eigenvalues of the F + ΣCCSD(ω) and F + ΣB−CCSD(ω) matrices, respectively, for the N2 system in 3-21G basis set (all electrons were correlated). The active orbitals used to define block approximation of the CCSD Green function correspond to the 5th, 6th, 7th, 8th, 9th, and 10th Hartree-Fock molecular orbitals. Here, ϵHF’s denote Hartree-Fock orbital energies. The boldface values are the eigenvalues corresponding to active orbitals.

MOϵHFσCC(0)σB−CC(0)σCC(−0.1)σB−CC(−0.1)σCC(−0.15)σB−CC(−0.15)
−15.61 −15.52 −15.61 −15.42 −15.61 −15.37 −15.61 
−15.61 −15.52 −15.61 −15.42 −15.61 −15.37 −15.61 
−1.51 −1.46 −1.51 −1.36 −1.51 −1.31 −1.51 
−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)
−15.61 −15.52 −15.61 −15.42 −15.61 −15.37 −15.61 
−15.61 −15.52 −15.61 −15.42 −15.61 −15.37 −15.61 
−1.51 −1.46 −1.51 −1.36 −1.51 −1.31 −1.51 
−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 
TABLE X.

The ω-dependence (ω = 0.0, − 0.1, − 0.15) of σCC(ω) and σB−CC(ω) eigenvalues of the F + ΣCCSD(ω) and F + ΣB−CCSD(ω) matrices, respectively, for the H2O system in cc-pVDZ basis set (all electrons were correlated). The active orbitals used to define block approximation of the CCSD Green function that correspond to the 4th, 5th, 6th, and 7th Hartree-Fock molecular orbitals are shown here. The only eigenvalues shown in the table are those corresponding to the active orbitals. Here, ϵHF ’s denote Hartree-Fock orbital energies.

MOϵHFσCC(0)σB−CC(0)σCC(−0.1)σB−CC(−0.1)σCC(−0.15)σB−CC(−0.15)
−0.57 −0.54 −0.54 −0.44 −0.44 −0.39 −0.39 
−0.49 −0.46 −0.46 −0.35 −0.35 −0.30 −0.30 
0.19 0.17 0.17 0.07 0.07 0.02 0.02 
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)
−0.57 −0.54 −0.54 −0.44 −0.44 −0.39 −0.39 
−0.49 −0.46 −0.46 −0.35 −0.35 −0.30 −0.30 
0.19 0.17 0.17 0.07 0.07 0.02 0.02 
0.26 0.24 0.24 0.14 0.14 0.09 0.09 
TABLE XI.

The ω-dependence (ω = 0.0, − 0.1, − 0.15) of σCC(ω) and σB−CC(ω) eigenvalues of the F + ΣCCSD(ω) and F + ΣB−CCSD(ω) matrices, respectively, for the N2 system in cc-pVDZ basis set (all electrons were correlated). The active orbitals used to define block approximation of the CCSD Green function that correspond to the 5th, 6th, 7th, 8th, 9th, and 10th Hartree-Fock molecular orbitals are shown here. The only eigenvalues shown in the table are those corresponding to the active orbitals. Here, ϵHF’s denote Hartree-Fock orbital energies.

MOϵHFσCC(0)σB−CC(0)σCC(−0.1)σB−CC(−0.1)σCC(−0.15)σB−CC(−0.15)
−0.63 −0.64 −0.64 −0.54 −0.54 −0.48 −0.48 
−0.61 −0.64 −0.64 −0.54 −0.54 −0.48 −0.48 
−0.61 −0.59 −0.59 −0.48 −0.48 −0.43 −0.43 
0.18 0.16 0.16 0.05 0.05 0.00 0.00 
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)
−0.63 −0.64 −0.64 −0.54 −0.54 −0.48 −0.48 
−0.61 −0.64 −0.64 −0.54 −0.54 −0.48 −0.48 
−0.61 −0.59 −0.59 −0.48 −0.48 −0.43 −0.43 
0.18 0.16 0.16 0.05 0.05 0.00 0.00 
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 

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 Xp(ω) and Yp(ω) operators. We demonstrated that our approach analytically includes all poles corresponding to the eigenvalues of the CCSD similarity transformed Hamiltonian in the 1h/2h − 1p (IP) and 1p/1h − 2p (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 Xp(ω) and Yq(ω) 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.

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.

1.
L.
Hedin
,
Phys. Rev.
139
,
A796
(
1965
).
2.
L. S.
Cederbaum
,
J. Phys. B: At., Mol. Opt. Phys.
8
,
290
(
1975
).
3.
J.
Schirmer
,
Phys. Rev. A
26
,
2395
(
1982
).
4.
W.
von Niessen
,
J.
Schirmer
, and
L. S.
Cederbaum
,
Comput. Phys. Rep.
1
,
57
(
1984
).
5.
H. G.
Weikert
,
H. D.
Meyer
,
L. S.
Cederbaum
, and
F.
Tarantelli
,
J. Chem. Phys.
104
,
7122
(
1996
).
6.
W.
von Niessen
,
G. H. F.
Diercksen
, and
L. S.
Cederbaum
,
J. Chem. Phys.
67
,
4124
(
1977
).
7.
E.
Andersen
and
J.
Simons
,
J. Chem. Phys.
66
,
2427
(
1977
).
8.
Y.
Öhrn
and
G.
Born
, in
Molecular Electron Propagator Theory and Calculations
, edited by
P.-O.
Löwdin
,
Advances in Chemical Physics
(
Academic Press
,
1981
), Vol.
13
, pp.
1
88
.
9.
M. F.
Herman
,
K. F.
Freed
, and
D. L.
Yeager
,
Advances in Chemical Physics
(
John Wiley & Sons, Inc.
,
2007
), pp.
1
69
.
10.
J.
Geertsen
,
J.
Oddershede
, and
G. E.
Scuseria
,
Int. J. Quantum Chem.
32
,
475
(
1987
).
11.
P.
Albertsen
,
P.
Jorgensen
, and
D. L.
Yeager
,
Mol. Phys.
41
,
409
(
1980
).
12.
J. V.
Ortiz
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
3
,
123
(
2013
).
13.
Applied Many-Body Methods in Spectroscopy and Electronic Structure
, edited by
D.
Mukherjee
(
Plenum Publishing Co.
,
New York and London
,
1992
).
14.
A.
Ankudinov
,
B.
Ravel
,
J.
Rehr
, and
S.
Conradson
,
Phys. Rev. B
58
,
7565
(
1998
).
15.
J.
Rehr
and
R.
Albers
,
Rev. Mod. Phys.
72
,
621
(
2000
).
16.
J. J.
Kas
,
J. J.
Rehr
, and
L.
Reining
,
Phys. Rev. B
90
,
085112
(
2014
).
17.
S.
Faleev
,
M.
van Schilfgaarde
, and
T.
Kotani
,
Phys. Rev. Lett.
93
,
126406
(
2004
).
18.
M.
van Schilfgaarde
,
T.
Kotani
, and
S.
Faleev
,
Phys. Rev. Lett.
96
,
226402
(
2006
).
19.
J. B.
Neaton
,
M. S.
Hybertsen
, and
S. G.
Louie
,
Phys. Rev. Lett.
97
,
216405
(
2006
).
20.
G.
Samsonidze
,
M.
Jain
,
J.
Deslippe
,
M. L.
Cohen
, and
S. G.
Louie
,
Phys. Rev. Lett.
107
,
186404
(
2011
).
21.
M. J.
van Setten
,
F.
Weigend
, and
F.
Evers
,
J. Chem. Theory Comput.
9
,
232
(
2013
).
22.
S.
Hirata
,
M. R.
Hermes
,
J.
Simons
, and
J. V.
Ortiz
,
J. Chem. Theory Comput.
11
,
1595
(
2015
).
23.
J. J.
Phillips
and
D.
Zgid
,
J. Chem. Phys.
140
,
241101
(
2014
).
24.
T. N.
Lan
,
A. A.
Kananenka
, and
D.
Zgid
,
J. Chem. Phys.
143
,
241102
(
2015
).
25.
A. J.
Heeger
,
Rev. Mod. Phys.
73
,
681
(
2001
).
26.
S.
Barlow
,
Q.
Zhang
,
B.
Kaafarani
,
C.
Risko
,
F.
Amy
,
C.
Chan
,
B.
Domercq
,
Z.
Starikova
,
M.
Antipin
,
T.
Timofeeva
,
B.
Kippelen
,
J.-L.
Brédas
,
A.
Kahn
, and
S.
Marder
,
Chem. - Eur. J.
13
,
3537
(
2007
).
27.
S.
Günes
,
H.
Neugebauer
, and
N. S.
Sariciftci
,
Chem. Rev.
107
,
1324
(
2007
).
28.
X.
Zhan
,
C.
Risko
,
F.
Amy
,
C.
Chan
,
W.
Zhao
,
S.
Barlow
,
A.
Kahn
,
J.-L.
Brédas
, and
S. R.
Marder
,
J. Am. Chem. Soc.
127
,
9021
(
2005
).
29.
M. D.
Perez
,
C.
Borek
,
S. R.
Forrest
, and
M. E.
Thompson
,
J. Am. Chem. Soc.
131
,
9281
(
2009
).
30.
J.-L.
Brédas
,
J. E.
Norton
,
J.
Cornil
, and
V.
Coropceanu
,
Acc. Chem. Res.
42
,
1691
(
2009
).
31.
J. L.
Delgado
,
E.
Espildora
,
M.
Liedtke
,
A.
Sperlich
,
D.
Rauh
,
A.
Baumann
,
C.
Deibel
,
V.
Dyakonov
, and
N.
Martin
,
Chem. - Eur. J.
15
,
13474
(
2009
).
32.
S.
Baroni
,
R.
Resta
,
A.
Baldereschi
, and
M.
Peressi
,
Spectroscopy of Semiconductor Microstructures
(
Springer
,
1989
), pp.
251
271
.
33.
W.
Mönch
,
Semiconductor Surfaces and Interfaces
(
Springer
,
2001
), Vol.
26
.
34.
S. M.
Sze
and
K. K.
Ng
,
Physics of Semiconductor Devices
(
John Wiley & Sons
,
2006
).
35.
S.
Okamoto
and
T. A.
Maier
,
Phys. Rev. Lett.
101
,
156401
(
2008
).
36.
H.
Luo
,
Z.
Yamani
,
Y.
Chen
,
X.
Lu
,
M.
Wang
,
S.
Li
,
T. A.
Maier
,
S.
Danilkin
,
D. T.
Adroja
, and
P.
Dai
,
Phys. Rev. B
86
,
024508
(
2012
).
37.
G.
Kotliar
,
S. Y.
Savrasov
,
K.
Haule
,
V. S.
Oudovenko
,
O.
Parcollet
, and
C. A.
Marianetti
,
Rev. Mod. Phys.
78
,
865
(
2006
).
38.
F.
Coester
,
Nucl. Phys.
7
,
421
(
1958
).
39.
F.
Coester
and
H.
Kümmel
,
Nucl. Phys.
17
,
477
(
1960
).
40.
J.
Čížek
,
J. Chem. Phys.
45
,
4256
(
1966
).
41.
J.
Paldus
,
I.
Shavitt
, and
J.
Čížek
,
Phys. Rev. A
5
,
50
(
1972
).
42.
G.
Purvis
and
R.
Bartlett
,
J. Chem. Phys.
76
,
1910
(
1982
).
43.
K.
Kowalski
,
K.
Bhaskaran-Nair
, and
W. A.
Shelton
,
J. Chem. Phys.
141
,
094102
(
2014
).
44.
K.
Kowalski
,
K.
Bhaskaran-Nair
, and
W. A.
Shelton
, “
Coupled Cluster Green Function and Accurate Self-Energies: A 1st Principles Approach for Electron Transport
” (unpublished).
45.
K.
Bhaskaran-Nair
,
K.
Kowalski
,
J.
Moreno
,
M.
Jarrell
, and
W. A.
Shelton
,
J. Chem. Phys.
141
,
074304
(
2014
).
46.
K.
Bhaskaran-Nair
,
K.
Kowalski
,
M.
Jarrell
,
J.
Moreno
, and
W. A.
Shelton
,
Chem. Phys. Lett.
641
,
146
(
2015
).
47.
K.
Bhaskaran-Nair
,
M.
Valiev
,
S. H. M.
Deng
,
W. A.
Shelton
,
K.
Kowalski
, and
X.-B.
Wang
,
J. Chem. Phys.
143
,
224301
(
2015
).
48.
J.
Noga
and
R.
Bartlett
,
J. Chem. Phys.
86
,
7041
(
1987
).
49.
J.
Noga
and
R.
Bartlett
,
J. Chem. Phys.
89
,
3041
(
1988
).
50.
G.
Scuseria
and
H.
Schaefer
,
Chem. Phys. Lett.
152
,
382
(
1988
).
51.
M.
Nooijen
and
J.
Snijders
,
Int. J. Quantum Chem.
44
,
55
(
1992
).
52.
M.
Nooijen
and
J.
Snijders
,
Int. J. Quantum Chem.
48
,
15
(
1993
).
53.
M.
Nooijen
and
J.
Snijders
,
J. Chem. Phys.
102
,
1681
(
1995
).
54.
L.
Meissner
and
R.
Bartlett
,
Int. J. Quantum Chem. Symp.
48
(
S27
),
67
(
1993
).
55.
J.
Stanton
and
J.
Gauss
,
J. Chem. Phys.
103
,
1064
(
1995
).
56.
M.
Nooijen
and
R. J.
Bartlett
,
J. Chem. Phys.
102
,
3629
(
1995
).
57.
M.
Musiał
,
J. Chem. Phys.
136
,
134111
(
2012
).
58.
M.
Musiał
,
K.
Kowalska-Szojda
,
D. I.
Lyakh
, and
R. J.
Bartlett
,
J. Chem. Phys.
138
,
194103
(
2013
).
59.
O.
Demel
,
K. R.
Shamasundar
,
L.
Kong
, and
M.
Nooijen
,
J. Phys. Chem. A
112
,
11895
(
2008
).
60.
T.
Kuś
and
A. I.
Krylov
,
J. Chem. Phys.
135
,
084109
(
2011
).
61.
J.
Shen
and
P.
Piecuch
,
J. Chem. Phys.
138
,
194102
(
2013
).
62.
J.
Arponen
,
Ann. Phys.
151
,
311
(
1983
).
63.
J. F.
Stanton
and
R. J.
Bartlett
,
J. Chem. Phys.
98
,
7029
(
1993
).
64.
J.
Stanton
and
R.
Bartlett
,
J. Chem. Phys.
99
,
5178
(
1993
).
65.
R. J.
Bartlett
and
M.
Musiał
,
Rev. Mod. Phys.
79
,
291
(
2007
).
66.
B.
Jeziorski
and
H.
Monkhorst
,
Phys. Rev. A
24
,
1668
(
1981
).
67.
U.
Mahapatra
,
B.
Datta
, and
D.
Mukherjee
, in
Coupled Cluster Theory Electron Correlation Workshop on 50 Years of the Correlation Problem, Cedar Key, FL, 15–19 June 1997
[
Mol. Phys.
94
,
157
(
1998
)].
68.
U.
Mahapatra
,
B.
Datta
, and
D.
Mukherjee
,
J. Chem. Phys.
110
,
6171
(
1999
).
69.
J.
Pittner
,
P.
Nachtigall
,
P.
Carsky
,
J.
Masik
, and
I.
Hubac
,
J. Chem. Phys.
110
,
10275
(
1999
).
70.
J.
Pittner
,
J. Chem. Phys.
118
,
10876
(
2003
).
71.
K.
Bhaskaran-Nair
and
K.
Kowalski
,
J. Chem. Phys.
138
,
204114
(
2013
).
72.
J.
Brabec
,
J.
Pittner
,
H. J. J.
van Dam
,
E.
Apra
, and
K.
Kowalski
,
J. Chem. Theory Comput.
8
,
487
(
2012
).
73.
K.
Bhaskaran-Nair
,
J.
Brabec
,
E.
Apra
,
H. J. J.
van Dam
,
J.
Pittner
, and
K.
Kowalski
,
J. Chem. Phys.
137
,
094112
(
2012
).
74.
J.
Brabec
,
K.
Bhaskaran-Nair
,
K.
Kowalski
,
J.
Pittner
, and
H. J. J.
van Dam
,
Chem. Phys. Lett.
542
,
128
(
2012
).
75.
K.
Bhaskaran-Nair
,
W.
Ma
,
S.
Krishnamoorthy
,
O.
Villa
,
H. J. J.
van Dam
,
E.
Apra
, and
K.
Kowalski
,
J. Chem. Theory Comput.
9
,
1949
(
2013
).
76.
S.
Hirata
,
J. Phys. Chem. A
107
,
9887
(
2003
).
77.
K.
Kowalski
,
S.
Krishnamoorthy
,
R.
Olson
,
V.
Tipparaju
, and
E.
Apra
, in
2011 International Conference for High Performance Computing, Networking, Storage and Analysis (SC)
(
IEEE
,
2011
), article No. 72, pp.
1
10
.
78.
M.
Valiev
,
E.
Bylaska
,
N.
Govind
,
K.
Kowalski
,
T.
Straatsma
,
H. V.
Dam
,
D.
Wang
,
J.
Nieplocha
,
E.
Apra
,
T.
Windus
, and
W.
de Jong
,
Comput. Phys. Commun.
181
,
1477
(
2010
).
79.
E. R.
Davidson
,
J. Comput. Phys.
17
,
87
(
1975
).
80.
B.
Rosen
,
Spectroscopic Data Relative to Diatomic Molecules, International Tables of Selected Constants
(
Pergamon Press
,
1970
).
81.
P.
Jensen
,
S.
Tashkun
, and
V.
Tyuterev
,
J. Mol. Spectrosc.
168
,
271
(
1994
).
82.
D.
Roca-Sanjuàn
,
M.
Rubio
,
M.
Merchàn
, and
L.
Serrano-Andrés
,
J. Chem. Phys.
125
,
084302
(
2006
).
83.
K. B.
Bravaya
,
O.
Kostko
,
S.
Dolgikh
,
A.
Landau
,
M.
Ahmed
, and
A. I.
Krylov
,
J. Phys. Chem. A
114
,
12305
(
2010
).
84.
X.
Qian
,
P.
Umari
, and
N.
Marzari
,
Phys. Rev. B
84
,
075103
(
2011
).
85.
N.
Hush
and
A. S.
Cheung
,
Chem. Phys. Lett.
34
,
11
(
1975
).
86.
G.
Lauer
,
W.
Schäfer
, and
A.
Schweig
,
Tetrahedron Lett.
16
,
3939
(
1975
).
87.
V.
Orlov
,
A.
Smirnov
, and
Y.
Varshavsky
,
Tetrahedron Lett.
17
,
4377
(
1976
).
88.
D.
Dougherty
,
E.
Younathan
,
R.
Voll
,
S.
Abdulnur
, and
S.
McGlynn
,
J. Electron Spectrosc. Relat. Phenom.
13
,
379
(
1978
).
89.
A.
Padva
,
T. O.
Donnell
, and
P.
Lebreton
,
Chem. Phys. Lett.
41
,
278
(
1976
).
90.
J.
Lin
,
C.
Yu
,
S.
Peng
,
I.
Akiyama
,
K.
Li
,
L. K.
Lee
, and
P. R.
LeBreton
,
J. Am. Chem. Soc.
102
,
4627
(
1980
).
91.
S.
Urano
,
X.
Yang
, and
P. R.
LeBreton
,
J. Mol. Struct.
214
,
315
(
1989
).
92.
D.
Dougherty
and
S. P.
McGlynn
,
J. Chem. Phys.
67
,
1289
(
1977
).
93.
J.
Lin
,
C.
Yu
,
S.
Peng
,
I.
Akiyama
,
K.
Li
,
L. K.
Lee
, and
P. R.
LeBreton
,
J. Phys. Chem.
84
,
1006
(
1980
).
94.
C.
Yu
,
S.
Peng
,
I.
Akiyama
,
J.
Lin
, and
P. R.
LeBreton
,
J. Am. Chem. Soc.
100
,
2303
(
1978
).
95.
D.
Dougherty
,
K.
Wittel
,
J.
Meeks
, and
S. P.
McGlynn
,
J. Am. Chem. Soc.
98
,
3815
(
1976
).
96.
M. E.
Jacox
,
J. Phys. Chem. Ref. Data
17
,
269
(
1988
).
97.
B. I.
Verkin
,
L. F.
Sukhodub
, and
I. K.
Yanson
,
Dokl. Akad. Nauk SSSR
228
,
1452
(
1976
).
98.
T. H.
Dunning
,
J. Chem. Phys.
90
,
1007
(
1989
).
99.
A. J.
Sadlej
,
Collect. Czech. Chem. Commun.
53
,
1995
(
1988
).
100.
V. G.
Zakrzewski
,
O.
Dolgounitcheva
,
A. V.
Zakjevskii
, and
J. V.
Orti
,
Int. J. Quantum Chem.
110
,
2918
(
2010
).
101.
J.
Baker
and
B. T.
Pickup
,
Chem. Phys. Lett.
76
,
537
(
1980
).
102.
J.
Schirmer
,
L. S.
Cederbaum
, and
O.
Walter
,
Phys. Rev. A
28
,
1237
(
1983
).
103.
A. B.
Trofimov
,
G.
Stetler
, and
J.
Schirmer
,
J. Chem. Phys.
111
,
9982
(
1999
).
104.
J. V.
Ortiz
, in
Computational Chemistry: Reviews of Current Trends
, edited by
J.
Leszczynski
(
World Scientific
,
Singapore
,
1997
), Vol.
2
, p.
1
.
105.
J. V.
Ortiz
,
J. Chem. Phys.
108
,
1008
(
1998
).
106.
J. V.
Ortiz
,
J. Chem. Phys.
104
,
7599
(
1996
).
107.
P.
Umari
,
G.
Stenuit
, and
S.
Baroni
,
Phys. Rev. B
79
,
201104
(
2009
).
108.
P.
Umari
,
G.
Stenuit
, and
S.
Baroni
,
Phys. Rev. B
81
,
115104
(
2010
).
109.
C.
Faber
,
C.
Attaccalite
,
V.
Olevano
,
E.
Runge
, and
X.
Blase
,
Phys. Rev. B
83
,
115123
(
2011
).
110.
O.
Dolgounitcheva
,
M.
Diaz-Tinoco
,
V. G.
Zakrzewski
,
R. M.
Richard
,
N.
Marom
,
C. D.
Sherrill
, and
J. V.
Ortiz
,
J. Chem. Theory Comput.
12
,
627
(
2016
).
111.
V. G.
Zakrzewski
,
O.
Dolgounitcheva
, and
J. V.
Ortiz
,
J. Chem. Phys.
105
,
8748
(
1996
).
112.
R.
Boschi
,
E.
Clar
, and
W.
Schmidt
,
J. Chem. Phys.
60
,
4406
(
1974
).
113.
S.
Millefiori
,
A.
Gulino
, and
M.
Casarin
,
J. Chim. Phys. Phys.-Chim. Biol.
87
,
317
(
1990
).
114.
C. R.
Brundle
,
H.
Basch
,
N. A.
Kuebler
, and
M. B.
Robin
,
Inorg. Chem.
11
,
20
(
1972
).
115.
H.
Müther
,
T.
Taigel
, and
T.
Kuo
,
Nucl. Phys. A
482
,
601
(
1988
).
116.
H.
Müther
and
L.
Skouras
,
Phys. Lett. B
306
,
201
(
1993
).
117.
J. S.
Binkley
,
J. A.
Pople
, and
W. J.
Hehre
,
J. Am. Chem. Soc.
102
,
939
(
1980
).