Feature Article pubs.acs.org/JPCA
Second-Order Many-Body Perturbation Theory: An Eternal Frontier So Hirata,* Xiao He,† Matthew R. Hermes, and Soohaeng Y. Willow‡ Department of Chemistry, University of Illinois at Urbana−Champaign, 600 South Mathews Avenue, Urbana, Illinois 61801, United States ABSTRACT: Second-order many-body perturbation theory [MBPT(2)] is the lowest-ranked member of a systematic series of approximations convergent at the exact solutions of the Schrödinger equations. It has served and continues to serve as the testing ground for new approximations, algorithms, and even theories. This article introduces this basic theory from a variety of viewpoints including the Rayleigh− Schrödinger perturbation theory, the many-body Green’s function theory based on the Dyson equation, and the related Feynman−Goldstone diagrams. It also explains the important properties of MBPT(2) such as size consistency, its ability to describe dispersion interactions, and divergence in metals. On this basis, this article surveys three major advances made recently by the authors to this theory. They are a finite-temperature extension of MBPT(2) and the resolution of the Kohn−Luttinger conundrum, a stochastic evaluation of the correlation and self-energies of MBPT(2) using the Monte Carlo integration of their Laplace-transformed expressions, and an extension to anharmonic vibrational zero-point energies and transition frequencies based on the Dyson equation.
■
INTRODUCTION Second-order many-body perturbation theory,1,2 MBPT(2), is the lowest-ranked member of a systematic series of approximations formally convergent at the exact solutions of the electronic Schrödinger equations. It is also the simplest ab initio molecular orbital (MO) theory3−5 that can describe covalent, ionic, hydrogen-bond, and dispersion interactions accurately and from first principles. MBPT(2), therefore, enjoys special importance in electronic structure theory and increasingly in other areas of computational chemistry such as vibrational structure theory,6−8 serving as the testing ground for new approximations, algorithms, and even theories. There are many facets to this simple and fundamental theory, which are reviewed in the introduction. In the subsequent sections, we describe three major advances we have recently made to MBPT(2), namely, a finite-temperature extension of MBPT(2) and the resolution of the Kohn−Luttinger conundrum,9 a stochastic evaluation of MBPT(2) correlation energies and self-energies,10−12 and anharmonic vibrational MBPT(2) based on the Dyson equation.13 Rayleigh−Schödinger Perturbation Theory. The easiest way to introduce MBPT(2) is to use the Rayleigh−Schrödinger perturbation theory (RSPT).3 We partition the Hamiltonian (Ĥ ) into the zeroth-order part (Ĥ 0) and the perturbation (V̂ ) and expand the exact energy and wave function as infinite sums of progressively smaller, higher-order perturbation corrections: (1) (2) Ψ0 = Φ(0) 0 + Φ0 + Φ0 + ...
E0 =
E0(0)
+
E0(1)
+
E0(2)
+ ...
(0) (0) Ĥ 0 Φ(0) n = En Φn
Substituting these in the Schrödinger equation (with Ĥ ) and equating terms of comparable sizes, we obtain the general expression for E(2) 0 , which reads E0(2) =
∑
̂ (0) (0) ̂ (0) ⟨Φ(0) 0 |V |Φn ⟩⟨Φn |V |Φ0 ⟩ E0(0) − En(0)
n>0
(4)
By applying this to ab initio MO theory with the Møller− Plesset (MP) partitioning1 of the exact electronic Hamiltonian, we obtain E0(2) =
1 4
occ. vir.
∑∑ i,j
a,b
⟨Φ0|V̂ |Φijab⟩⟨Φijab|V̂ |Φ0⟩ ϵi +ϵj −ϵa −ϵb
(5)
as the second-order correction to the energy of the Hartree− Fock (HF) or self-consistent-field (SCF) theory. The sum of E(0) 0 and E(1) 0 is equal to the HF energy. Here, i and j run over all occupied HF spin−orbitals, a and b over all virtual spin−orbitals, Φ0 denotes the HF wave function, Φab ij is the two-electronexcited Slater determinant thereof, and ϵp is the energy of the pth canonical spin−orbital. Evaluating the factors in the numerator with, say, the Slater−Condon rules,3 we obtain E0(2) =
(1)
1 4
occ. vir.
∑∑ i,j
a,b
⟨ij || ab⟩⟨ab || ij⟩ ϵi +ϵj −ϵa −ϵb
(6)
where ⟨pq∥rs⟩ is an antisymmetrized two-electron integral in physicists’ notation.3 This is the canonical, spin−orbital energy expression of MBPT(2), which is familiar to many quantum
(2)
The subscript labels a state, with zero being the ground state, the parenthesized superscript denotes the perturbation order, and (0) Φ(0) n and En are the solutions of the zeroth-order Schrödinger equation with Ĥ 0, i.e., © 2013 American Chemical Society
(3)
Received: October 25, 2013 Revised: December 6, 2013 Published: December 12, 2013 655
dx.doi.org/10.1021/jp410587b | J. Phys. Chem. A 2014, 118, 655−672
The Journal of Physical Chemistry A
Feature Article
Hence, the modern derivation of MBPT(2) is simply to draw the one diagram in Figure 1 and directly interpret it algebraically, according to the established rules found in the literature.2,3,5 This goes as follows: The upper and lower vertexes in Figure 1 are interpreted as the antisymmetrized two-electron integrals, ⟨ij∥ab⟩ and ⟨ab∥ij⟩, respectively. A resolvent, (ϵi + ϵj − ϵa − ϵb)−1, is attached to each pair of adjacent vertexes. The orbital energies in it (with the plus and minus signs for occupied and virtual orbitals) are associated with the edges found in between the pair. Summations must be taken over all internal edges (edges whose two ends are attached to vertexes). Finally, a factor of (−1)h+l(1/2)e must be multiplied, where h, l, and e are the numbers of occupied-orbital (hole) edges, loops, and pairs of equivalent edges, respectively (two edges are equivalent when they start and end with the same vertexes and have the same occupied/virtual attribute). In this diagram, h = 2, l = 2, and e = 2. This interpretation again leads to eq 6, as it should. Quantum chemists at the forefront of MBPT and CC developments can, in fact, conceive new theories and devise efficient algorithms entirely diagrammatically.15 Spin Integration. For a closed-shell molecule (with an even number of electrons), we can rewrite the above spin−orbital expression into spatial−orbital one. It then becomes
chemists. However, derived in this way using RSPT, it may be more suitably referred to as the second-order Møller−Plesset perturbation (MP2) theory.3 Many-Body Perturbation Theory. The reliance on the Slater−Condon rules and the first quantization becomes almost an insurmountable obstacle for treating anything beyond MBPT(2), including the related coupled-cluster (CC) theory.5,14,15 We must switch to more powerful mathematical tools furnished by quantum field theory (QFT).5,16−18 The use of these tools distinguishes MBPT(2) from MP2 or, generally, MBPT from RSPT. A derivation based on QFT can ensure size consistency, which is a crucial requirement of any many-body theory applicable to atoms, molecules, and solids as well as everything in between in size. There are three essential QFT tools, which are, in the increasing degree of sophistication and expediency, the second quantization,3 the normal-ordered second quantization and Wick’s theorem,5,16 and the Feynman−Goldstone diagrams.5,16,18 Using the normal-ordered second quantization, the task of evaluating the factors in the numerator of eq 5 is replaced by a systematic process of drawing all possible Wick’s contraction lines between normal-ordered creation and annihilation operators. This process can even be automated.19 For instance,
occ. vir.
E0(2) = 2 ∑ ∑ ∼ ∼ i ̃, j a∼ , b
̃ ab̃ |̃ i j̃ ̃ ⟩ ⟨i j̃ ̃ |ab̃ ⟩⟨ − ϵ i ̃+ϵ j ̃−ϵa ̃−ϵb ̃
occ. vir.
∑∑ ∼ ∼ i ̃, j a∼ , b
⟨i j̃ ̃ |bã ⟩̃ ⟨ab̃ |̃ i j̃ ̃ ⟩ ϵ i ̃+ϵ j ̃−ϵa ̃−ϵb ̃ (8)
where ⟨p̃q̃|r̃s⟩̃ is a two-electron integral (not antisymmetrized) over spatial MOs, which are distinguished from spin−orbitals by tildes. This is the programmable equation of MBPT(2) using the spin-restricted HF reference wave function.3 It can be obtained by a spin-integration of eq 6 or by direct algebraic interpretation of Goldstone diagrams3 shown in Figure 2: the left and right diagrams correspond to the first and second terms of eq 8, respectively.
This leads back to the same canonical MBPT(2) expression given by eq 6, as it should. A few exercises of drawing Wick’s contraction lines in MBPT or CC derivations tend to lead to the realization that the act of drawing a contraction line between two operators is equivalent to the act of connecting two vertexes with an edge (a line) in the corresponding Feynman−Goldstone diagram. 5,16 The MBPT(2) diagram in the Hugenholtz style is shown in Figure 1. The four contraction lines in each term of eq 7 correspond to the four edges in this diagram. Other than the fact that the single diagram can encompass all four patterns of contraction in eq 7, there is one-to-one mapping between the normal-ordered second quantized and diagrammatic derivations of MBPT(2).
Figure 2. Second-order Goldstone diagrams (type N).
Integral Transformation. This expression (eq 8) may appear to suggest O(n4) dependence of the operation cost of MBPT(2) on the number of orbitals (n). This is not the case. The true bottleneck of an MBPT(2) calculation is the transformation of two-electron integrals from the atomic orbital (AO) to MO basis,3 namely, ⟨i j̃ ̃ |ab̃ ⟩̃ =
̃
̃
̃
̃ b Cμi *Cνj *CκaC λ ⟨μν|κλ⟩
∑ μ,ν ,κ ,λ
=
∑ Cμi ̃*(∑ Cνj ̃ *(∑ Cκa(̃ ∑ Cλb⟨̃ μν|κλ⟩))) μ
ν
κ
(9)
λ
where ⟨μν|κλ⟩ is a two-electron integral over AOs and is an expansion coefficient of the μth AO in the p̃th spatial MO. As indicated by the parentheses, the summation can be performed stepwise, but the operation cost for generating all MO-integrals nevertheless increases as staggering O(n5). The memory cost is Cp̃μ
Figure 1. Second-order Hugenholtz diagram (type N). 656
dx.doi.org/10.1021/jp410587b | J. Phys. Chem. A 2014, 118, 655−672
The Journal of Physical Chemistry A
Feature Article
also large and can scale as O(n4) in the most primitive implementation. Neither the AO- nor MO-integrals can usually fit in a core memory, and therefore, they need to be stored on disks or across many processors’ memories in a parallel machine. Worse still, the AO- and MO-integrals have completely different data layouts: the former according to distances between spatially localized AOs and the latter based on symmetries of delocalized MOs; the data locality is hard to achieve. In short, this is one of the most nightmarish parts of an entire electronic structure program. In a later section, we will describe a whole new algorithm of MBPT(2) that eliminates this integral-transformation step altogether.10−12 Many-Body Green’s Function Theory. Another, closely related viewpoint to MBPT is provided by many-body Green’s function theory (MBGF).16−18,20−24 Here, the central equation is the Dyson equation. In the matrix form, it reads3,20,22 G(ω) = G0(ω) + G0(ω)Σ(ω)G(ω)
Figure 3. Dyson equation. In the diagonal approximation, the doublelined edge represents {G}pp and the single-lined edge {G0}pp.
⎛ 1 ⎞e −iE0(2) = ( −1)l ⎜ ⎟ ⎝2⎠ ∞
∫−∞
(10)
(11)
dω2 2π
dω3 {−i⟨ij || ab⟩}{−i⟨ab || ij⟩} 2π
(14)
{G0(τ )}pq = ∓δpq exp{−(ϵp −μ)τ }
(15)
where μ is the chemical potential, τ is the time difference, and the minus sign is taken if ϵp > μ and the plus sign otherwise. The origin of time is set at the lower vertex of Figure 1 and the upper vertex occurs at time τ > 0. In this picture,18 the MBPT(2) energy is the product of the algebraic interpretations of diagrammatic components summed and integrated over the intermediate quantities, which now include τ:
(12)
which is known as the inverse Dyson equation.3 Substituting an exact one-electron energy ϵ for ω and neglecting off-diagonal elements (i.e., invoking the diagonal approximation), we find ϵ−ϵp = {Σ(ϵ)}pp
∞
∫−∞
where l = 2 and e = 2. It may be noticed that the sum of the frequencies (multiplied by −1 when associated with a hole edge) is conserved at any vertical position of the diagram. The integrations over ω in the above expression can be carried out analytically, as shown in the Appendix. The result is, of course, the same MBPT(2) expression, eq 6. Laplace-Transform MP2. The derivation based on MBGF can also be done using zeroth-order Green’s functions in the realor imaginary-time domain or even in the imaginary-frequency domain.18,26 In the imaginary-time domain, their elements are defined by
This element has a pole at the energy of the pth HF orbital (with the plus sign taken when pth orbital is virtual and the minus sign when it is occupied) and δ is a positive infinitesimal that nevertheless satisfies ∞·δ = ∞.18 The exact and HF Green’s functions are related to each other by the Dyson equation through the self-energy denoted by Σ(ω). All the matrices are dependent on frequency ω. Multiplying G0−1 from the left and G−1 from the right on the Dyson equation, we obtain G0−1(ω) = G−1(ω) + Σ(ω)
dω1 2π
{iG0(ω3)}jj
δpq ω − ϵp ± iδ
i ,j,a,b
∫−∞
{iG0(ω1 + ω2)}aa{iG0(ω2)}ii {iG0( −ω1 + ω3)}bb
where G is the exact Green’s function, which has poles at exact one-electron energies (correlated ionization or electron-attachment energies) and G0 is the zeroth-order Green’s function (in our case, the HF Green’s function) given by its elements18 as {G0(ω)}pq =
∞
∑
(13)
⎛ 1 ⎞e −E(2) = lim ( −1)l ⎜ ⎟ ⎝2⎠ β →∞
where we have used the fact that G has a pole at ω = ϵ, where G−1 vanishes.3 One can, therefore, obtain the correlated orbital energy ϵ directly from the self-energy evaluated at ω = ϵ and not as a small difference between two large quantities.3,22,25 Notice the recursive structure of this equation in which the unknown ϵ appears both in the left- and right-hand sides. We will show, in a later section, how this structure can be exploited to avoid divergence of a perturbation series in the presence of strong correlation and to obtain multiple physically meaningful roots from a single Dyson equation.13 The self-energy is subject to a diagrammatic many-body perturbation expansion18 and hence the Dyson equation itself, eq 10, can be expressed diagrammatically as in Figure 3. From this figure, we now see that each single-lined edge in the Hugenholtz diagram of Figure 1 has a more profound physical meaning16−18 as the zeroth-order Green’s function {G0}pp than just a Wick’s contraction. In this view, the same MBPT(2) diagram is interpreted18 as the product of the algebraic interpretations of diagrammatic components (vertexes and edges) summed and integrated over the intermediate quantities (the edge indexes and frequencies ω):
∑∫ 0
β
dτ {−⟨ij || ab⟩}
(16)
i ,j,a,b
{−⟨ab || ij⟩}{−G0(τ )}aa{−G0( −τ )}ii {−G0(τ )}bb{−G0( −t )}jj =
1 4
∑
⟨ij || ab⟩⟨ab || ij⟩
i ,j,a,b
∫0
∞
dτ
(17)
exp{(ϵi +ϵj −ϵa −ϵb)τ }
where β is the inverse temperature and is equal to infinity at zero temperature. This last expression is known as Laplace-transform MP2.27 This has been introduced by Almlöf27 as a purely mathematical transformation with potential algorithmic advantages, but it can be seen that τ has the physical meaning of imaginary time from the MBGF derivation. A use of this transformation10−12 with the most dramatic and unexpected advantage will be discussed in a later section. Dispersion Interactions. In the imaginary-frequency domain,18,26 the zeroth-order Green’s function is given by 657
dx.doi.org/10.1021/jp410587b | J. Phys. Chem. A 2014, 118, 655−672
The Journal of Physical Chemistry A {G0(ω)}pq =
Feature Article
δpq iω − ϵp +μ
(18)
With this, the MBPT(2) diagram is interpreted algebraically as ⎛ 1 ⎞e − E0(2) = (− 1)l ⎜ ⎟ ⎝2⎠
∞
∑∫ −∞
i ,j ,a,b
dω1 2π
∞
∫−∞
d ω2 2π
(19)
∞
dω3 {−⟨ij || ab⟩}{−⟨ab || ij⟩} 2π {− G0(ω1 + ω2)}aa{− G0(ω2)}ii {− G0(− ω1 + ω3)}bb{− G0(ω3)}jj
Figure 4. Bubble diagram.
∫−∞
=
1 4
∞
∑∫ −∞
i ,j ,a,b
⟨ij || ab⟩⟨ab || ij⟩ dω1 2π (iω1 + ϵi −ϵa)(− iω1 + ϵj −ϵb)
by the Feynman−Goldstone diagrams is their ability to serve as a rigorous criterion for size consistency.16,31,32 By just looking at the topology (or more specifically, connectedness) of the diagrams that define a method, one can tell whether the method yields an extensive energy and is, therefore, a size-consistent “many-body” method applicable to condensed matter. In the authors’ view,31,32 the only valid size-consistency criterion that exists today is the one based on diagrams or the equivalent kvector counting method in crystal orbital theory (not the one based on fragmentation); size consistency3,33 is synonymous with size extensivity2 because there is only one valid criterion. There are different views, criteria, and nomenclature of size consistency put forth by others.2,34,35 Size consistency of MBPT(2) can be proven as follows. Each vertex in the diagram of Figure 1 represents a two-electron integral and is inversely proportional (V−1) to the volume (V) of the system (see below for an explanation). When applied to crystals, each edge (or its algebraic interpretation as the zerothorder Green’s function) is associated with a momentum (k vector) in addition to a frequency (ω). In the MBPT(2) diagram, there are four edges and thus four k vectors, but only three are linearly independent because of the momentum conservation law. An integration over each of the independent k vectors (which are intermediate quantities) gives rise to a factor of V since
(20)
By carrying out the integration over ω1 in the last expression, which can be done analytically using the same contour in the Appendix, we once again recover the canonical MBPT(2) energy expression of eq 6. However, the unintegrated expression of eq 20 is instructive in showing how MBPT(2) describes dispersion interactions. Suppose that the system consists of two well-separated subsystems and that i and a label orbitals in subsystem A and j and b label those in subsystem B. We can then approximate each antisymmetrized two-electron integral by a Coulomb energy between two dipoles: ⟨ij || ab⟩ ≈ ⟨ij|ab⟩ ≈
xiaAx jbB + yiaA yjbB − 2ziaAz jbB r3
(21)
zAia
where is the z-dipole moment (along the A−B axis) of the iaorbital product density centered at subsystem A and r is the distance between the two subsystems. With this approximation, the summation over i and a and that over j and b in eq 20 can be carried out before the integration over ω1, yielding the wellknown Casimir−Polder formula28−30 of the van der Waals C6 coefficient: (2) EAB
3 ≈− 6 πr
∫0
∞
A
∑= k
∫ dk
(25)
To summarize (see Figure 5), the MBPT(2) diagram has two V−1 vertexes and three V1 edges (k vectors) and thus is proportional
B
dω1α (iω1)α (iω1)
V (2π )3
(22)
(2) Here, E(2) AB denotes the portion of E0 (or more specifically, of the
left diagram in Figure 2) that describes the dispersion interaction between the two subsystems and αA is the imaginary-frequencydependent polarizability in the isotropic approximation α A(iω1) =
1 A {αxx(iω1) + αyyA(iω1) + αzzA(iω1)} 3
(23)
with the individual polarizabilities given in the sum-over-state formula ⎫ ⎧ |xiaA|2 |xiaA|2 ⎬ αxxA(iω1) = −∑ ⎨ + iω + ϵi −ϵa −iω1 + ϵi −ϵa ⎭ i,a ⎩ 1
Figure 5. Volume dependence of diagrammatic components. The k vector of one of the four edges (dashed) is determined by the other three with the momentum conservation law and does not contribute to the volume dependence of the diagram.
(24)
Note the parallels between this expression and eq 20. In other words, the so-called “bubble” diagram (Figure 4),17,18 obtained by cutting in half the left diagram in Figure 2, can be understood to represent the sum-over-state polarizability of a subsystem at an imaginary frequency, i.e., an instantaneous or induced dipole moment. The correlation energy in the long-range is, therefore, nothing but a sum of dispersion interactions between moieties. Size Consistency. One of the important conceptual advantages (in addition to the expediency of derivations) offered
to V−2·V3 = V. Therefore, its correlation energy is extensive and MBPT(2) is size consistent. This analysis can be generalized to prove the extensivity of an energy described by any connected diagram composed of extensive vertexes such as the Hamiltonian vertexes.31,32 This constitutes the extensive diagram theorem,31,32 which, in turn, is a part of the celebrated linked-diagram theorem of Goldstone.36 658
dx.doi.org/10.1021/jp410587b | J. Phys. Chem. A 2014, 118, 655−672
The Journal of Physical Chemistry A
Feature Article
diverges at Δk = 0 inverse-quadratically.37,38 As seen in Figure 6, at Δk = 0, virtual orbital ã and occupied orbital i ̃ have the same momentum k1 and become identical to each other according to eqs 26 and 28. Likewise, orbitals b̃ and j ̃ also become the same. Hence, this is the situation in which there are degenerate highest occupied and lowest unoccupied MOs, such as on the Fermi surface of a metal. Let us analyze the effect of the divergence of two-electron integrals on MBPT(2), again, taking a HEG as an example of metals. Following the convention in solid-state physics, we partition the Hamiltonian such that Ĥ 0 is the kinetic-energy operator. With this partitioning and the orbitals in eqs 26−29, we have
An analogous theorem for intensive quantities and another for sums of extensive and intensive quantities have been derived by one of the authors.31,32 Once a method is judged size consistent by the diagrammatic criterion based on an analysis of a periodic solid, size consistency becomes an intrinsic attribute of the method regardless of whether it is applied to atoms, molecules, periodic, or nonperiodic solids. Size Dependence of Two-Electron Integrals. Why does a two-electron integral scale as V−1? The easiest way to answer this is to evaluate one for the simplest orbitals, i.e., the planewaves, which are the orbitals of a homogeneous electron gas (HEG):16,37 1 φa(̃ r ) = exp{i(k1 + Δk) ·r } (26) V φb(̃ r ) =
1 exp{−i(k 2 + Δk) ·r } V
(27)
φ i (̃ r ) =
1 exp{ik1·r } V
(28)
φ j (̃ r ) =
1 exp{−ik 2·r } V
(29)
|k1 + Δk|2 |k + Δk|2 k2 k 2 , ϵb ̃ = 2 , ϵi ̃ = 1 , ϵ j ̃ = 2 2 2 2 2
ϵa ̃ =
(32)
and ⟨bã |̃ i j̃ ̃ ⟩ =
4π V |k1 + k 2 + Δk|2
(33)
Substituting these as well as ⟨ãb̃|ij̃ ⟩̃ of eq 31 into the spatialorbital expression of eq 8, we obtain
The factor of V−1/2 in each orbital is to normalize it in the whole volume V. A two-electron integral over these orbitals, ⟨ãb̃|ij̃ ⟩,̃ is depicted diagrammatically in Figure 6. It represents the physical
E0(2) = −
V 2π
4 7
∫ dΔk ∫k
Second-Order Many-Body Perturbation Theory: An Eternal Frontier So Hirata,* Xiao He,† Matthew R. Hermes, and Soohaeng Y. Willow‡ Department of Chemistry, University of Illinois at Urbana−Champaign, 600 South Mathews Avenue, Urbana, Illinois 61801, United States ABSTRACT: Second-order many-body perturbation theory [MBPT(2)] is the lowest-ranked member of a systematic series of approximations convergent at the exact solutions of the Schrödinger equations. It has served and continues to serve as the testing ground for new approximations, algorithms, and even theories. This article introduces this basic theory from a variety of viewpoints including the Rayleigh− Schrödinger perturbation theory, the many-body Green’s function theory based on the Dyson equation, and the related Feynman−Goldstone diagrams. It also explains the important properties of MBPT(2) such as size consistency, its ability to describe dispersion interactions, and divergence in metals. On this basis, this article surveys three major advances made recently by the authors to this theory. They are a finite-temperature extension of MBPT(2) and the resolution of the Kohn−Luttinger conundrum, a stochastic evaluation of the correlation and self-energies of MBPT(2) using the Monte Carlo integration of their Laplace-transformed expressions, and an extension to anharmonic vibrational zero-point energies and transition frequencies based on the Dyson equation.
■
INTRODUCTION Second-order many-body perturbation theory,1,2 MBPT(2), is the lowest-ranked member of a systematic series of approximations formally convergent at the exact solutions of the electronic Schrödinger equations. It is also the simplest ab initio molecular orbital (MO) theory3−5 that can describe covalent, ionic, hydrogen-bond, and dispersion interactions accurately and from first principles. MBPT(2), therefore, enjoys special importance in electronic structure theory and increasingly in other areas of computational chemistry such as vibrational structure theory,6−8 serving as the testing ground for new approximations, algorithms, and even theories. There are many facets to this simple and fundamental theory, which are reviewed in the introduction. In the subsequent sections, we describe three major advances we have recently made to MBPT(2), namely, a finite-temperature extension of MBPT(2) and the resolution of the Kohn−Luttinger conundrum,9 a stochastic evaluation of MBPT(2) correlation energies and self-energies,10−12 and anharmonic vibrational MBPT(2) based on the Dyson equation.13 Rayleigh−Schödinger Perturbation Theory. The easiest way to introduce MBPT(2) is to use the Rayleigh−Schrödinger perturbation theory (RSPT).3 We partition the Hamiltonian (Ĥ ) into the zeroth-order part (Ĥ 0) and the perturbation (V̂ ) and expand the exact energy and wave function as infinite sums of progressively smaller, higher-order perturbation corrections: (1) (2) Ψ0 = Φ(0) 0 + Φ0 + Φ0 + ...
E0 =
E0(0)
+
E0(1)
+
E0(2)
+ ...
(0) (0) Ĥ 0 Φ(0) n = En Φn
Substituting these in the Schrödinger equation (with Ĥ ) and equating terms of comparable sizes, we obtain the general expression for E(2) 0 , which reads E0(2) =
∑
̂ (0) (0) ̂ (0) ⟨Φ(0) 0 |V |Φn ⟩⟨Φn |V |Φ0 ⟩ E0(0) − En(0)
n>0
(4)
By applying this to ab initio MO theory with the Møller− Plesset (MP) partitioning1 of the exact electronic Hamiltonian, we obtain E0(2) =
1 4
occ. vir.
∑∑ i,j
a,b
⟨Φ0|V̂ |Φijab⟩⟨Φijab|V̂ |Φ0⟩ ϵi +ϵj −ϵa −ϵb
(5)
as the second-order correction to the energy of the Hartree− Fock (HF) or self-consistent-field (SCF) theory. The sum of E(0) 0 and E(1) 0 is equal to the HF energy. Here, i and j run over all occupied HF spin−orbitals, a and b over all virtual spin−orbitals, Φ0 denotes the HF wave function, Φab ij is the two-electronexcited Slater determinant thereof, and ϵp is the energy of the pth canonical spin−orbital. Evaluating the factors in the numerator with, say, the Slater−Condon rules,3 we obtain E0(2) =
(1)
1 4
occ. vir.
∑∑ i,j
a,b
⟨ij || ab⟩⟨ab || ij⟩ ϵi +ϵj −ϵa −ϵb
(6)
where ⟨pq∥rs⟩ is an antisymmetrized two-electron integral in physicists’ notation.3 This is the canonical, spin−orbital energy expression of MBPT(2), which is familiar to many quantum
(2)
The subscript labels a state, with zero being the ground state, the parenthesized superscript denotes the perturbation order, and (0) Φ(0) n and En are the solutions of the zeroth-order Schrödinger equation with Ĥ 0, i.e., © 2013 American Chemical Society
(3)
Received: October 25, 2013 Revised: December 6, 2013 Published: December 12, 2013 655
dx.doi.org/10.1021/jp410587b | J. Phys. Chem. A 2014, 118, 655−672
The Journal of Physical Chemistry A
Feature Article
Hence, the modern derivation of MBPT(2) is simply to draw the one diagram in Figure 1 and directly interpret it algebraically, according to the established rules found in the literature.2,3,5 This goes as follows: The upper and lower vertexes in Figure 1 are interpreted as the antisymmetrized two-electron integrals, ⟨ij∥ab⟩ and ⟨ab∥ij⟩, respectively. A resolvent, (ϵi + ϵj − ϵa − ϵb)−1, is attached to each pair of adjacent vertexes. The orbital energies in it (with the plus and minus signs for occupied and virtual orbitals) are associated with the edges found in between the pair. Summations must be taken over all internal edges (edges whose two ends are attached to vertexes). Finally, a factor of (−1)h+l(1/2)e must be multiplied, where h, l, and e are the numbers of occupied-orbital (hole) edges, loops, and pairs of equivalent edges, respectively (two edges are equivalent when they start and end with the same vertexes and have the same occupied/virtual attribute). In this diagram, h = 2, l = 2, and e = 2. This interpretation again leads to eq 6, as it should. Quantum chemists at the forefront of MBPT and CC developments can, in fact, conceive new theories and devise efficient algorithms entirely diagrammatically.15 Spin Integration. For a closed-shell molecule (with an even number of electrons), we can rewrite the above spin−orbital expression into spatial−orbital one. It then becomes
chemists. However, derived in this way using RSPT, it may be more suitably referred to as the second-order Møller−Plesset perturbation (MP2) theory.3 Many-Body Perturbation Theory. The reliance on the Slater−Condon rules and the first quantization becomes almost an insurmountable obstacle for treating anything beyond MBPT(2), including the related coupled-cluster (CC) theory.5,14,15 We must switch to more powerful mathematical tools furnished by quantum field theory (QFT).5,16−18 The use of these tools distinguishes MBPT(2) from MP2 or, generally, MBPT from RSPT. A derivation based on QFT can ensure size consistency, which is a crucial requirement of any many-body theory applicable to atoms, molecules, and solids as well as everything in between in size. There are three essential QFT tools, which are, in the increasing degree of sophistication and expediency, the second quantization,3 the normal-ordered second quantization and Wick’s theorem,5,16 and the Feynman−Goldstone diagrams.5,16,18 Using the normal-ordered second quantization, the task of evaluating the factors in the numerator of eq 5 is replaced by a systematic process of drawing all possible Wick’s contraction lines between normal-ordered creation and annihilation operators. This process can even be automated.19 For instance,
occ. vir.
E0(2) = 2 ∑ ∑ ∼ ∼ i ̃, j a∼ , b
̃ ab̃ |̃ i j̃ ̃ ⟩ ⟨i j̃ ̃ |ab̃ ⟩⟨ − ϵ i ̃+ϵ j ̃−ϵa ̃−ϵb ̃
occ. vir.
∑∑ ∼ ∼ i ̃, j a∼ , b
⟨i j̃ ̃ |bã ⟩̃ ⟨ab̃ |̃ i j̃ ̃ ⟩ ϵ i ̃+ϵ j ̃−ϵa ̃−ϵb ̃ (8)
where ⟨p̃q̃|r̃s⟩̃ is a two-electron integral (not antisymmetrized) over spatial MOs, which are distinguished from spin−orbitals by tildes. This is the programmable equation of MBPT(2) using the spin-restricted HF reference wave function.3 It can be obtained by a spin-integration of eq 6 or by direct algebraic interpretation of Goldstone diagrams3 shown in Figure 2: the left and right diagrams correspond to the first and second terms of eq 8, respectively.
This leads back to the same canonical MBPT(2) expression given by eq 6, as it should. A few exercises of drawing Wick’s contraction lines in MBPT or CC derivations tend to lead to the realization that the act of drawing a contraction line between two operators is equivalent to the act of connecting two vertexes with an edge (a line) in the corresponding Feynman−Goldstone diagram. 5,16 The MBPT(2) diagram in the Hugenholtz style is shown in Figure 1. The four contraction lines in each term of eq 7 correspond to the four edges in this diagram. Other than the fact that the single diagram can encompass all four patterns of contraction in eq 7, there is one-to-one mapping between the normal-ordered second quantized and diagrammatic derivations of MBPT(2).
Figure 2. Second-order Goldstone diagrams (type N).
Integral Transformation. This expression (eq 8) may appear to suggest O(n4) dependence of the operation cost of MBPT(2) on the number of orbitals (n). This is not the case. The true bottleneck of an MBPT(2) calculation is the transformation of two-electron integrals from the atomic orbital (AO) to MO basis,3 namely, ⟨i j̃ ̃ |ab̃ ⟩̃ =
̃
̃
̃
̃ b Cμi *Cνj *CκaC λ ⟨μν|κλ⟩
∑ μ,ν ,κ ,λ
=
∑ Cμi ̃*(∑ Cνj ̃ *(∑ Cκa(̃ ∑ Cλb⟨̃ μν|κλ⟩))) μ
ν
κ
(9)
λ
where ⟨μν|κλ⟩ is a two-electron integral over AOs and is an expansion coefficient of the μth AO in the p̃th spatial MO. As indicated by the parentheses, the summation can be performed stepwise, but the operation cost for generating all MO-integrals nevertheless increases as staggering O(n5). The memory cost is Cp̃μ
Figure 1. Second-order Hugenholtz diagram (type N). 656
dx.doi.org/10.1021/jp410587b | J. Phys. Chem. A 2014, 118, 655−672
The Journal of Physical Chemistry A
Feature Article
also large and can scale as O(n4) in the most primitive implementation. Neither the AO- nor MO-integrals can usually fit in a core memory, and therefore, they need to be stored on disks or across many processors’ memories in a parallel machine. Worse still, the AO- and MO-integrals have completely different data layouts: the former according to distances between spatially localized AOs and the latter based on symmetries of delocalized MOs; the data locality is hard to achieve. In short, this is one of the most nightmarish parts of an entire electronic structure program. In a later section, we will describe a whole new algorithm of MBPT(2) that eliminates this integral-transformation step altogether.10−12 Many-Body Green’s Function Theory. Another, closely related viewpoint to MBPT is provided by many-body Green’s function theory (MBGF).16−18,20−24 Here, the central equation is the Dyson equation. In the matrix form, it reads3,20,22 G(ω) = G0(ω) + G0(ω)Σ(ω)G(ω)
Figure 3. Dyson equation. In the diagonal approximation, the doublelined edge represents {G}pp and the single-lined edge {G0}pp.
⎛ 1 ⎞e −iE0(2) = ( −1)l ⎜ ⎟ ⎝2⎠ ∞
∫−∞
(10)
(11)
dω2 2π
dω3 {−i⟨ij || ab⟩}{−i⟨ab || ij⟩} 2π
(14)
{G0(τ )}pq = ∓δpq exp{−(ϵp −μ)τ }
(15)
where μ is the chemical potential, τ is the time difference, and the minus sign is taken if ϵp > μ and the plus sign otherwise. The origin of time is set at the lower vertex of Figure 1 and the upper vertex occurs at time τ > 0. In this picture,18 the MBPT(2) energy is the product of the algebraic interpretations of diagrammatic components summed and integrated over the intermediate quantities, which now include τ:
(12)
which is known as the inverse Dyson equation.3 Substituting an exact one-electron energy ϵ for ω and neglecting off-diagonal elements (i.e., invoking the diagonal approximation), we find ϵ−ϵp = {Σ(ϵ)}pp
∞
∫−∞
where l = 2 and e = 2. It may be noticed that the sum of the frequencies (multiplied by −1 when associated with a hole edge) is conserved at any vertical position of the diagram. The integrations over ω in the above expression can be carried out analytically, as shown in the Appendix. The result is, of course, the same MBPT(2) expression, eq 6. Laplace-Transform MP2. The derivation based on MBGF can also be done using zeroth-order Green’s functions in the realor imaginary-time domain or even in the imaginary-frequency domain.18,26 In the imaginary-time domain, their elements are defined by
This element has a pole at the energy of the pth HF orbital (with the plus sign taken when pth orbital is virtual and the minus sign when it is occupied) and δ is a positive infinitesimal that nevertheless satisfies ∞·δ = ∞.18 The exact and HF Green’s functions are related to each other by the Dyson equation through the self-energy denoted by Σ(ω). All the matrices are dependent on frequency ω. Multiplying G0−1 from the left and G−1 from the right on the Dyson equation, we obtain G0−1(ω) = G−1(ω) + Σ(ω)
dω1 2π
{iG0(ω3)}jj
δpq ω − ϵp ± iδ
i ,j,a,b
∫−∞
{iG0(ω1 + ω2)}aa{iG0(ω2)}ii {iG0( −ω1 + ω3)}bb
where G is the exact Green’s function, which has poles at exact one-electron energies (correlated ionization or electron-attachment energies) and G0 is the zeroth-order Green’s function (in our case, the HF Green’s function) given by its elements18 as {G0(ω)}pq =
∞
∑
(13)
⎛ 1 ⎞e −E(2) = lim ( −1)l ⎜ ⎟ ⎝2⎠ β →∞
where we have used the fact that G has a pole at ω = ϵ, where G−1 vanishes.3 One can, therefore, obtain the correlated orbital energy ϵ directly from the self-energy evaluated at ω = ϵ and not as a small difference between two large quantities.3,22,25 Notice the recursive structure of this equation in which the unknown ϵ appears both in the left- and right-hand sides. We will show, in a later section, how this structure can be exploited to avoid divergence of a perturbation series in the presence of strong correlation and to obtain multiple physically meaningful roots from a single Dyson equation.13 The self-energy is subject to a diagrammatic many-body perturbation expansion18 and hence the Dyson equation itself, eq 10, can be expressed diagrammatically as in Figure 3. From this figure, we now see that each single-lined edge in the Hugenholtz diagram of Figure 1 has a more profound physical meaning16−18 as the zeroth-order Green’s function {G0}pp than just a Wick’s contraction. In this view, the same MBPT(2) diagram is interpreted18 as the product of the algebraic interpretations of diagrammatic components (vertexes and edges) summed and integrated over the intermediate quantities (the edge indexes and frequencies ω):
∑∫ 0
β
dτ {−⟨ij || ab⟩}
(16)
i ,j,a,b
{−⟨ab || ij⟩}{−G0(τ )}aa{−G0( −τ )}ii {−G0(τ )}bb{−G0( −t )}jj =
1 4
∑
⟨ij || ab⟩⟨ab || ij⟩
i ,j,a,b
∫0
∞
dτ
(17)
exp{(ϵi +ϵj −ϵa −ϵb)τ }
where β is the inverse temperature and is equal to infinity at zero temperature. This last expression is known as Laplace-transform MP2.27 This has been introduced by Almlöf27 as a purely mathematical transformation with potential algorithmic advantages, but it can be seen that τ has the physical meaning of imaginary time from the MBGF derivation. A use of this transformation10−12 with the most dramatic and unexpected advantage will be discussed in a later section. Dispersion Interactions. In the imaginary-frequency domain,18,26 the zeroth-order Green’s function is given by 657
dx.doi.org/10.1021/jp410587b | J. Phys. Chem. A 2014, 118, 655−672
The Journal of Physical Chemistry A {G0(ω)}pq =
Feature Article
δpq iω − ϵp +μ
(18)
With this, the MBPT(2) diagram is interpreted algebraically as ⎛ 1 ⎞e − E0(2) = (− 1)l ⎜ ⎟ ⎝2⎠
∞
∑∫ −∞
i ,j ,a,b
dω1 2π
∞
∫−∞
d ω2 2π
(19)
∞
dω3 {−⟨ij || ab⟩}{−⟨ab || ij⟩} 2π {− G0(ω1 + ω2)}aa{− G0(ω2)}ii {− G0(− ω1 + ω3)}bb{− G0(ω3)}jj
Figure 4. Bubble diagram.
∫−∞
=
1 4
∞
∑∫ −∞
i ,j ,a,b
⟨ij || ab⟩⟨ab || ij⟩ dω1 2π (iω1 + ϵi −ϵa)(− iω1 + ϵj −ϵb)
by the Feynman−Goldstone diagrams is their ability to serve as a rigorous criterion for size consistency.16,31,32 By just looking at the topology (or more specifically, connectedness) of the diagrams that define a method, one can tell whether the method yields an extensive energy and is, therefore, a size-consistent “many-body” method applicable to condensed matter. In the authors’ view,31,32 the only valid size-consistency criterion that exists today is the one based on diagrams or the equivalent kvector counting method in crystal orbital theory (not the one based on fragmentation); size consistency3,33 is synonymous with size extensivity2 because there is only one valid criterion. There are different views, criteria, and nomenclature of size consistency put forth by others.2,34,35 Size consistency of MBPT(2) can be proven as follows. Each vertex in the diagram of Figure 1 represents a two-electron integral and is inversely proportional (V−1) to the volume (V) of the system (see below for an explanation). When applied to crystals, each edge (or its algebraic interpretation as the zerothorder Green’s function) is associated with a momentum (k vector) in addition to a frequency (ω). In the MBPT(2) diagram, there are four edges and thus four k vectors, but only three are linearly independent because of the momentum conservation law. An integration over each of the independent k vectors (which are intermediate quantities) gives rise to a factor of V since
(20)
By carrying out the integration over ω1 in the last expression, which can be done analytically using the same contour in the Appendix, we once again recover the canonical MBPT(2) energy expression of eq 6. However, the unintegrated expression of eq 20 is instructive in showing how MBPT(2) describes dispersion interactions. Suppose that the system consists of two well-separated subsystems and that i and a label orbitals in subsystem A and j and b label those in subsystem B. We can then approximate each antisymmetrized two-electron integral by a Coulomb energy between two dipoles: ⟨ij || ab⟩ ≈ ⟨ij|ab⟩ ≈
xiaAx jbB + yiaA yjbB − 2ziaAz jbB r3
(21)
zAia
where is the z-dipole moment (along the A−B axis) of the iaorbital product density centered at subsystem A and r is the distance between the two subsystems. With this approximation, the summation over i and a and that over j and b in eq 20 can be carried out before the integration over ω1, yielding the wellknown Casimir−Polder formula28−30 of the van der Waals C6 coefficient: (2) EAB
3 ≈− 6 πr
∫0
∞
A
∑= k
∫ dk
(25)
To summarize (see Figure 5), the MBPT(2) diagram has two V−1 vertexes and three V1 edges (k vectors) and thus is proportional
B
dω1α (iω1)α (iω1)
V (2π )3
(22)
(2) Here, E(2) AB denotes the portion of E0 (or more specifically, of the
left diagram in Figure 2) that describes the dispersion interaction between the two subsystems and αA is the imaginary-frequencydependent polarizability in the isotropic approximation α A(iω1) =
1 A {αxx(iω1) + αyyA(iω1) + αzzA(iω1)} 3
(23)
with the individual polarizabilities given in the sum-over-state formula ⎫ ⎧ |xiaA|2 |xiaA|2 ⎬ αxxA(iω1) = −∑ ⎨ + iω + ϵi −ϵa −iω1 + ϵi −ϵa ⎭ i,a ⎩ 1
Figure 5. Volume dependence of diagrammatic components. The k vector of one of the four edges (dashed) is determined by the other three with the momentum conservation law and does not contribute to the volume dependence of the diagram.
(24)
Note the parallels between this expression and eq 20. In other words, the so-called “bubble” diagram (Figure 4),17,18 obtained by cutting in half the left diagram in Figure 2, can be understood to represent the sum-over-state polarizability of a subsystem at an imaginary frequency, i.e., an instantaneous or induced dipole moment. The correlation energy in the long-range is, therefore, nothing but a sum of dispersion interactions between moieties. Size Consistency. One of the important conceptual advantages (in addition to the expediency of derivations) offered
to V−2·V3 = V. Therefore, its correlation energy is extensive and MBPT(2) is size consistent. This analysis can be generalized to prove the extensivity of an energy described by any connected diagram composed of extensive vertexes such as the Hamiltonian vertexes.31,32 This constitutes the extensive diagram theorem,31,32 which, in turn, is a part of the celebrated linked-diagram theorem of Goldstone.36 658
dx.doi.org/10.1021/jp410587b | J. Phys. Chem. A 2014, 118, 655−672
The Journal of Physical Chemistry A
Feature Article
diverges at Δk = 0 inverse-quadratically.37,38 As seen in Figure 6, at Δk = 0, virtual orbital ã and occupied orbital i ̃ have the same momentum k1 and become identical to each other according to eqs 26 and 28. Likewise, orbitals b̃ and j ̃ also become the same. Hence, this is the situation in which there are degenerate highest occupied and lowest unoccupied MOs, such as on the Fermi surface of a metal. Let us analyze the effect of the divergence of two-electron integrals on MBPT(2), again, taking a HEG as an example of metals. Following the convention in solid-state physics, we partition the Hamiltonian such that Ĥ 0 is the kinetic-energy operator. With this partitioning and the orbitals in eqs 26−29, we have
An analogous theorem for intensive quantities and another for sums of extensive and intensive quantities have been derived by one of the authors.31,32 Once a method is judged size consistent by the diagrammatic criterion based on an analysis of a periodic solid, size consistency becomes an intrinsic attribute of the method regardless of whether it is applied to atoms, molecules, periodic, or nonperiodic solids. Size Dependence of Two-Electron Integrals. Why does a two-electron integral scale as V−1? The easiest way to answer this is to evaluate one for the simplest orbitals, i.e., the planewaves, which are the orbitals of a homogeneous electron gas (HEG):16,37 1 φa(̃ r ) = exp{i(k1 + Δk) ·r } (26) V φb(̃ r ) =
1 exp{−i(k 2 + Δk) ·r } V
(27)
φ i (̃ r ) =
1 exp{ik1·r } V
(28)
φ j (̃ r ) =
1 exp{−ik 2·r } V
(29)
|k1 + Δk|2 |k + Δk|2 k2 k 2 , ϵb ̃ = 2 , ϵi ̃ = 1 , ϵ j ̃ = 2 2 2 2 2
ϵa ̃ =
(32)
and ⟨bã |̃ i j̃ ̃ ⟩ =
4π V |k1 + k 2 + Δk|2
(33)
Substituting these as well as ⟨ãb̃|ij̃ ⟩̃ of eq 31 into the spatialorbital expression of eq 8, we obtain
The factor of V−1/2 in each orbital is to normalize it in the whole volume V. A two-electron integral over these orbitals, ⟨ãb̃|ij̃ ⟩,̃ is depicted diagrammatically in Figure 6. It represents the physical
E0(2) = −
V 2π
4 7
∫ dΔk ∫k
