An Alternative Hierarchy of Electron Correlation beyond the Electron

An Alternative Hierarchy of Electron Correlation beyond the Electron Pair Approximation. Werner Kutzelnigg*. Lehrstuhl für Theoretische Chemie, ...
0 downloads 0 Views 395KB Size
J. Phys. Chem. A 2010, 114, 8913–8922

8913

An Alternative Hierarchy of Electron Correlation beyond the Electron Pair Approximation† Werner Kutzelnigg* Lehrstuhl fu¨r Theoretische Chemie, Ruhr-UniVersita¨t Bochum, D-44780 Bochum, Germany ReceiVed: May 19, 2010; ReVised Manuscript ReceiVed: June 15, 2010

Starting from the Nakatsuji theorem, a hierarchy of approximations is considered that begins with traditional coupled cluster theory with singles and doubles (CCSD) and proceeds via the ansatz of semigeneralized singles and doubles (CCSGSD), with operators of the types a˜ab ˜ ka ic and a ij included, to the generalized singles pq p and doubles (CCGSD) ansatz with the full basis {a˜rs ; a˜q}, and if necessary, beyond. The simplest realization of CCSGSD is the constrained coupled cluster ansatz with singles, doubles, and triples (CCCSDT). It is related to traditional coupled cluster with singles, doubles, and triples (CCSDT), but it requires (for sufficiently large n) a smaller number of parameters of O(nm3) vs O(n3m3) in CCSDT, n being the number of electrons and m the size of the basis. A linear system of equations for the CCCD[T] approximation is derived and strategies toward its solution are discussed. I. Introduction Electron correlation is mainly related to electron pairs, and whatever formalism one uses, the total correlation energy can always be written as a sum of pair energies (and possibly orbital contributions). Virtual excitations of three or more particles at a time affect the correlation energy only indirectly. Nevertheless, accurate correlation energies are obtained only if one takes care of n-particle correlations in a hierarchical way. In coupled cluster (CC) theory (for a recent review that contains an introduction to the formalism used here; see ref 1), this hierarchy is manifest in the use of k-particle excitation operators. For a closed-shell reference state Φ with the spin orbitals ψi, ψj occupied and ψa, ψb virtual, one uses the basis operators

aia;

aijab;

abc aijk ;

(1)

etc.

These can be expressed through the creation and annihilation operators for a spin orbital ψp, ap† ) ap, and ap, respectively, as

apq ) apaq;

pq ars ) aqaparas;

etc.

(2)

The use of this basis is very systematic but has one drawback, namely, an unfortunate scaling with the number n of electrons and the size m of the one-electron spin-orbital basis. At the level of k-particle excitations the number of parameters is of the order O[nk(m - n)k], where we ignore combinatorial factors. If we further assume that m . n, the dependence is of O[nkmk]. We know, on the other hand, that all information on a quantum mechanical system is in the Fock space Hamiltonian H (11), and all information on a state in its 2-particle density pq . The number of parameters matrix γ2 with matrix elements γrs that specify H or γ2 is of the order m4, which means that CC theory is oVerparametrized, i.e., describes the state by more parameters than are ideally needed to characterize the state. There is obviously a substantial redundancy in the CC ampli†

Part of the “Klaus Ruedenberg Festschrift”. * Corresponding author. Electronic address: werner.kutzelnigg@ ruhr-uni-bochum.de.

tudes, but even more so, e.g., in the coefficients of a full-CI wave function. The concept of minimal parametrization has been outlined elsewhere2 and we refer the reader for details to this study. Attempts to formulate many-electron theory in a minimal parametrization, i.e., in terms of O(m4) parameters2 have not yet been very satisfactory. The minimization of the energy as a functional of the two-particle density matrix γ2 meets the problem that the so-called n-representability conditions are not expressible in terms of γ2 alone and require the γk of higher particle rank k.3 The formulation of the theory in terms of γ1 and the two-particle cumulant λ2 looks somewhat more promising.4 Another interesting approach is that based on the Nakatsuji theorem5 and the Nooijen conjecture,5 where the wave operator is formulated in minimal parametrization.2 Although the Nooijen conjecture is not true (at least not proven), in the sense that it would lead to an exact solution in minimal parametrization,2 there is evidence that it can be the basis for very accurate results, though, so far, not in a very cheap way.7 The approach presented here is related to the Nooijen conjecture and the Nakatsuji theorem, but it proceeds in a hierarchical way, which makes a direct comparison with CC theory straightforward. II. General Formulation The most general formulation of an n electron wave function Ψ, compatible with the separability theorem1,8,9 is the exponential ansatz

Ψ ) eSΦ;

S)

∑ sµΞµ

(3)

µ

where W ) eS is the wave operator, Φ is a separable normalized reference function, e.g., a Slater determinant or an MC-SCF function, {Ξµ} is a basis of particle-number conserving Fock space operators, e.g., of the type (1), and sµ is an expansion coefficients. The energy expectation value is then

10.1021/jp104568g  2010 American Chemical Society Published on Web 06/29/2010

8914

J. Phys. Chem. A, Vol. 114, No. 33, 2010

Kutzelnigg



E)

〈Φ|eS HeS |Φ〉 S†

〈Φ|e e |Φ〉 S

(4)

and the condition for stationarity with respect to variation of the expansion coefficients sµ is

{ } †

0 ) Re〈Φ|

deS (H - E)eS |Φ〉 dsµ

(5)

Let us, for a moment, assume that the basis operators Ξµ commute. Then we get two alternative expressions for the stationarity condition (5) †

0 ) Re〈Φ|Ξµ† eS (H - E)eS |Φ〉 †

0 ) Re〈Φ|eS Ξµ† (H - E)eS |Φ〉

(6) (7)

In the first of these formulations (6) it is recommended to decompose the operator basis Ξ into two subsets, an actiVe one and an inactiVe one, such that inactive basis operators annihilate Φ, while active operators have a nonzero action on Φ. Equation 6 is automatically fulfilled for inactive operators; the energy is hence stationary with respect to variation of their coefficients, such that one need only consider actiVe operators. For a closed shell reference state Φ as specified in the beginning of section I, the operators aai , aab ij , etc. are active, while, e.g., aba or aab ic is inactive. We arrive so at the operator basis (1) familiar from coupled-cluster (CC) theory, although here in the context of variational CC, rather than in the simpler traditional CC (TCC).10 Since the basis operators (1) commute, the derivation of (6) is justified a posteriori. Before we come to the second formulation (7), for which the distinction between active and inactive operators is less helpful, let us note that the derivation of both (6) and (7) is triVial for a basis of commuting operators but requires only that the basis is Lie algebra.2 This is shown in Appendix A. III. Nakatsuji Theorem and Nooijen Conjecture Assume that an operator basis {Ωµ} is chosen such that the Hamiltonian H can be expanded in this basis. (At the moment the wave operator and its expansion in a basis does not concern us.)

H)

∑ hµΩµ

(8)

µ

We assume further that the contracted Schro¨dinger equation5,11

〈Ψ|Ωµ† (H - E)|Ψ〉 ) 0

(9)

is satisfied for all Ωµ in the expansion (8). Then it is trivial to show that

〈Ψ|(H - E)2 |Ψ〉 ) 0

(10)

which is only possible if (H - E)Ψ ) 0. This is the Nakatsuji theorem5 in its most general form. For the Fock space Hamiltonian9,12,13

1 pq rs H ) hpqaqp + Vrs apq 2

(11)

the operator basis {Ωµ} to be chosen is obviously that of all one-particle and two-particle excitation operators

apq,

pq ars

(12)

with all labels of the chosen spin-orbital basis, irrespectively of whether these spin-orbitals are occupied or unoccupied in Φ. If then (9) is fulfilled, Ψ is exact, which means in the present context that it is a full-CI wave function. The contracted Schro¨dinger equation (9) is N conditions, if N is the dimension of the basis {Ωµ}. There is the challenge to formulate the wave operator W in terms of the same basis, with N parameters associated with the individual basis operators, to arrive at a system of N equations for N unknowns. The simplest possible separable ansatz for Ψ in this philosophy was considered by Nooijen6 and Nakatsuji,14 namely, the exponential formulation in analogy with (3), but in terms of the basis (12), rather than that (1) of CC theory. This ansatz has been called CCGSD,14 for coupled cluster with generalized singles and doubles. This name is somewhat misleading, because this ansatz is fundamentally different from that of CC, and rather represents a step in an alternatiVe hierarchy of approximations. However, this hierarchy has not yet been realized, and its formulation is the main concern of the present paper. Nooijen6 has inserted this ansatz into the contraced Schro¨dinger equation (9) and has so arrived at at a nonlinear system of N equations for N unknowns. Nooijen has observed correctly that, if this system had a solution, this would be the exact (i.e., full-CI) wave function, but he was not able to prove that a solution exists. When one refers to the Nooijen conjecture, one usually means the claim that eq 9 with a CCGSD ansatz for Ψ has a solution. Existence of a solution implies then, by virtue of the Nakatsuji theorem, that this solution is exact. After the pioneering studies by Nakatsuji and Nooijen, some formal theoretical papers were published, in view of proving or disproving the Nooijen conjecture.15-17,19 The lack of definite conclusions has, in part, been due to differences in formulating the conjecture, in particular with respect to the role of the reference function Φ. It is relatively easy to show, but physically not meaningful, that one cannot expect the Nooijen conjecture to be true if one formulates it with respect to any completely arbitrary reference function.15,16 It is also obvious, but meaningless as well, that the conjecture is true for a reference function chosen as Φ ) e-SΨ with Ψ exact. One must rather consider ordinary reference functions Φ of the same kind as those that one uses in CC theory, e.g., as a Hartree-Fock type Slater determinant. Then the Nooijen conjecture becomes nontrivial. Some numerical studies.7,18,19 lead to the conclusion that a variational ansatz in the spirit of the Nooijen conjecture, based on an ordinary reference, can lead to, if not to exact, at least to very accurate results. A detailed and conclusive theoretical analysis was given by Mukherjee and the present author,2,20 where all previous studies15-17 were analyzed critically. The main messages2,20 are (a) Equation 7 gives conditions for the stationarity of the energy expectation value 〈Ψ|H|Ψ〉/〈Ψ|Ψ〉 with respect to the variations Ψ f exp(Ωµ)Ψ, but not conditions for the stationarity of (4) with respect to variation of S, because the Ωµ of the CCGSD ansatz do not constitute a Lie algebra. Only in this case does (7) follow from (5).

Alternative Hierarchy of Electron Correlation

J. Phys. Chem. A, Vol. 114, No. 33, 2010 8915

(b) Equation 7 has an undesired singular solution of the type19

Ψ ) WΦ;

W ) lim eλH

(13)

λf-∞

in which the wave operator becomes a projection operator, and where the reference function can be completely arbitrary. There are ways to avoid this undesired solution, e.g., in terms of perturbation theory2 (because the singular solution has no perturbation expansion) or a variational ansatz.7,19 (c) If one is able to discard singular solutions, the existence of other solutions of eq 7 cannot be proven. In a perturbation expansion only the leading order exists.2 This does, however, not exclude that approximate solutions based on the variational method with the CCGSD ansatz can be rather accurate. (d) A key to the disproof of the Nooijen conjecture is2 that the Nooijen ansatz cannot be exact to 4th in perturbation theory, because this would require it to satisfy a linear system of O(n3m3) equations for O(nm3) unknowns, while there is no obvious indication that the rank of this system is only of O(nm3). (e) It is important that the basis of the CCGSD ansatz is not closed under commutation. If one can expand the Hamiltonian into an operator basis that constitutes a Lie algebra, an analog of the Nooijen conjecture would hold. The question whether or not the Nooijen conjecture is exact or only very accurate does not affect the content of the present paper. However, a numerical study of the ansatz proposed here could, in view of a possible equivalence of CCD[T] and CCCD[T] (that has not yet been excluded rigorously), have some impact on the understanding of the Nooijen conjecture.

pq 〈Φ|a˜pq |Φ〉 ) 〈Φ|a˜rs |Φ〉 ) 0

(18)

For the special case of Φ a single Slater determinant this normal order is equivalent to that of the traditional particle-hole picture. The present formulation has the advantage that it is easily generalizable to arbitrary reference states.21 In a graphical representation of operators or their expectation values in the particle-hole picture (see Figure 1) a matrix element is represented by a vertex (a point vertex for Hugenholtz type diagrams9,22) and spin orbitals are symbolized by straight or curved, dominantly vertical lines, with upgoing arrows for particles (virtual spin orbitals, labels a, b, c, ...) and downgoing arrows for holes (occupied spin orbitals, lables i, j, k, ...). The order of operators in a product is from right to left in the formula, and from bottom to top in a diagram. Lines that connect vertices represent contractions. Sometimes it is convenient to draw diagrams without arrows, so-called skeletons. They represent all diagrams with the same topology, but with different arrangements of the arrows. A product of two a˜ operators is the normal product plus all normal ordered particle, hole, or combined contractions, e.g., pr a˜pqa˜sr ) a˜qs + (δrq - κrq)a˜sp - κspa˜rq + (δrq - κrq)κsp

(19) Only full contractions, which are simple scalars, contribute to an expectation Value with respect to the reference function Φ. The corresponding diagrams have no open fermion lines.

〈Φ|a˜pqa˜sr |Φ〉 ) (δrq - κrq)κsp

IV. Normal Ordering in the Particle-Hole Picture

(20)

κqp,

etc. of the reduced density We need the matrix elements matrices κk corresponding to the reference state Φ, which we choose normalized to unity.

κpq ) 〈Φ|apq |Φ〉 ) npδpq

κ1:

(14)

(Note that Φ is normalized to unity.) Normal products do not contribute to the commutator of normal-ordered operators.21

[a˜pq,a˜sr] ) [apq,asr] ) δrqasp - δsparq ) δrqa˜sp - δspa˜rq + δrqκsp δspκrq

pq κrs

κ2:

)

pq 〈Φ|ars |Φ〉

)

npnq(δrpδsq

-

(15) where np ()0 or 1) is the occupation number of ψp in Φ. δqp is a Kronecker delta. The corresponding quantities referring to the correlated wave function Ψ are the γk with elements γqp etc. From now on we use excitation operators such as a˜pq in normal order with respect to the reference function Φ and symbolize these with a tilde, while excitation operators without a tilde such as aqp are in normal order with respect to the genuine vacuum. The relation between the two types of operators is21

a˜pq ) apq - κpq ) apq - npδpq pq a˜rs

)

pq ars

-

κrpasq

)

pq ars

-

npδrpasq

+

κsparq +

nqδsqarp

+

κrqasp

npδsparq

+

+

-

nqδrqasp

npnq(δrpδsq

The Fock space Hamiltonian in normal order with respect to Φ is

1 pq rs H ) E(0) + fqpa˜qp + Vrs a˜pq 2

(22)

1 pq rs E(0) ) 〈Φ|H|Φ〉 ) hpqκqp + Vrs κpq 2

(23)

pr s j qs fqp ) hpq + V κr;

-

+

pq κrs

-

δspδrq)

(17) Expectation values of a˜ operators with respect to Φ vanish.

pr pr rp j qs V ) Vqs - Vqs ;

H0 ) E(0) + fqpa˜qp

(16) κsqarp

(21)

δspδrq)

(24)

The Einstein summation convention is implied.12 E(0) is the energy expection value with respect to Φ. fqp is a matrix element of the Fock operator. The normal order Hamiltonian is automatically in Møller-Plesset partition.1,9 We shall preferentially use a one electron basis in terms of which the Fock operator is diagonal.

fqp ) δpqεp;

H0 ) εpa˜pp ) εia˜ii + εaa˜aa

(25)

8916

J. Phys. Chem. A, Vol. 114, No. 33, 2010

Kutzelnigg

Figure 2. Illustration of constrained triple excitations.

One must further include operators of the typea Y† and U† at an appropriate step. So we arrive eventually at CCGSD. It has then to be decided whether the CCGSD basis can be regarded as sufficiently close to complete, or whether one has to augment it further. One can conjecture that if the basis is not good enough, this is so because it does not constitute a Lie algebra, so one must try to diminish the deviations from a closure under commutation, e.g., include the commutators between important basis operators. This question becomes relevant only after one has succeeded in solving the CCSGSD or the CCGSD equations in a satisfactory way. It is, at this moment, somewhat premature. In the philosophy of the introductory part of this paper one should use the Variational formulation. However, this is tedious already in the context of CCD or CCSD, from which we want to start. Although one will possibly follow this route in the future, we now present a formalism in intermediate normalization, which allows us to take advantage of all the knowhow familiar from traditional coupled-cluster theory. This leads us to an approach proposed a few years ago, but not yet applied, namely, constrained CC theory (CCC).2

Figure 1. Types of two-particle operators.

VI. Formulation in Intermediate Normalization V. The New Hierarchy Previous numerical studies based on the CCGSD ansatz were performed in a brute force way,7,18,19 which is very laborious and has some pitfalls.2 Here we propose a hierarchical approach instead. Rather than expanding the operator S in the basis (12) from the very beginning, we start with the basis a˜ab ij familiar from CCD (coupled-cluster with doubles) and continue by including a˜ai as in CCSD (coupled-cluster with singles and doubles). We begin on safe grounds. Next we augment the operator basis by considering

a˜icab;

a˜ijka

(26)

We use the letter Y for these operators, because their Hugenholtz diagrams (see Figure 1) look like the letter Y.2 For the same reason operators of type (1) are associated with the letter U, and operators (27) with the letter X. One may call the formulation in terms of U and Y a coupledcluster theory with semigeneralized singles and doubles (CCSGSD). In a later step we also allow for basis operators of the type X such as

a˜ab cd ;

a˜ijkl;

jb a˜ia ;

a˜ab;

a˜ij

(27)

Remember that we consider the special case that the reference function Φ is a closed shell single Slater determinant, normalized as 〈Φ|Φ〉 ) 1 with the spin orbitals ψi occupied and ψa unoccupied. We want to apply the CCSGSD wave operator eS with

S ) U + Y;

U ) Uijaba˜ijab + Uiaa˜ia; Y ) Yicaba˜icab + Yijkaa˜ijka (28)

to Φ. One can show (see Appendix B and ref 2) that, since YΦ ) 0 (see below)

eSΦ ) eRΦ;

1 1 R ) U + [Y,U] + [Y,[Y,U]] + ...; 2 6 〈Φ|eR |Φ〉 ) 1 (29)

Obviously, eR looks like a CC ansatz, in intermediate normalization, in which U plays the same role as S2 in traditional CC, while R3 ) 1/2[Y,U] obviously represents triple excitations like S3. The operators that sum up to R, such as U and R3 all commute with each other (even if R3 is constructed from two noncommuting operators).

Alternative Hierarchy of Electron Correlation

J. Phys. Chem. A, Vol. 114, No. 33, 2010 8917

1 R4′ ) [[X,U],U] 6

(37)

VII. Coupled Cluster with Doubles and Triples

Figure 3. Four-particle contributions to R.

One sees easily, as is illustrated in Figure 2, that typical members of R3 are ebd aed [a˜de ˜ ijab] ) δac aijk + δbc aijk kc ,a

(30)

abd adb [a˜klmd,a˜ijab] ) -δimakjl - δjmailk

(31)

However, while S3 is expressed in terms of O(n3m3) parameters, R3 is (for U given) counted by only O(nm3) parameters, which matters for n sufficiently large (see section VIII). It has to be determined to which extent this more compact parametrization simplifies the calculation without affecting too much the final result. For the derivation of (29) it is essential that YΦ ) 0. This is the case if the Y type operators are in normal order with respect to Φ. The operator a˜ab ic does not differ from the corresponding operator without a tilde, and we get

Let us, to simplify the formalism, ignore single excitations, and treat triple excitations only noniteratively, i.e., first CCD and then CCCD[T], the CCC counterpart of CCD[T]. We formulate traditional CCD and CCD[T] in terms of S2 and S3 in a way that is convenient for a generalization to CCCD[T]. The energy is obtained as

E ) E(0) + 〈V + [V,U]〉;

U)

∑ UµΩUµ

(38)

µ

We use the notation

〈...〉 ) 〈Φ|...|Φ〉

(39)

ΩUµ is a double excitation operator of the set U, i.e., of the type a˜ab ij , with µ collecting the labels i, j, a, b (without double counting). If the Hamiltonian H0 of eq 25 is specified, there is, for each basis operator Ωµ, an energy difference, e.g.,

e(a˜ijab) ) εa + εb - εi - εj

(40)

At CCD level we get U from

a˜icabΦ ) abaaaiacΦ ) 0

(32)

because ψc is not occupied in Φ and, hence, cannot be annihilated. For a˜kb ij the situation is a little more complicated.

a˜ijkbΦ ) (abakaiaj - δikabaj + δjkabai)Φ

(33)

If k * i and k * j, the second and third terms in parentheses are absent, the action of the first term vanishes, because

abakaiajΦ ) -abaiakajΦ ) abaiajakΦ ) 0

(34)

since ψk is occupied in Φ and cannot be created again. If k ) i, we get

a˜ijibΦ

) (a a aiaj - a aj)Φ ) 0 b i

b

(35)

since the action of the two terms cancels. The case k ) j is analogous. It is suggestive to refer to a CC ansatz based on eR as a constrained coupled cluster (CCC) ansatz2 and specifically the ansatz truncated at R3 as CCCSDT. The next step would be to include R4. However, there is also a 4-particle term that originates if one goes one step further in the new hierarchy,2 i.e., for

S ) U + Y + X; kl ij jb ia X ) Xcd ˜ ab ˜ kl + Xia a˜jb + Xbaa˜ab + Xija˜ji aba cd + Xij a

〈 (

1 † 0 ) ΩUµ [H0,U] + V + [V,U] + [[V,U],U] 2

)〉 (41)

We can rewrite eq 41 as † † -〈[[ΩUµ ,H0],U]〉 ) e(ΩUµ)〈[ΩUµ ,U]〉 ) e(ΩUµ)Uµ ) 1 † ΩUµ V + [V,U] + [[V,U],U] (42) 2

〈 (

Equation 42 is a nonlinear equation for the coefficients Uµ, which we solve iteratively, starting usually with the lowest order of perturbation theory † e(ΩUµ)UUµ ) 〈ΩUµ V〉 ) Vµ;

U ) VH

(43)

where we have expanded V in the full basis of two-electron operators which contains the a˜ab ij as a subset. The solution of (42) can be written in compact form as

1 U ) V + [V,U] + [[V,U],U] 2 H

(

)

(44)

where the subscript H means inversion of the commutator with H0.1,2 A good approximation is often the linear CEPA0 system.

(36) U ) (V + [V,U])H;

The new term is then (see Figure 3)

)〉

† eµUµ ) 〈ΩUµ (V + [V,U])〉

(45)

8918

J. Phys. Chem. A, Vol. 114, No. 33, 2010

Kutzelnigg

Figure 4. Comparison of skeletons for E3 and E˜3. A square dot represents S2, U, and Y, a circular dot represents V, and a broken line is an energy denominator. Arrows are not marked.

Figure 5. Different ways to 3-particle excitations, Type 1.

At the CCD[T] level we have to solve two linear systems of equations † 0 ) 〈Ω3µ ([H0,S3] + [V,S2]3)〉

(46)

〈 (

1 † 0 ) Ω2µ [H0,S2] + [V,S2]2 + [V,S3]2 + [[V,S2],S2] 2 (47)

)〉

with the formal solution Figure 6. Different ways to 3-particle excitations, Type 2.

S3 ) [V,S2]3H

(48) VIII. Combinatorial Considerations

1 S2 ) [V,S2]2 + [V,S3]2 + [[V,S2],S2] 2 H

(

)

(49)

A number as a subscript indicates the particle rank of an operator. The respective second equations (47) and (49) of the two sets, which describe the updating of U to S2, are as important as the corresponding first equations (46) and (48), since S3 enters the energy (38) only via S2. Let us, therefore, pay special attention to the equation for

∆S2 ) S2 - U

(50)

We have, so far, only considered the case that n is in the asymptotic regime in which it matters whether the number of parameters is of O(nm3) or O(n3m3). If one wants to include small n as well, the situation changes. The number of parameters to characterize S3 for arbitrary n is

∆S2 ) [V,S3]2H

(51)

)( )

()

(54)

while

NR )

namely

(

m-n n m3 n ) + O(m2) 3 3 6 3

NS )

(

)

m-n m3n (m - n)n ) + O(m2) 2 2

(55)

parameters are needed for Y. The ratio between these factors is

or

e(Ω2µ)∆S2µ )

† 〈Ω2µ [V,S3]2〉

) [V,S2]3H

(52)

This leads to the triples contribution to the energy

E3 ) 〈[V,∆S2]〉 ) 〈[V,[V,S3]2H]〉 ) 〈[V[V,[V,S2]3H]2H]〉 (53) An important contribution to E3 is shown on Figure 4. Each skeleton displayed in Figure 4 corresponds to two Hugenholtz type diagrams with different positions of the arrows.1,9 Note the presence of two energy denominators. The content of this section has been a recapitulation of traditional CCD[T] in a compact, somewhat unconventional formulation.1

NR /NS )

18 + O(m-1) (n - 1)(n - 2)

(56)

it goes to 18/n2 for large n but has the values 9, 3, 3/2, and 9/10 for n ) 3, 4, 5, 6, respectively. This means that for n < 6 one has even more parameters in CCCD[T] than in CCD[T], and that one cannot gain by CCCDT with respect to CCDT for small n. This is not too surprising, since there are many ways to construct an operator of type S3 in the form [Y, U], as shown on Figures 5 and 6, for the same final upper and lower labels. There are actually 18 possibilities, of which 9 on Figure 5 correspond to the first line of Figures 2 and 9 on Figure 6 to the second line of Figure 2. In CCD[T] one adds all 18 combinations of VY and U, which lead to the same 3-particle excitation, and divides then by the

Alternative Hierarchy of Electron Correlation

J. Phys. Chem. A, Vol. 114, No. 33, 2010 8919

cumulative energy denominator, to get S3. In CCCD[T] one chooses a basis of operators ΩYµ and constructs from these the [ΩYµ, U]. Each of these, for one µ, is a linear combination of generally different 3-particle excitation operators. For the special case of n ) 3, m ) 6 only one 3-particle excitation operator is possible. The 18 operators [ΩYµ, U] are hence all proportional to each other. It would be sufficient to take a single [ΩYµ, U], e.g., for µ ) 1, to expand R3. It is obvious that, at least for small n, the basis operators of CCCD[T] are linearly dependent, and that one must find a systematic way to avoid linear dependencies. In preliminary attempts use of a singular value decomposition can be helpful, but this should be avoided in a practicable approach.

1 ˆ 0,Rˆ3] + [Vˆ,U ˆ ]3)〉 0 ) 〈Zˆ†κ([H 2 To simplify this, we note that

ˆ ]3 ) [Vˆ,U

∆S˜2 ) [V,R3]2H

1 ) [V,[Y,U]]2H 2

(57)

∑ VYλ[Ωˆ Yλ,Uˆ] ) 2VYλZˆλ

(63)

λ

ˆ 0,Rˆ3] ) [H

IX. Constrained Coupled Cluster Theory with Doubles and Triples To go from CCDT to CCCDT, we must replace S3 by R3 ) 1 /2[Y,U]. Let us first assume that we know Y and use this to evaluate ∆S˜2 and E˜3, the counterparts of ∆S2 and E3 studied at the end of the section VIII. We get

(62)

∑ Yλ[Hˆ0,Zˆλ] ) 21 ∑ Yλ[Hˆ0,[Ωˆ Yλ,Uˆ]] λ

λ

1 )2

∑ Yλ([[Hˆ0,Ωˆ Yλ],Uˆ] - [[Hˆ0,Uˆ],Ωˆ Yλ])

)

∑ Yλ(e(ΩYλ)Zˆλ - 21 [[Hˆ0,Uˆ],Ωˆ Yλ])

λ

λ

(64) The condition to be satisfied by the coefficients Yκ is then the linear system

∑ Yλe(ΩYλ)〈Zˆ†κZˆλ〉 - 21 ∑ Yλ〈Zˆ†κ[[Hˆ0,Uˆ],Ωˆ Yλ]〉 ) λ

λ

-2

∑ VYλ〈Zˆ†κZˆλ〉

(65)

λ

1 E˜3 ) 〈[V,∆S˜2]〉 ) 〈[V,[V,[Y,U]]2H]〉 2

(58)

An important contribution to E˜3 is compared with the counterpart of E3 on Figure 4. There is only one explicit energy denominator in E˜3, the other is implicit in Y. One must also consider contributions of a different topology, as displayed in Figures 8-11. Let us, for the rest of this paper, designate operators by a hat, to distinguish them from expansion coefficients. To construct the linear system satisfied by Y, i.e., the counterpart of (46), we expand

(59)

κ

ˆ Yκ basis operators of the type (26) and define with Ω

1 ˆ ˆ Zˆκ ) [Ω ,U] 2 Yκ

(60)

then Rˆ3 is expanded as

∑ Yκ[Ωˆ Yκ,Uˆ] ) ∑ YκZˆκ κ

1 ˆ† ˆ† ˆ ˆ ,ΩYκ][ΩYλ,U]〉 〈Zˆ†κZˆλ〉 ) 〈[U 4 1 ˆ† ˆ† ˆ ˆ U*U ) 4 µ ν〈[ΩUµ,ΩYκ][ΩYλ,ΩUν]〉 µ,ν



(66)

ˆ† ˆ† ˆ ˆ 〈Zˆ†κ[[Hˆ0,Uˆ]Ωˆ Yλ]〉 ) 21 ∑ U*U µ νe(ΩUµ)〈[ΩUµ,ΩYκ][ΩYλ,ΩUν]〉 µ,ν

(67) X. Discussion of the Linear System Eq 65

ˆ Yκ Yˆ) ∑ YκΩ

1 ˆ 1 Rˆ3 ) [Yˆ,U ]) 2 2

with

(61)

κ

The basis operators Zˆκ span a subset (or the full set) of the threeˆ 3µ. particle excitation operators Ω We can, of course, not simply replace Sˆ3 in (46) by Rˆ3, because the number of O(n3m3) equations to be solved would be larger than the number of available coefficients Yκ of O(nm3), at least for large n and m. We must rather require

Both the linear system (46) of CCD[T] theory and that (65) of CCCD[T] theory can be regarded as conditions for stationarity of a Hylleraas functional, as is illustrated on Figure 7. In the case of (46) this functional consists of the diagrams (a) and twice the real part of (b), for (65) analogously of (c) and (d). Diagrams of a different topology (see Figures 8-11) must also be considered. In Figure 7 vertices of type U or S2 are symbolized by open square dots, vertices for V by full spherical dots, and vertices for Y or S3 by full square dots. A thick horizontal line means multiplication by an orbital energy difference. The algebraic equivalents of the diagrams in Figure 7 are (a) 〈S3†[H0 S3]〉, (b) 〈S3†[V,U]3〉, (c) 〈R3†[H0,R3]〉, and (d) 〈R3†[V,U]3〉. The linear system (46) of CCD[T] theory is of the simple form

Ax b)b b

(68)

with A a diagonal matrix, which makes the solution particularly simple. Just the dimension is very large. While the system (46) is regular, the linear system (65) is likely to be ill-conditioned, as discussed in the previous section.

8920

J. Phys. Chem. A, Vol. 114, No. 33, 2010

Kutzelnigg

Figure 7. Illustration of the Hylleraas functionals: (a) and (b) for CCD[T]; (c) and (d) for CCCD[T]. Figure 10. Diagrams for expectation values 〈[U†, Y†][Y, U]〉, skeleton type b.

Figure 8. Skeletons for expectation values 〈[U†, Y†][Y, U]〉.

Figure 11. Diagrams for expectation values 〈[U†, Y†][Y, U]〉, skeleton type d.

Figure 9. Diagrams for expectation values 〈[U , Y ][Y, U]〉, skeleton types a and c. †



Let us now assume that we have somehow solved this problem and achieved that the basis operators are linearly independent. Although the matrix of the linear system (65) certainly has also nondiagonal elements, it is reasonable to assume that the diagonal terms dominate. Let us, therefore, for a moment ignore the terms with κ * λ. If we further neglect the second term on the left-hand side of (65), we get a very simple solution

Yλ ) 2{e(ΩYλ)}-1VYλ;

Yˆ ) 2(VˆY)H

(69)

We can compare this with the result for the first-order approximation U(1) for U

Uµ(1) ) {(ΩUµ)}e-1VUµ;

ˆ ) (VˆU)H U

(70)

which is very similar, apart from the factor 2. At this level of approximation we get

1 R3 ) [Y,U] ) 2

∑ VYκ{e(ΩYκ)}-1 ) [VH,U]3

(71)

κ

instead of the correct CCD[T] result

S3 ) [V,U]3H

(72)

In this very crude approximation we replace the cumulatiVe energy denominator present in S3, which makes the energy nonfactorizable, by a local energy denominator, associated with only the vertex VY, such that the diagram for R3 becomes factorizable. It is plausible, but also good news, that Yˆ is mainly the direct response to VˆY. In a second step we can, keeping κ ) λ, take care of the neglected term in (65) iteratively. We need to consider

1 〈Zˆ†κZˆκ〉 ) 4 1 ) 4

ˆ† ˆ† ˆ ˆ ∑ U*U µ µ〈[ΩUµ,ΩYκ][ΩYκ,ΩUµ]〉 µ

∑ U*U µ µ µ

(73)

Alternative Hierarchy of Electron Correlation

〈Zκ†[[Hˆ0,Uˆ]ΩYκ]〉

) )

1 2

∑ U*U e(Ω

1 2



µ

µ

J. Phys. Chem. A, Vol. 114, No. 33, 2010 8921

ˆ† ˆ† ˆ ˆ Uµ)〈[ΩUµ,ΩYκ][ΩYκ,ΩUµ]〉

µ

† U* µ Uµe(ΩUµ)

µ

(74) Let us define an average pair excitation energy

κ )

ˆ 0,U ˆ ],Ω ˆ Yκ]〉 〈Zˆ†κ[[H 〈Zˆ†κZˆκ〉

)

† ∑ U*U µ µe(ΩUµ) µ

∑ U*U µ µ

(75)

µ

Then we can write (65), for κ ) λ as

Yκ(e(ΩYκ) + εκ)〈Zˆ†κZˆκ〉 ) -2VYκ〈Zˆ†κZˆκ〉

(76)

We get a shift of the energy denominator of Yˆκ, which simulates to some extent the cumulative denominator, which arises in CCD[T], but keeps R3 factorizable. After this, one will have to take care of the terms with κ * λ. It is too early to devise an optimal strategy. It is, however, useful to show the diagrams that illustrate the contributions to the expressions 〈Zκ†Zκ〉. In Figure 8 we show the 4 possible skeletons, of which (a) is probably the most important, because it contains the diagonal terms with κ ) λ. There are 2 Hugenholtz diagrams both for skeletons (a) and (c) and 6 Hugenholtz diagrams both for skeletons (b) and (d), i.e., a total of 16. One has to find out which diagrams can possibly be neglected, or whether one can use the Hylleraas functional, but make only the bulk of it stationary. The following observation is useful. A typical contribution in CCD[T] or MP4 is of the type (a) lc UijabVjckmVkm Uilab(εa + εb + εc - εi - εk - εm)-1

(77)

XI. Conclusions This paper is a step toward minimal parametrization of an n-electron state. It presents an alternative to the coupled-cluster hierarchy and is based on the Nakatsuji theorem.5 While the change of paradigm is novel, the lowest-order approximation in the new context has already been suggested a few years ago as constrained coupled-cluster (CCC),2 though in a somewhat hidden and probably too modest form. None of the experts of computer implementation has so far been challenged to consider CCC. There are strong physical arguments in favor of CCC. Triple excitations with respect to the reference function are mathematical artifacts and can hardly be interpreted in simple physical terms. Correlation is the response to the interaction. Since there are no 3-particle interactions in the Hamiltonian, a 3-particle correlation is hard to understand, moreover, since it affects the energy only indirectly. It is, however, physically plausible that, in a first step, namely, in electron pair theory, one adds pair correlation upon the reference function, while in a second step one puts additional pair correlations onto the lowest order pair correlated function, as one does in CCCDT. Here triple excitations appear as double excitations to double excitations. So CCCDT should not miss much physics. On the other hand, the computational effort of CCCD[T] should ideally scale as n5 with the system size, while CCD[T] scales with n7. It is a worthwhile challenge to implement CCCD[T] such that its scaling comes close to the ideal one. It must also be admitted that there are some unsolved problems, mainly those related to possible linear dependencies of excitation operators. Further steps in the hierarchy proposed here should be postponed until sufficient experience has been gained with CCCD[T] and CCCDT. Acknowledgment. This paper is dedicated to Klaus Ruedenberg, the great pioneer of quantum chemistry, who has remained so young. The author thanks Jozef Noga and Volker Staemmler for valuable comments, and Peter Knowles for discussions on this topic. Appendix A. The GCCSD Stationarity Condition and the Nakatsuji Theorem

We must first sum

∑ VkmlcUilab

(78)

Consider the expression

l

e-S



which is an n7 step (if we assume that m is proportional to n), before we can divide by the cumulative energy denominator. If, however, we have to evaluate an expression without a cumulative denominator such as lc UijabVjckmYkm Uilab

∑ UijabUilab;

∑ VjckmYkmlc

)

1 1 d + Ωµ† + [Ωµ† ,S†] + [[Ωµ† ,S†],S†] + ... dsµ 2 6

)

d ˜ µ† +Ω dsµ

(79)

we can first sum

i,a,b

[ ] [[ ] ] [[[ ] ] ]

d 1 d † † d S† d † e ) + ,S + ,S ,S + dsµ dsµ dsµ 2 dsµ 1 d † † † ,S ,S ,S + ... 6 dsµ

(80)

(81)

We conclude that

k,m,c

which are two n5 steps (again if m is proportional to n). One should take advantage of this as much as possible.

[

]

† d S† ˜ µ† ,e ) eS Ω dsµ

(82)

8922

J. Phys. Chem. A, Vol. 114, No. 33, 2010

Kutzelnigg





0 ) 〈eS |eS〉

† d 〈eS |H|eS〉 ) 〈eS Ωµ† |(H - E)|eS〉 dsµ 〈eS† |eS〉

One concludes immediately that

˜ µ† ) Ωµ† + [Ωµ† ,S†] + 1 [[Ωµ† ,S†],S†] + Ω 2 1 [[[Ωµ† ,S†],S†],S†] + ... (84) 6 ˜ µ† are If the operators Ωµ† span a Lie algebra, the operators Ω members of this Lie algebra and hence linear combinations of the Ωµ† . Therefore, eqs 6 and 7 follow from eq 5. Appendix B. Re-formulation of the Wave Operator The re-formulation of the wave operator from the form eS to eR is a special case of the application of a generalized Wick theorem to a product of operators, which leads to a sum of the respective normal product and all normal-ordered contractions.21 Inactive operators, such as Y, which annihilate Φ can only survive in active contractions. For the present purpose a pedestrian approach is possible. Expand exp(U + Y) in a power series and consider the first terms: (U + Y)Φ ) UΦ 1 2 1 1 2 (U + Y) Φ ) (U + UY + YU + Y2)Φ ) (U2 + [Y,U])Φ 2 2 2 1 (U + Y)3Φ ) 1 3 6 (U + U2Y + UYU + YU2 + UY2 + 6 YUY + Y2U + Y3)Φ )

1 3 (U + 3U[Y,U] + [Y,[Y,U]])Φ 6

1 1 R ) U + [Y,U] + [Y,[Y,U]] + ... 2 6

(83)

(85)

We have frequently used that YΦ ) 0 and noted that U commutes with [Y, U].

(86)

References and Notes (1) Kutzelnigg, W. Unconventional aspects of coupled-cluster theory. In Recent Progress in Coupled-Cluster Methods; Carski, P., et al., Eds.; Springer, Berlin 2010. (2) Kutzelnigg, W.; Mukherjee, D. Phys. ReV. A 2005, 41, 022502. (3) Nakata, M.; Nakatsuji, H.; Ehara, M.; Fukuda, M.; Nakata, K.; Fujisawa, K. J. Chem. Phys. 2001, 114, 8282. Mazziotti, D. A. Phys. ReV. Lett. 2004, 93, 213001. Mazziotti, D. A. J. Chem. Phys. 2004, 121, 10957. Zhao, Z. J.; Braams, B. J.; Fukuda, M.; Overton, M. L.; Percus, J. K. J. Chem. Phys. 2004, 120, 2095. (4) Kutzelnigg, W. J. Chem. Phys. 2006, 125, 171101. (5) Nakatsuji, H. Phys. ReV. A 1976, 14, 41. (6) Nooijen, M. Phys. ReV. Lett. 2000, 84, 2108. Nooijen, M. J. Chem. Phys. 2003 118, 4832. (7) Piecuch, P.; Fan, P.-D. J. Mol. Struct. THEOCHEM 2006, 768, 3. (8) Primas, H. In Modern Quantum Chemistry; Sinanoglu, O., Ed.; Academic Press, New York, 1965; Vol 2, p 45. (9) Kutzelnigg, W. Int. J. Quantum Chem. 2009, 109, 3858. (10) Kutzelnigg, W. Theor. Chim. Acta 1991, 80, 349. KutzelniggW. Chem. Phys. 1998, 94, 65. (11) Cohen, L.; Freshberg, C. Phys. ReV. 1976, A13, 927. (12) Kutzelnigg, W. J. Chem. Phys. 1982, 77, 3081. Kutzelnigg, W.; Koch, S. J. Chem. Phys. 1983, 79, 4315. Kutzelnigg, W. J. Chem. Phys. 1984, 80, 822. Kutzelnigg, W. J. Chem. Phys. 1985, 82, 4166. (13) Kutzelnigg, W. Quantum Chemistry in Fock Space. In Aspects of Many-Body Effects in Molecules and Extended Systems; Mukherjee, D., Ed.; Lecture Notes in Chemistry 50; Springer: Berlin, 1989. (14) Nakatsuji, H. J. Chem. Phys. 2000, 113, 2949. Nakatsuji, H. J. Chem. Phys. 2001, 115, 2465. (15) Ronen, S. Phys. ReV. Lett. 2003, 91, 123002. (16) Davidson, E. R. Phys. ReV. Lett. 2003, 91, 123001. (17) Mazziotti, D. A. Phys. ReV. A 2004, 69, 012507. (18) van Voorhis, T.; Head-Gordon, M. J. Chem. Phys. 2001, 115, 5033. (19) Piecuch, P.; Kowalski, K.; Fan, P.-D.; Jedziniak, K. Phys. ReV. Lett. 2003, 90, 113001. (20) Mukherjee, D.; Kutzelnigg, W. Chem. Phys. Lett. 2004, 397, 174. (21) Kutzelnigg, W.; Mukherjee, D. J. Chem. Phys. 1997, 107, 432. (22) Hugenholtz, N. M. Physica 1957, 23, 481.

JP104568G