Bilinear Constraints upon the Correlation Contribution to the Electron

Subscriber access provided by UNIV OF SOUTHERN INDIANA

Quantum Electronic Structure

Bilinear Constraints Upon the Correlation Contribution to the Electron-Electron Repulsion Energy as a Functional of the One-Electron Reduced Density Matrix Jerzy Cioslowski, Zsuzsanna E. Mihalka, and Agnes Szabados J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.9b00443 • Publication Date (Web): 11 Jul 2019 Downloaded from pubs.acs.org on July 23, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

BILINEAR CONSTRAINTS UPON THE CORRELATION CONTRIBUTION TO THE ELECTRON-ELECTRON REPULSION ENERGY AS A FUNCTIONAL OF THE ONE-ELECTRON REDUCED DENSITY MATRIX

Jerzy Cioslowski

a)

´ Mihálka *, Zsuzsanna E.

b,c)

´ , and Agnes Szabados

b)

a)

Institute of Physics, University of Szczecin, Wielkopolska 15, 70-451 Szczecin, Poland b)

Laboratory of Theoretical Chemistry, Institute of Chemistry, Faculty of Science, ELTE Eötvös Loránd University, POB 32, H-1518 Budapest, Hungary c)

Hevesy György Ph.D. School of Chemistry, Faculty of Science, ELTE Eötvös Loránd University, POB 32, H-1518 Budapest, Hungary Abstract: Perturbative analysis of the functional U [n, ψ] that yields the correlation component U of the electron-electron repulsion energy in terms of the vectors ψ(1) and n of the natural spinorbitals and their occupation numbers (the 1-matrix functional) facilitates examination of the flaws inherent to the present implementations of the density matrix functional theory. Recognizing that the practical usefulness of any approximate 1-matrix functional hinges upon its capability of exactly reproducing the leading contribution to U at the limit of vanishing electron-electron interactions gives rise to asymptotic bilinear constraints for the (exact or model) 2cumulant 2 G that enters the expression for U . The asymptotic behavior of certain blocks of 2 G is found to be equally important. These identities, which are obtained for both the single-determinantal and a model multi-determinantal cases, take precedence over the linear constraints commonly enforced in the course of approximate construction of such functionals. This observation reveals the futility of designing sophisticated approximations tailored for the second-order contribution to 2 G while neglecting proper formulation of the respective first-order contribution that in the case of the so-called JKL-only functionals requires abandoning the JKdependence altogether. It has its repercussions not only for the functionals of the PNOF family but also for the expressions involving only the L-type two-electron repulsion integrals (in the guise of their exchange counterparts) that account only for the correlation effects due to electrons with antiparallel spins and are welldefined only for spin-unpolarized and high-spin systems (yielding vanishing U for the latter). *To whom all the correspondence should be addressed

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

I. INTRODUCTION Witnessing equally impressive gains in software and hardware performance, the many decades of sustained development of approximate approaches to the electron correlation problem have produced quantum-chemical methods of remarkable diversity. The spectrum of these methods spans the semiempirical formalisms that rely upon sets of parameters specific to each element of the periodic table, those that involve general parameters or heuristic terms in the relevant energy expressions, and those that, thanks to being based upon mathematically well-defined approximations to the electronic wavefunction, are devoid of any ad hoc assumptions or parameters other than the specifics of the basis sets employed. At present, the methods of the first category are rarely used and hardly subject to any attempts at improvement, whereas the approaches of the third category have been perfected over the years, leaving relatively little room for their further development. Consequently, most (though not all1 ) of the novel practical implementations of the electronic structure theory that have emerged in the recent decade are of somewhat less rigorous nature than their more conventional counterparts of the ”true ab initio” provenance. The vast majority of these new implementations derive from functionals for the electronic energy, the extent of the empiricism involved being directly related to the argument of a particular functional. Thus, computationally viable incarnations of the density functional theory (DFT),2 in which the electronic energy is expressed in terms of the electron density, are usually associated with extensive parameterization. In fact, the recent proliferation of approximate DFT approaches is rapidly acquiring an uncanny resemblance to the plethora of semiempirical formalisms developed and employed during the bygone era of early quantum chemistry. At the other end of the spectrum, one finds functionals of the two-electron reduced density matrix (a.k.a. the 2-matrix), in which the N -representability of functional’s argument is partially enforced by incorporation of a set of necessary conditions or by its (often quite heuristic) reconstruction.3 The empiricism of such approaches stems from the multitude of the possible choices for the conditions and/or reconstruction schemes that are defined in a less systematic fashion than e.g. the truncation levels of the conventional configuration interaction (CI) or coupled cluster formalisms. This arbitrariness hampers ordering the levels of theory into systematic hierarchies that allow for increasingly accurate predictions of electronic properties. In particular, this problem afflicts the density cumulant functional theory (DCFT)4–7 that appears to be the most promising approach of this type. In the DFT formalism, the electron density with readily enforceable N 2


Page 2 of 30

Page 3 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


representability enters the energy expression comprising five components, of which two (the Coulomb component of the electron-electron repulsion energy W and the energy of interaction with the external potential) are exact and three (the kinetic energy as well as the exchange and correlation components of W ) have to be approximated. In contrast, DCFT employs an exact functional of the two-electron density cumulant matrix (a.k.a. the 2-cumulant) derived from the 2-matrix for which the N -representability can be enforced only in an approximate fashion.8, 9 A reasonable trade-off between these two extremes is offered by the density matrix functional theory (DMFT), in which the one-electron reduced density matrix (a.k.a. the 1-matrix) serves as the argument of a functional with only one (relatively small) energy component (i.e. the correlation component U of W ) that requires approximate treatment.10–13 In practical implementations, the functional of the 1-matrix that yields U (referred to in the following simply as the 1-matrix functional) is formulated in terms of natural spinorbitals (NOs) ψ(1) ≡ {ψp (1)} and their occupation numbers n ≡ {np } as the minimum13, 14 U [n, ψ] ≡ U (n, g) 1 X X X IJ np nq (gpqpq − gpqqp ) = min Bpqrs grspq − CI∗ CJ C 2 pq pqrs IJ

(1)

with respect to the normalized vector C ≡ {CI } of the CI expansion coefficients subject to the constraint X CI∗ CJ AIJ . (2) pq = np δpq

∀ pq

IJ

IJ In the above equations, A ≡ {AIJ pq } and B ≡ {Bpqrs } are the systemindependent tensors of the one- and two-electron coupling coefficients pertinent to the Slater determinants built from the NOs, whereas the two-electron −1 repulsion integrals gpqrs = hψp (1)ψq (2)|r12 |ψr (1)ψs (2)i are the elements of the tensor g ≡ {gpqrs }. It should be emphasized that, since this formulation derives from the CI expansion that implies a pure state, the vector C exists only for the sets of occupation numbers that simultaneously satisfy the P trivial ensemble N -representability conditions (i.e. p np = N , where N is the number of electrons, and ∀p 0 ≤ np ≤ 1)15 and the generalized Pauli conditions.16–20 This narrowing of the domain of U [n, ψ] assures avoidance of the recently uncovered conceptual difficulties arising from the nonequivalence of the energy minimizations over ensemble and pure N -representable 1-matrices21 (which, however, are not of concern in the case of closed-shell systems that exhibit the time-reversal symmetry16 ).

3



Page 4 of 30

Although appearing deceptively simple (especially in comparison with its DFT counterpart that requires auxiliary functions and nonlinear integral transforms22 ), the minimization in the r.h.s. of Eq. (1) does not provide a path to facile construction of the exact 1-matrix functional. This is so because such minimization is equivalent to solution of the constrained eigenproblem X X X IJ CJ Bpqrs grspq − AIJ η (3) pq pq = 0 ,

∀ I

pqrs

J

pq

where the matrix η ≡ {ηpq } of Lagrange multipliers enforces satisfaction of Eq. (2), that cannot be handled with the standard techniques of linear algebra such as matrix inversion and diagonalization. Moreover, even if such a solution were found, it would simply constitute a two-step (first the constrained C, followed by ψ and n) variational procedure offering no computational advantage over the conventional one-step CI formalism. An alternative approach that opens the avenue to practical implementations of DMFT commences with the observation [which readily follows from inspection of Eq. (1)] that U (n, g) is first-degree homogeneous in g 14, 23 and thus (by virtue of Euler’s homogeneous function theorem) X 2 Gpqrs (n, g) grspq , (4) U (n, g) = pqrs

where 2

Gpqrs (n, g) =

∂U (n, g) ∂grspq

(5)

are the elements of the 2-cumulant 2 G(n, g). Consequently, approximate expressions for U (n, g) can be derived from model 2-cumulants with elements that conform at least to the linear constraints24

∀

Gpqrs (n, g) = 2 G∗rspq (n, g) ,

(6)

Gpqrs (n, g) = −2 Gpqsr (n, g) = 2 Gqpsr (n, g) = −2 Gqprs (n, g) ,

(7)

1 Gpqrq (n, g) = − (1 − np ) np δpr 2

(8)

2

pqrs

∀

pqrs

2

and

∀ pr

X q

2

(note that these constraints are occasionally augmented with inequalities stemming from the so-called D-, G-, and Q-conditions8, 9, 25 ). The promise of 4


Page 5 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


DMFT being a treatment of the electron correlation problem that combines predictive power with low computational cost stems from judicious construction of such matrices, those with sharply reduced numbers of nonzero elements being of particular interest.25 Regrettably, the recent abundance of approximate 1-matrix functionals developed along these lines has not translated into emergence of new quantumchemical methods competitive in the accuracy/cost ratio with their conventional counterparts. In fact, surveys employing benchmarks based upon fewelectron harmonium atoms have uncovered substandard performance of such functionals with none of the 14 approximate expressions for U (n, g) tested being capable of reproducing with reasonable accuracy the exact electronelectron repulsion energies within three distinct correlation regimes of electronic states with five different spin multiplicities.26, 27 At first glance, one could be tempted to attribute these failures to the neglect of all the elements of 2 G(n, g) with more than two indices pertaining to NOs with different spatial parts (the JK-only approximation28 ). However, the neglected elements are known to be usually small in comparison with those retained in the functionals in question29 that are in fact parameterizations of the doubly-occupied (a.k.a. pair-occupied or seniority-0) CI (DOCI), which is known to perform rather well (especially for the non-dynamical correlation).30–33 Consequently, one can state with confidence that the source of the deficiencies in such approximate expressions for U (n, g) lies elsewhere. Prompted by these observations, we have recently embarked upon a detailed investigation of the behavior of U [n, ψ] within the weak correlation regime. The results of this investigation, which thoroughly elucidates the reasons behind the observed failures of the present practical implementations of DMFT and provides clear pathways to their alleviation, are reported in this paper. II. THE WEAK-CORRELATION REGIME Consider a nondegenerate eigenfunction of a nonrelativistic Hamiltonian describing a system of N electrons. The partial derivatives of the corresponding eigenenergy, which is given by X X 1 2 E(h, g) = Γpq hqp + Γpqrs grspq , (9) pq

pqrs

where both the previously defined tensor g and the matrix h ≡ {hpq } of the ˆ one-electron integrals hpq = hφp (1)|h(1)|φ q (1)i involve arbitrary one-electron functions (spinorbitals) φ(1) ≡ {φp (1)}, yield the elements of the respective 5



Page 6 of 30

1- and 2-matrices 1 Γ ≡ {1 Γpq } and 2 Γ ≡ {2 Γpqrs } as Γpq ≡ 1 Γpq (h, g) =

∂E(h, g) ∂hqp

Γpqrs ≡ 2 Γpqrs (h, g) =

∂E(h, g) ∂grspq

1

(10)

and 2

.

(11)

These elements possess the well-known properties 1

Γpq = 1 Γ∗qp ,

2

Γpqrs = 2 Γ∗rspq ,

2

Γpqrs = −2 Γpqsr = 2 Γqpsr = −2 Γqprs , (12)

and X

2

Γpqrq =

q

N −1 1 Γpr 2

,

X

1

Γpp = N

,

(13)

p

whereas those of the 2-cumulant 2 G ≡ {2 Gpqrs } , defined as 2

Gpqrs ≡ 2 Gpqrs (h, g) = 2 Γpqrs −

1 2

1

Γpr 1 Γqs − 1 Γps 1 Γqr

,

(14)

satisfy the conditions24 2

Gpqrs = 2 G∗rspq

,

2

Gpqrs = −2 Gpqsr = 2 Gqpsr = −2 Gqprs

,

(15)

and X

2

Gpqrq =

q

1 2

1

Γ2 − 1 Γ

pr

.

(16)

The energy E(h, λ g) (where λ ≥ 0) that pertains to scaled electronelectron interactions has the power-series expansion E(h, λ g) =

∞ X

E (m) (h, g) λm

(17)

m=0

around λ = 0 provided the state in question remains nondegenerate for all λ > 0. Consequently (note that here and in the following, individual quantities employed in different contexts are distinguished by their arguments), 1

∞ ∂E(h, λg) X 1 (m) m Γpq (λ) ≡ Γpq (h, λg) = = Γpq λ ∂hqp m=0 1

6


(18)

Page 7 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and 2

2

Γpqrs (λ) ≡ Γpqrs (h, λg) = λ

−1

∞ ∂E(h, λg) X 2 (m) m = Γpqrs λ ∂grspq m=0

,

(19)

where 1 (m) Γpq

(m)

≡ 1 Γpq (h, g) =

∂E (m) (h, g) ∂hqp

(20)

and 2 (m) Γpqrs

(m)

≡ 2 Γpqrs (h, g) =

∂E (m+1) (h, g) ∂grspq

.

(21)

The analogous expansions for the elements of the respective 2-cumulant read 2

2

Gpqrs (λ) ≡ Gpqrs (h, λg) =

∞ X

2

m G(m) pqrs λ

,

(22)

m=0

where 2

(m)

2 G(m) pqrs ≡ Gpqrq (h, g) m 1 X 1 (m−m0 ) 1 (m0 ) 1 (m−m0 ) 1 (m0 ) (m) = 2 Γpqrs − Γpr Γqs − Γps Γqr 2 m0 =0

.

(23)

P g) 2 = Combining the identity λ ∂E(h,λ pqrs Γpqrs (h, λ g) (λ grspq ) [which fol∂λ lows from the Hellmann-Feynman and Euler’s homogeneous function theorems applied to the electron-electron repulsion energy component of E(h, λ g)] with Eqs. (17)-(21) yields the sum rules X (m) 1 Γpq hqp = (1 − m) E (m) (h, g) (24) pq

and X

2 (m−1) Γpqrs

grspq = m E (m) (h, g) .

(25)

pqrs

In the following, the spinorbitals φ(1) are assumed to be eigenfunctions ˆ of the core Hamiltonian h(1) with the corresponding eigenvalues ≡ {p }. Under such assumption, combining Eqs. (13), (24), and (25) produces the identities X (m) X (m) X (m−1) 1 2 2 Γpp p = 2 (N − 1)−1 Γpqpq p = (m−1 − 1) Γpqrs grspq (26) p

pq

pqrs

7



Page 8 of 30

that serve as order-by-order consistency checks for the derivations employed in the following sections of this paper. II.1. The single-determinantal case In the case of the wavefunction at λ = 0 being a single Slater determinant, G(0) vanishes and the zeroth-order contribution to the 1-matrix has the elements 2

1 (0) Γpq

= n(0) p δpq

,

(27) (0)

where the zeroth-order contributions n(0) ≡ {np } to the occupation numbers equal either 0 or 1. Application of Eqs. (20), (21), and (23) to the well-known expressions for E (1) (h, g) and E (2) (h, g)34, 35 yields 1 (1) Γpq

=

νp − νq Gpq p − q

(28)

and 2

G(1) pqrs =

1 νp νq ηr ηs − ηp ηq νr νs gpqrs 2 p + q − r − s

.

(29)

In these equations, νp = δn(0) and ηp = δn(0) are the indicator functions p ,1 p ,0 for the ν (core) and η (virtual) subsets of the spinorbitals {φp (1)}, whereas P gpqrs = gpqrs − gpqsr and Gpq = ν g i i piqi are the elements of the tensor of antisymmetrized two-electron repulsion integrals and its weighted partial contraction, respectively. Here and in the following, all the fractions with vanishing numerators are assumed to equal zero. Since the expressions for {1 Γ(m) } and {2 G(m) } result from the perturbation theory based upon the core Hamiltonian rather than the MøllerPlesset partitioning, it is convenient at this point to introduce a new vec˜ tor φ(1) of spinorbitals related to φ(1) through a unitary transformation ˜ φ(1) = exp(−λX + ) φ(1), where the antihermitian matrix X has the ele(1) ments Xpq = (νp − νq ) 1 Γqp . While zeroing 1 Γ(1) , this transformation affects ˜ (2) to the 1-matrix neither 1 Γ(0) nor 2 G(1) . The second-order contribution 1 Γ in the transformed basis turns out to have the elements (see Appendix A for

8


Page 9 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


details of derivation) 1 ˜ (2) Γpq

X η ν ν − ν η η G G q p i q p i pi iq q − i p − q i X νp − νq νi − νj Gij gjpiq + − − p q i j ij 1 ˆ † X ηq νp νi νj ηk − νq ηp ηi ηj νk gkpij gijkq + Ppq 2 q + k − i − j p − q ijk X 2 ˜ (1) 2 ˜ (1) + 2 (ηp ηq − νp νq ) Gipjk Gjkiq , † = Pˆpq

(30)

ijk † where the operator Pˆpq acts on two-index quantities according to the pre† ˆ scription Ppq Apq = Apq + A∗qp . As the division of {φp (1)} into the ν and η subsets carries over to the ˜ (0) + λ 1 Γ ˜ (1) + λ2 1 Γ ˜ (2) comprises the rotated spinorbitals {φ˜p (1)}, the sum 1 Γ ηη, ην, and νη blocks whose all elements are proportional to λ2 , and the νν block that equals a unit matrix perturbed by such elements. Consequently, by virtue of Eqs. (27) and (30) combined with simple arguments based upon ˜ [2] (λ) with the elements perturbation theory, diagonalization of the matrix 1 Γ X 1 ˜ [2] 2 ˜ (1) 2 ˜ (1) Gipjk Gjkiq (31) Γpq (λ) = νp δpq + 2 λ2 (ηp ηq − νp νq ) ijk

(note the absence of the ”off-diagonal” blocks νη and ην) produces the eige˜ nvalues of the 1-matrix in the basis of φ(1) that are correct through the second order in λ. This connection between the second-order contribution to the 1-matrix and the first-order contribution to the 2-cumulant, which also follows from the irreducible Brillouin conditions that are obtained from rather convoluted derivations,5, 36 has two important properties. First, it conserves the number of particles and preserves the spin-blocked structure (if any) of the 1-matrix. Second, it is invariant to any unitary transformation of the spinorbitals that retains the distinction between the ν and η subsets. One ˜ such transformation, ψ(1) = Q φ(1), which relates the vector ψ(1) of the ˜ ˜ turns Eq. (31) into natural spinorbitals (NOs) to φ(1) and diagonalizes 1 Γ, the identity X 2 ¯ (1) 2 ¯ (1) n(2) Gipjk Gjkiq (32) p δpq = 2 (ηp ηq − νp νq ) ijk (2)

for the elements of the second-order contribution n(2) ≡ {np } that enters 9



Page 10 of 30

the power-series expansion n=n

(0)

+

∞ X

n(m) λm

(33)

m=2 (1)

¯ , which for the vector n ≡ n(λ) of the occupation numbers. The tensor 2 G ˜ (1) through the appropriate transformation involving Q, is obtains from 2 G readily identified as the leading (in this case, the first-order) term in the power-series expansion for the 2-cumulant in the NO basis. Retaining the structure of 2 G(1) given by Eq. (29), this term has six unique [i.e. unrelated by the permutational symmetries spelled out in Eq. (15)] blocks, namely νννν, νηνν, ννηη, νηνη, ηηνη and ηηηη, out of which only the third one comprises nonvanishing elements. Since (1 − np ) np = (np − νp ) (ηp − np ) =

np − νp − (np − νp )2 ηp − νp

,

(34)

where (np − νp )2 is of the order of λ4 rather than λ2 , it turns out that, for any pair (p, q) of indices pointing to the same subset of NOs, Eq. (32) gives rise to the bilinear constraint with manifest particle-hole symmetry X 1 1 2 Gipjk (n, g) 2 Gjkiq (n, g) = δpq (1 − np ) np ijk 2

(35)

that has to be satisfied at the limit of all the occupation numbers approaching either 0 or 1 by the elements of any exact or model 2-cumulant that enter Eqs. (4) and (5) for the correlation component of the electron-electron repulsion energy expressed in terms of n and g. In addition, all the elements { (1−n1p ) np 2 Gpqrs (n, g)} other than those pertaining to the pairs (p, q) and (r, s) pointing to different subsets of NOs have to remain finite at this limit. II.2. An illustrative example: two-electron systems Before proceeding further, it is instructive to discuss briefly the constraint (35) in the context of two-electron systems that are described by the nonrelativistic wavefunction37 X Ψ(1, 2) = Cp ψp (1) ψp• (2) (36) p

involving pairs { ψp (1), ψp• (1) } of NOs that for singlet states share the spatial components while differing in the spin functions, the converse being 10


Page 11 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


P 2 true for high-spin triplet states. Note that ∀p Cp• = −Cp , p |Cp | = 1, • • • • ∀p (p ) = p, ∀pq p 6= q ⇒ p 6= q , and the sum that enters Eq. (36) [in which the standard notation of i ≡ (~ri , σi ) has been used for the combined spatial and spin coordinates] encompasses both p and p• . The 1- and 2-matrices that correspond to this wavefunction have the elements 1

Γpq = 2 |Cp |2 δpq = np δpq

(37)

and 2

Γpqrs = Cr∗ Cp δp• q δr• s

.

(38)

The resulting matrix elements of the 2-cumulant, which read [compare Eq. (14)] 2

Gpqrs = Cr∗ Cp δp• q δr• s − 2 |Cp |2 |Cq |2 (δpr δqs − δps δqr ) ,

(39)

are easily confirmed to satisfy the conditions (15) and (16). Inserting the above result into the triple-index sum that enters the constraint (35) produces X 1 2 Gipjk 2 Gjkiq (1 − np ) np ijk X i 1h np 2 n2k δpq = 1 + 3 np + np − 2− 2 1 − np k

.

(40)

At the limit of vanishing electron-electron interactions, n0 = n0• → 1 and • } np = np• → 0. It is readily shown that, in conformance with the ∀p∈{0,0 / constraint (35), the r.h.s. of Eq. (40) indeed tends to 21 δpq at this limit. II.3. The multi-determinantal case The driving force behind the development of new electronic structure methods is the quest for an approach capable of handling both the dynamical and non-dynamical correlation effects with equal ease and accuracy. Depending explicitly on the occupation numbers, even approximate 1-matrix functionals hold the potential of fulfilling this goal. Therefore, derivation of bilinear constraints analogous to Eq. (35) for systems with spinorbitals that are degenerate in the absence of the electron-electron interactions is of much importance. Consider the N -electron wavefunction at λ = 0 given by Ψ(0; A) =

1 X ˆ+ Apq ξp ξq (φˆ+ q φp Ψ0 ) , 2 pq 11


(41)


Page 12 of 30

where A ≡ {Apq } is a skew-symmetric (a.k.a antisymmetric) matrix with the normalization 12 Tr A† A = 1, {φˆ+ p } are the creation operators, ξp = 1 − (νp + ηp ) is the indicator function for the ξ (active) subset of the spinorbitals {φp (1)} whose elements share the same eigenvalue ¯ of the core Hamiltonian, and Ψ0 is the Slater determinant built from the spinorbitals belonging to the ν (core) subset. The simplest case, which is considered here, corresponds to two P P electrons on four degenerate spinorbitals, i.e. p ξp = 4 and p νp = N −2. In this instance, the Youla decomposition38 of A (conveniently written as a 4×4 block) reads   0 a1 0 0 −a1 0 0 0  UT , A=U (42)  0 0 0 a2  0 0 −a2 0 where U is a unitary matrix, giving rise to the ξξ block of 1 Γ(0) equal to  2  |a1 | 0 0 0  0 |a1 |2 0 0   UT . A† A = U ∗  (43)  0 0 |a2 |2 0  0 0 0 |a2 |2 Thus, requesting the entire zeroth-order contribution 1 Γ(0) to the 1-matrix to be diagonal implies U being a unit matrix (or more generally a permutation matrix, which is of no relevance here), the occupation numbers of the spinorbitals {φp (1)} belonging to the subset ξ coming in pairs: two equaling |a1 |2 and the other two equaling |a2 |2 = 1 − |a1 |2 (which is a vivid example of the generalized Pauli condition). Unlike in the single-determinantal case, the zeroth-order contribution 2 G(0) to the 2-cumulant does not vanish entirely, its elements being instead given by 1 (0) 2 (0) Apq A∗rs − n(0) n (δ δ − δ δ ) ξp ξq ξr ξs . (44) Gpqrs = pr qs ps qr p q 2 ˜ [2] after Its first-order counterpart 2 G(1) , which is obtained together with 1 Γ rather tedious algebra (see Appendix B for details), consists of 34 = 81 blocks, out of which 21 are unique. Its nonvanishing elements comprise just eight unique blocks, namely ννηη, ννηξ, ννξξ, νξηη, νξηξ, νξξξ, ηηξξ, and ηξξξ. These elements turn out to conform to the analogue of the constraint (35) that reads X 2 1 1 (0) (0) (0) 2 Gjkiq (n, g) = δpq Ω n(0) p , ni , nj , nk ) Gipjk (n, g (1 − np ) np ijk 2 (45) 12


Page 13 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


for any pair (p, q) of indices pointing to either the ν or η subset of NOs. The weight function Ω np , nq , nr , ns ) = − (np + nq − 1) (nr + ns − 1) + |(np − nq ) (nr − ns )| −1 + δ(np +nq −1) (nr +ns −1)+|(np −nq ) (nr −ns )|,0 [(np + nq + nr + ns − 2)2 − 1] (46) . that enters the triple sum of Eq. (45) exhibits the particle hole symmetry [i.e. Ω np , nq , nr , ns ) = Ω 1 − np , 1 − nq , 1 − nr , 1 − ns )] and is invariant to the permutations p ↔ q, r ↔ s, and (p, q) ↔ (r, s). III. BILINEAR CONSTRAINTS FOR APPROXIMATE 1-MATRIX FUNCTIONALS In general, it is desirable that any new approach to the electron correlation problem improves upon the accuracy/cost ratio of the available formalisms. In particular, it means that at the limit of vanishing electron-electron interactions it produces correlation energies that coincide with those furnished by the second-order perturbation theory with the Møller-Plesset partitioning [MP2 a.k.a. MBPT(2)] while surpassing it in accuracy at finite correlation strengths. In the case of DMFT, this expectation can be fulfilled only when the bilinear constraints presented in this paper are satisfied by the relevant 1-matrix functionals. With a single exception,39 the 1-matrix functionals published thus far25, 26 involve gross approximations introduced in the course of construction of the model 2-cumulant. First of all, they overlook the nonequivalence of the firstdegree homogeneity and linearity of U (n, g) in g, the linearity [or, in other words, the independence of G(n, g) of g] being exhibited only in particular cases, e.g. the two-electron systems37 and the systems described by wavefunctions obtained within the APSG approximation.28 The fact that several of the known approximate 1-matrix functionals have been obtained by modification of their two-electron archetype29, 40 explains the ubiquity of this unfortunate oversimplification that is readily invalidated by inspection of the variational construction codified in Eqs. (1) and (2). It is also worth noting that the nonlinear dependence of U (n, g) on g is a prerequisite for a proper description of dispersion interactions.41 The second approximation, introduced for the sake of computational efficiency, is the neglect of all the elements of the 2-cumulant with more than two independent indices. In other words, the elements of G(n, g) are assumed to have the form • ◦ Gpqrs (n, g) = γpr (n) δp• q δr• s + γpq (n) (δp◦ r δq◦ s − δp◦ s δq◦ r ) ,

13


(47)


Page 14 of 30

where p• and p◦ are the indices of the NOs that are two companions of the NO with the index p. The pairing schemes that assign these companion NOs have the properties ∀p (p• )• = (p◦ )◦ = p, ∀pq p 6= q ⇒ p• 6= q • , and ∀pq p 6= q ⇒ p◦ 6= q ◦ , i.e. there are one-to one correspondences between p and p• , and ◦ • (n)} (n)} and γ ◦ (n) ≡ {γpq between p and p◦ . The matrices γ • (n) ≡ {γpq have the properties [compare Eqs. (6) and (7)] • • • • • ∗ (n) = −γpq γpq • (n) = −γp• q (n) = γp• q • (n) = [γqp (n)]

(48)

◦ ◦ γpq (n) = γqp (n) = [γp◦◦ q◦ (n)]∗

(49)

and .

The model 2-cumulant of the general form (47) is not completely defined until the relationships between the indices p, p• , and p◦ are specified. The choice ˆ ψp• (1) = made in the published 1-matrix functionals, namely ∀p ψp◦ (1) = S ˆ is the spin-flip operator, gives rise to the approximate expresψp (1), where S sion X • ◦ U (n, g) = γpq (n) Lpq + γpq (n) (Jpq − Kpq ) (50) pq

for U in terms of the matrices J ≡ {Jpq }, K ≡ {Kpq }, and L ≡ {Lpq } whose elements are the Coulomb Jpq = gpqpq , exchange Kpq = gpqqp , and timeinversion/exchange Lpq = gppqq integrals (note that L = K for real-valued NOs). The ”JKL-only” approximation (50), which is ill-defined for systems that are neither spin-unpolarized nor high-spin,26, 42 can be expected to satisfy the minimum requirement for the 1-matrix functionals spelled out at the opening paragraph of this section, i.e. the reproduction of the exact second-order contribution to U , only when the elements of the individual blocks of the firstorder contribution to the model 2-cumulant (47) either vanish or conform to the aforederived bilinear constraints. In particular, the vanishing of the νννν, νηνν, νηνη, ηηνη and ηηηη unique blocks in the single-determinantal case implies that the leading contribution to γ ◦ (n) equals zero. In other words, as far as its leading contribution is concerned, the JKL-only approximation is in reality the L-only approximation. Consequently, since the (p, p• ) pair of indices refers to NOs with opposite spins, the 1-matrix functional (50) does not account at all for the second-order contribution to the correlation component of the energy of repulsion between same-spin electrons, which means that the leading contribution to U it produces vanishes for high-spin systems. For spin-unpolarized systems, inserting the remaining part of the 14


Page 15 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


model 2-cumulant (47) into Eq. (35) furnishes the asymptotics (note the disappearance of the Kronecker delta) X 1 1 • |γip (n)|2 = (51) (1 − np ) np i 2 that constraints the leading contribution to the elements of γ • (n) with indices pointing to different subsets of NOs, the contribution to the other elements equaling zero. These conclusions carry over to the multi-determinantal case considered in the subsection II.3 of this paper, the only difference being the asymptotic constraint [compare Eq. (45)] X 1 (0) (0) (0) • 2 Ω ni , ni , n(0) p , np ) |γip (n)| (1 − np ) np i = −

• X |γip (n)|2 1 1 = (0) (0) (1 − np ) np i 2ni − 1 2np − 1 2

,

(52)

valid for any p that points to either the ν or η subset of NOs, replacing Eq. (51). As expected, the same contributions from the elements of the νη unique block of the matrix γ • (n) enter the sums that appear in Eqs. (51) and (52), whereas in the latter case those stemming from the νξ and ηξ unique blocks are enhanced by factors that are either smaller than −1 or greater than 1. IV. CONCLUSIONS Perturbative analysis of the functional U [n, ψ] that yields the correlation component U of the electron-electron repulsion energy in terms of the vectors ψ(1) and n of the natural spinorbitals and their occupation numbers facilitates examination of the flaws inherent to the present implementations of the DMFT approach. In turn, elucidating the sources of these flaws allows one to uncover the pathways to their remediation. As already mentioned in this paper, the practical usefulness of any approximate 1-matrix functional hinges upon its capability of exactly reproducing the leading contribution to U . For systems whose electronic wavefunctions become single Slater determinants at the limit of vanishing electron-electron interactions, such contribution arises from 2 G(1) . Consequently, the asymptotic bilinear constraint (35) takes precedence over its linear counterpart (8) that, despite being relevant only to the third-order contribution to U , is commonly enforced in the course of approximate reconstruction of the 2cumulant from which the 1-matrix functional follows.25, 26 In fact, the proper 15



linear constraint for 2 G(1) is simply X 2 (1) ∀ Gpqrq (n, g) = 0 . pr

Page 16 of 30

(53)

q

Equally important is assuring that all the elements of 2 G(1) that do not belong to the ννηη unique block vanish. These observations are slightly altered for systems that at the limit of vanishing electron-electron interactions are described by linear combinations of Slater determinants. In such instances, the zeroth-order contribution 2 G(0) to the 2-cumulant does not vanish, giving rise to the part of U that is due to the non-dynamical correlation effects. However, at least for the 4 orbitals/2 electrons active space considered in this paper, the elements of 2 G(1) remain subject to a bilinear constraint valid for the pairs of indices that point to the νν and ηη blocks, the only difference from the single-determinantal case being the appearance of the weighting function (note that, involving 2 G(0) and 2 G(2) in addition to 2 G(1) , the analogous bilinear constraint for the ξξ block is of no interest in the present context). Again, all the elements of 2 G(1) that do not belong to certain unique blocks have to vanish. The particular form of the constraint (35) and its generalization (45) leads to further observations concerning 2 G(1) and hence the leading contribution to the part of U that is due to the dynamical correlation effects. First of all, because of their bilinear nature, these constrains are not expected to alleviate the phase dilemma that plagues the construction of approximate 1-matrix functionals.43 Second, the occupation numbers can appear in such functionals only as the combinations {(1 − np ) np } that are compatible with the particle-hole symmetry. The results of the present investigation exert heavy toll on the existing approximate 1-matrix functionals.25, 26 Most importantly, they reveal the futility of constructing sophisticated approximations tailored for 2 G(2) while neglecting the proper formulation of 2 G(1) that in the case of the JKL-only functionals requires abandoning the JK-dependence altogether. While primarily impacting the functionals of the PNOF family,25 this conclusion has also repercussions for the expressions involving only the L-type two-electron repulsion integrals (in the guise of their exchange counterparts). First of all, such expressions turn out to account only for the correlation effects due to electrons with antiparallel spins and are well-defined only for spin-unpolarized and high-spin systems42 (note, however, that they have been widely employed for general open-shell systems as well26, 27, 44, 45 because of the incorrectly assumed equivalence between the {gpp• qq• } and {gpqqp } integrals). Consequently, they yield vanishing U for the latter systems (a failure that is shared 16


Page 17 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in various degree by the PNOF functionals26 ). Second, the exponents other than one-half in the power-functional45, 46 turn out to be precluded (as confirmed by the asymptotic behavior observed in benchmark calculations26, 27 ), reducing it to the M¨ uller-Buijse-Baerends expression47 1 X √ np nq − np nq Lpq , U (n, g) = (54) 2 pq which in the single-determinantal case has the leading contribution consistent with [compare Eq. (50)] q q 1 • (1 − np ) np ηp νq + (1 − nq ) nq νp ηq (55) γpq (n) = − 2 that, although exhibiting both proper block structure and explicit dependence on {(1 − np ) np }, obviously violates the bilinear constraint (51). Even worse, employment of Eq. (54) in conjunction with the multi-determinantal case considered in this paper gives rise to spurious blocks in the both the zeroth- and first-order contributions to the 2-cumulant . In light of these findings, any attempts at further development of 1-matrix functionals linear in J , K, and L are slated to be met with much skepticism. However, there is another general class of approximate expressions for U (n, g) that departs from the false assumption of linearity in g. These expressions are obtained by employing Eq. (4) in conjunction with 2-cumulants whose elements are implicitly defined as gpqrs 2 , (56) Gpqrs (n, g) = D(e, n) where the pseudoenergies e are determined by solving the bilinear constraints X 1 Ω np , ni , nj , nk ) 2 Gipjk (n, g 2 Gjkiq (n, g) = F (1 − np ) np δpq (57) 2 ijk inspired by Eq. (45). The functions D(e) (first-degree homogeneous in e), Ω np , nq , nr , ns ), and F(x) (with the property limx→0 F(x)/x = 1) are subject to parameterization within the confines of the proper behavior within the weak-correlation regime. Depending on the whole set of the two-electron repulsion integrals, such functionals are computationally more expensive than their JKL-only counterparts, but they do not suffer from any of their aforediscussed shortcomings. An attempt in this direction has been recently made39 but ultimately failed due to the divergence of the resulting functional at the complete-basis-set limit,48 underscoring the need for serious effort at development of this area of DMFT. In this context, the recently proposed OP-NOFT approach49 is also worth mentioning. 17



The final question to ponder here is whether the combination of Eqs. (45) and (46) is flexible enough to be valid for general active spaces, i.e. arbitrary numbers of electrons assigned to arbitrary numbers of NOs with degenerate eigenvalues of the core Hamiltonian. Should the answer to this question turn out to be affirmative, it would open an avenue to extensions of the DCFT approach equally applicable to single- and multi-determinantal systems. ACKNOWLEDGMENT The research described in this publication has been funded by the National Science Center (Poland) under the grant 2016/21/B/ST4/00597 and the National Research, Development and Innovation Office (NKFIH, Hun´ gary) under Grant No. K115744. In addition, support from the UNKP-18-3 New National Excellence Program of the Ministry of Human Capacities was ´ M. and support from the ELTE Institutional Excellence provided to Zs. E. Program (1783-3/2018/FEKUTSRAT) administered by the Hungarian Mi´ Sz. and Zs. E. ´ M.. The helpful nistry of Human Capacities was provided to A. comments of Prof. U˘gur Bozkaya (Hacettepe University) during the initial stages of this project are also acknowledged.

18


Page 18 of 30

Page 19 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


APPENDIX A: DERIVATIONS FOR THE SINGLE-DETERMINANTAL CASE Derivation of Eq. (30) commences with evaluation of 1 Γ(2) . The elements of this matrix are most conveniently obtained from partial derivatives of the expressions corresponding to the third-order energy diagrams already available in the literature.34, 35 For example, the second diagram on page 127 of the former reference (or the diagram D in Fig. 2 of the latter) gives rise to the energy contribution of ED =

X ijk

Gij Gjk Gki νi ηj ηk (i − j ) (i − k )

,

(A.1)

whose partial derivative with respect to Gqp , X Gpi Giq νi ηp ηq X νp η q 1 η p νq ∂ED Gpi Giq ηi , = + − ∂Gqp (p − i ) (q − i ) p − q i p − i q − i i (A.2) (2)

yields the respective contribution to 1 Γpq . Adding analogous contributions due to the diagrams labeled E-N35 produces X η ν ν − ν η η G G X ηp ηq νi − νp νq ηi q p i q p i pi iq 1 (2) † Γpq = Gpi Giq + Pˆpq ( − ) ( − ) − − p i q i q i p q i i X νp − νq νi − νj + Gij gjpiq p − q i − j ij 1 † X ηq νp νi νj ηk − νq ηp ηi ηj νk gkpij gijkq + Pˆpq 2 q + k − i − j p − q ijk X 2 (1) 2 (1) Gipjk Gjkiq (A.3) + 2 (ηp ηq − νp νq ) ijk

upon combining with Eq. (29). The orbital rotation removes the first term ˜ (1) = 2 G(1) ]. in the r.h.s. of Eq. (A.3), producing Eq. (30) [note that 2 G The detailed argument behind the assertion that the diagonalization of ˜ [2] (λ) with the elements given by Eq. (31) affords eigenvalues the matrix 1 Γ ˜ (0) +λ 1 Γ ˜ (1) + correct through the second order in λ runs as follows: The sum 1 Γ ˜ (2) differs from 1 Γ ˜ [2] (λ) by the ”off-diagonal” blocks νη and ην that λ2 1 Γ 2 are proportional to λ and arise from the first three terms in the r.h.s. of ˜ [2] (λ) as the Eq. (30). Treating this difference as the perturbation and 1 Γ unperturbed matrix leads to the conclusion that the leading correction to the eigenvalues (which arises in the second order of the perturbation theory) 19



Page 20 of 30

is proportional to λ4 (which comes from squares of the elements of the ”offdiagonal” blocks) divided by unity (the dominant term in the differences between the eigenvalues of the νν and ηη blocks), proving the assertion. Eq. (32) [which is simply Eq. (31) written in the basis of NOs] leads to the key identity (35) upon the observations that a) the difference between −νp is of the order of λ4 [compare Eq. (34)], b) np − νp (1 − np ) np and nηpp−ν p (2)

equals np λ2 , and c) ηp2 − νp2 is the same as ηp − νp . APPENDIX B: DERIVATIONS FOR THE MULTI-DETERMINANTAL CASE Derivation of Eq. (45) commences with the power-series expansion Ψ(λ) =

∞ X

Ψ(m) λm

(B.1)

m=0

ˆ + λW ˆ ˆ with scaled for the eigenfunction Ψ(λ) of the Hamiltonian H(λ) =h electron-electron interactions, where X ˆ= ˆ− h p φˆ+ (B.2) p φp p

and X ˆ+ ˆ− ˆ− ˆ =1 gpqrs φˆ+ W p φq φs φr 2 pqrs

.

(B.3)

(0) In the above equations, {φˆ− equals p } are the annihilation operators and Ψ Ψ(0; A) given by Eq. (41) with A resulting from Eq. (42) upon U set to a unit matrix. The elements of the first-order contributions to the respective 1- and 2-matrices read 1 (1) Γpq

(1) ˆ+ ˆ− (0) ˆ− (1) = hΨ(0) |φˆ+ q φp |Ψ i + hΨ |φq φp |Ψ i

(B.4)

and 2 (1) Γpqrs

=

1 1 ˆ+ ˆ− ˆ− (1) ˆ+ ˆ− ˆ− (0) hΨ(0) |φˆ+ hΨ(1) |φˆ+ r φs φq φp |Ψ i + r φs φq φp |Ψ i , (B.5) 2 2

whereas the elements of the second-order contributions to the respective 1matrix are given by 1 (2) Γpq

(2) ˆ+ ˆ− (0) ˆ− (2) = hΨ(0) |φˆ+ q φp |Ψ i + hΨ |φq φp |Ψ i (0) ˆ+ ˆ− (0) (1) (1) ˆ− (1) + hΨ(1) |φˆ+ q φp |Ψ i − hΨ |φq φp |Ψ i hΨ |Ψ i .

20


(B.6)

Page 21 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The Bloch formalism50 provides a convenient framework for computing Ψ and Ψ(2) . In the case considered here, it produces (1)

ˆW ˆ Ψ(0) Ψ(1) = Q

(B.7)

and ˆW ˆ Q ˆW ˆ −Q ˆ2 W ˆ O ˆW ˆ Ψ(0) Ψ(2) = Q

,

(B.8)

where the reduced resolvent ˆ= Q

X

|ΦK ihΦK |

K

ˆ (0) i − hΦK |h|Φ ˆ Ki hΨ(0) |h|Ψ

(B.9)

ˆ such involves the set {ΦK } of the single-determinantal eigenfunctions of h that ∀K ∀A hΨ(0; A)|ΦK i = 0, the orthogonality defining the projection opeˆ with the properties ∀A O ˆ Ψ(0; A) = Ψ(0; A) and ∀K O ˆ ΦK = 0 . The rator O denominator of Eq. (B.9) is given solely by the elements of the vector . Out of the eight terms in the expression 1 X gkpij + 2 Gpj δik ˆ+ Aij ξi ξj ξk ηp (φˆ+ p φk Ψ0 ) 2 ijkp ¯ − p gqpij 1 X ˆ+ Aij ξi ξj ηp ηq (φˆ+ + p φq Ψ0 ) 4 ijpq 2 ¯ − p − q gspir 1 X ˆ− ˆ+ ˆ+ + Aij ξi ξj ηp νr ηs (φˆ+ s φr φp φj Ψ0 ) 2 ijprs ¯ + r − p − s 1 X 2 gskjr − Gsr δjk ˆ− ˆ+ ˆ+ + Aij ξi ξj ξk νr ηs (φˆ+ s φr φi φk Ψ0 ) 2 ijkrs r − s 1 X gqpsr ˆ+ ˆ− ˆ− ˆ+ ˆ+ + Aij ξi ξj ηp ηq νr νs (φˆ+ p φq φs φr φj φi Ψ0 ) 4 ijpqrs r + s − p − q 1 X gkljr + Glr δjk ˆ− ˆ+ ˆ+ + Aij ξi ξj ξk ξl νr (φˆ+ i φr φl φk Ψ0 ) 2 ijklr r − ¯ 1 X gklsr ˆ+ ˆ− ˆ− ˆ+ ˆ+ + Aij ξi ξj ξk ξl νr νs (φˆ+ j φi φs φr φl φk Ψ0 ) 8 ijklrs r + s − 2 ¯ 1 X gkpsr ˆ+ ˆ− ˆ− ˆ+ ˆ+ + Aij ξi ξj ξk ηp νr νs (φˆ+ p φj φs φr φi φk Ψ0 ) 4 ijkprs r + s − p − ¯ Ψ(1) =

(B.10)

21



Page 22 of 30

that follows from Eqs. (B.2), (B.3), (B.7), and (B.9), only three contribute to 1 Γ(1) , yielding 1 (1) Γpq

Gpq [νq (ηp + ξp ) − νp (ηq + ξq )] q − p ∗ Bqp Bpq (νp + ηp ) ξq + (νq + ηq ) ξp + q − p p − q =

,

where the elements of the auxiliary matrices G and B are given by X (0) Gpq = ni gpiqi

(B.11)

(B.12)

i

and Bpq = n(0) q Gpq +

1 X Aij A∗kq gkpij ξi ξj ξk 2 ijk

.

(B.13)

An analogous treatment for 2 Γ(1) [where all the eight terms of Eq. (B.10) contribute] ultimately produces [compare Eq. (23)] ∗ ξp ξq νr νs 1 A∗rs θpq (2 ξp + ηp ) ηq ξr ξs + Apq θrs 2 (1) † − ˆ− ˆ ˆ Gpqrs = Ppqrs Ppq Prs 16 r + s − p − q ∗ 1 Ars κqp 1 Apq κ∗sr + ξp ηq ξr ξs + ξp ξq ξr νs 4 ¯ − q 4 ¯ − s 1 gpqrs + [2 n(0) r (ηp ηq + 2 ξp ηq + ξp ξq ) ξr νs 8 r + s − p − q (0) + (ηp ηq + ξp ξq − 2 n(0) p ξp ξq + 2 ηp ξq − 2 np ξp ηq ) νr νs ] 1 X giqjs ∗ + A Apj ξp (ηq + ξq ) ξr νs 2 ij s − q ri G −G 1 ˆ− ˆ− qs qs (0) [νs (ηq + ξq ) − νq (ηs + ξs )] − Ppq Prs δpr np ξp 2 s − q ∗ Bsq Bqs + ξs (νq + ηq ) + ξq (νs + ηs ) , s − q q − s (B.14) − † where the operators Pˆpq and Pˆpqrs act on two- and four-index quantities − † ˆ according to the prescriptions Ppq Apq = Apq − Aqp and Pˆpqrs Apqrs = Apqrs + ∗ Arspq , respectively, and the elements of the auxiliary matrices θ and κ read X θpq = Aij gpqij ξi ξj (B.15) ij

22


Page 23 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


κpq =

X

Aqi Gpi ξi

.

(B.16)

i

The elements of the first-order contribution to the 2-cumulant in the transformed basis are computed as 2 ˜ (1) Gpqrs

(1)

= 2 Gpqrs X (0) (0) (0) (0) + − Xip 2 Giqrs − Xiq 2 Gpirs + Xri 2 Gpqis + Xsi 2 Gpqri , (B.17) i

where the elements of the antihermitian matrix X are given by Xpq =

1 (1) Γqp (0) (0) np − nq

.

(B.18)

(0) (0) ˆ+ ˆ− 2 ˆ ˆ ˆ ˆ− ˆ ˆ ˆ ˆ (0) The brackets hΨ(0) |φˆ+ q φp Q W Q W |Ψ i and hΨ |φq φp Q W O W |Ψ i that arise upon inserting Eq. (B.8) into the expression (B.6) for 1 Γ(2) involve twenty-two and eight terms, respectively. Upon expansion into linear combinations of expectation values of products of creation and annihilation operators, tens of thousands primitive terms are generated, which are reduced with the help of algebraic manipulation software.51 The resulting elements of 1 (2) Γ enter the expression 1 X 1 (1) 1 (1) 1 1 1 ˜ (2) 1 (2) Γpq = Γpq + Γpi Γiq + (0) , (B.19) (0) (0) (0) 2 i np − ni nq − ni

for the elements of the second-order contribution to the 1-matrix in the transformed basis. In particular, for those of the νν and ηη blocks, one obtains, respectively 1 ˜ (2) Γpq

∗ ∗ X τpkij τqjik θpi θqi 1 X νi + ηi ξj ξk 4 i (p + i − 2 ¯) (q + i − 2 ¯) (p − i ) (q − i ) ijk gpijk gjkqi 1 X + 2 ijk (p + i − j − k ) (q + i − j − k ) (0) (0) (0) × [2 nk νi ξk − ni (ηk + ξk )] (ηj + ξj ) + ni ξi ξj ξk X gpjij gikqk (0) (0) + nj nk ηi ξj ξk (B.20) ( − ) ( − ) p i q i ijk

=−

23



Page 24 of 30

and 1 ˜ (2) Γpq

= +

∗ ∗ X τpkij τqjik θpi θqi 1 X ηi − νi ξj ξk 4 i (p + i − 2 ¯) (q + i − 2 ¯) ( p − i ) (q − i ) ijk 1 X gpijk gjkqi

2

ijk

(p + i − j − k ) (q + i − j − k ) (0)

(0)

× [2 nk (ηi + ξi ) ξk + (1 − ni ) νk ] νj X gpjij gikqk (0) (0) nj nk νi ξj ξk , − (p − i ) (q − i ) ijk where the elements of the auxiliary tensor τ read X τpqrs = Asi gpqri ξi .

(B.21)

(B.22)

i

˜ (2) are found to conform to the identity Upon inspection, these elements of 1 Γ X (0) 2 ˜ (1) 2 ˜ (1) (0) (0) 1 ˜ (2) , (B.23) Ω n(0) Γpq = 2 (ηp ηq − νp νq ) p , ni , nj , nk ) Gipjk Gjkiq ijk

from which the constraint (45) follows by virtue of the arguments presented in the subsection II.1 of this paper.

24


Page 25 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


References [1] See for example: Hättig, Ch.; Klopper, W.; Köhn, A.; Tew, D. P. Explicitly Correlated Electrons in Molecules. Chem. Rev. 2012, 112, 4-74; Mitroy, J.; Bubin, S.; Horiuchi, W.; Suzuki, Y.; Adamowicz, L.; Cencek, W.; Szalewicz, K.; Komasa, J.; Blume, D.; Varga, K. Theory and Application of Explicitly Correlated Gaussians. Rev. Mod. Phys. 2013, 85, 693-749, and the references cited therein. [2] Practically all of these approaches are based upon the Kohn-Sham formalism: Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133-A1138. [3] Mazziotti, D. A. Two-Electron Reduced Density Matrix as the Basic Variable in Many-Electron Quantum Chemistry and Physics. Chem. Rev. 2012, 112, 244-262. [4] Sokolov, A. Yu.; Schaefer, H. F.; Kutzelnigg, W. Density Cumulant Functional Theory from a Unitary Transformation: N-Representability, Three-Particle Correlation Effects, and Application to O4+ . J. Chem. Phys. 2014, 141, 074111. [5] Kutzelnigg, W. Density-Cumulant Functional Theory. J. Chem. Phys. 2006, 125, 171101. [6] Simmonett, A. C.; Wilke, J. J.; Schaefer, H. F.; Kutzelnigg, W. Density Cumulant Functional Theory: First Implementation and Benchmark Results for the DCFT-06 Model. J. Chem. Phys. 2010, 133, 174122. [7] Sokolov, A. Yu.; Simmonett, A. C.; Schaefer, H. F. Density Cumulant Functional Theory: The DC-12 Method, an Improved Description of the One-Particle Density Matrix. J. Chem. Phys. 2013, 138, 024107. [8] Kollmar, Ch. A Size Extensive Energy Functional Derived from a Double Configuration Interaction Approach: The Role of N Representability Conditions. J. Chem. Phys. 2006, 125, 084108. [9] DePrince III, A. E.; Mazziotti, D. A. Parametric Approach to Variational Two-Electron Reduced-Density-Matrix Theory. Phys. Rev. A 2007, 76, 042501. [10] Gilbert, T. L. Hohenberg-Kohn Theorem for Nonlocal External Potentials. Phys. Rev. B 1975, 12, 2111-2120.

25



[11] Donnelly R. A.; Parr, R. G. Elementary Properties of an Energy Functional of the First-Order Reduced Density Matrix. J. Chem. Phys. 1978, 69, 4431-4439. [12] Valone, S. M. Consequences of Extending 1-Matrix Energy Functionals from Pure–State Representable to all Ensemble Representable 1 Matrices. J. Chem. Phys. 1980, 73, 1344-1349. [13] Levy, M. Universal Variational Functionals of Electron Densities, FirstOrder Density Matrices, and Natural Spin-Orbitals and Solution of the v-Representability Problem. Proc. Natl. Acad. Sci. 1979, 76, 6062-6065. [14] Cioslowski, J.; Pernal, K.; Ziesche, P. Systematic Construction of Approximate One-Matrix Functionals for the Electron-Electron Repulsion Energy. J. Chem. Phys. 2002, 117, 9560-9566. [15] Coleman, A. J. Structure of Fermion Density Matrices. Rev. Mod. Phys. 1963, 35, 668-686. [16] Smith, D. W. N -Representability Problem for Fermion Density Matrices. II. The First-Order Density Matrix with N Even. Phys. Rev. 1966, 147, 896-898. [17] Borland, R. E.; Dennis, K. The Conditions on the One-Matrix for ThreeBody Fermion Wavefunctions with One-Rank Equal to Six. J. Phys. B: At. Mol. Phys. 1972, 5, 7-15. [18] Klyachko, A. Quantum Marginal Problem and N-Representability. J. Phys. Conf. Ser. 2006, 36, 72-86. [19] Schilling, C.; Gross, D.; Christandl, M. Pinning of Fermionic Occupation Numbers. Phys. Rev. Lett. 2013, 110, 040404. [20] Altunbulak, M.; Klyachko, A. The Pauli Principle Revisited. Commun. Math. Phys. 2008, 282, 287-322. [21] Schilling, C. Communication: Relating the Pure and Ensemble Density Matrix Functional. J. Chem. Phys. 2018, 149, 231102. [22] Cioslowski, J. Density Driven Self-Consistent Field Method. I. Derivation and Basic Properties. J. Chem. Phys. 1988, 89, 4871-4874. [23] Levy, M. in Density Matrices and Density Functionals; Erdahl, R., Smith, V. H., Jr., Eds.; Reidel: Dordrecht, 1987; p 479. 26


Page 26 of 30

Page 27 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


[24] Cioslowski, J. Many-Electron Densities and Reduced Density Matrices; Kluwer: New York, 2000. [25] Piris, M. A Natural Orbital Functional Based on an Explicit Approach of the Two-Electron Cumulant. Int. J. Quantum Chem. 2013, 113, 620630. [26] Cioslowski, J.; Piris, M.; Matito, E. Robust Validation of Approximate 1-Matrix Functionals with Few-Electron Harmonium Atoms. J. Chem. Phys. 2015, 143, 214101 and the references cited therein. [27] Cioslowski, J.; Strasburger, K. Five- and Six-Electron Harmonium Atoms: Highly Accurate Electronic Properties and Their Application to Benchmarking of Approximate 1-Matrix Functionals. J. Chem. Phys. 2018, 148, 144107. [28] Cioslowski, J.; Pernal, K.; Buchowiecki, M. Approximate 1-Matrix Functionals for the Electron-Electron Repulsion Energy from Geminal Theories. J. Chem. Phys. 2003, 119, 6443-6447. [29] Mentel, Ł. M.; van Meer, R.; Gritsenko, O. V.; Baerends, E. J. The Density Matrix Functional Approach to Electron Correlation: Dynamic and Nondynamic Correlation along the Full Dissociation Coordinate. J. Chem. Phys. 2014, 140, 214105. [30] Limacher, P. A.; Ayers, P. W.; Johnson, P. A.; De Baerdemacker, S.; Van Neck, D.; Bultinck, P. A New Mean-Field Method Suitable for Strongly Correlated Electrons: Computationally Facile Antisymmetric Products of Nonorthogonal Geminals. J. Chem. Theory Comp. 2013, 9, 1394-1401 and the references cited therein. [31] Bytautas, L.; Henderson, T. M.; Jiménez-Hoyos, C. A.; Ellis, J. K.; Scuseria, G. E. Seniority and Orbital Symmetry as Tools for Establishing a Full Configuration Interaction Hierarchy. J. Chem. Phys. 2011, 135, 044119. [32] Bytautas, L.; Scuseria, G. E.; Ruedenberg, K. Seniority Number Description of Potential Energy Surfaces: Symmetric Dissociation of Water, N2 , C2 , and Be2 . J. Chem. Phys. 2015, 143, 094105. [33] Alcoba, D. R.; Torre, A.; Lain, L.; Massaccesi, G. E.; O˜ na, O. B.; Ayers, P. W.; Van Raemdonck, M.; Bultinck, P.; Van Neck, D. Performance of Shannon-Entropy Compacted N-Electron Wave Functions for Configuration Interaction Methods. Theor. Chem. Acc. 2016, 135, 153. 27



[34] Shavitt, I.; Bartlett, R. J. Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory; Cambridge University Press: New York, 2009. [35] Wilson, S. Many-Body Perturbation Theory Using a Bare-Nucleus Reference Function: A Model Study. J. Phys. B: At. Mol. Phys. 1984, 17, 505-518; however, beware of the numerous errors in Eqs. (16)-(26). [36] Kutzelnigg, W.; Mukherjee, D. Irreducible Brillouin Conditions and Contracted Schrödinger Equations for N-Electron Systems. IV. Perturbative Analysis. J. Chem. Phys. 2004, 120, 7350-7368. [37] Löwdin, P.-O.; Shull, H. Natural Orbitals in the Quantum Theory of Two-Electron Systems. Phys. Rev. 1956, 101, 1730-1739. [38] Youla, D. C. A Normal Form for a Matrix under the Unitary Congruence Group. Can. J. Math. 1961, 13, 694-704. [39] Yasuda, K. Correlation Energy Functional in the Density-Matrix Functional Theory. Phys. Rev. A 2001, 63, 032517. [40] Gritsenko, O.; Pernal, K.; Baerends, E. J. An Improved Density Matrix Functional by Physically Motivated Repulsive Corrections. J. Chem. Phys. 2005, 122, 204102. Rohr, D. R.; Pernal, K.; Gritsenko, O. V.; Baerends, E. J. A Density Matrix Functional with Occupation Number Driven Treatment of Dynamical and Nondynamical Correlation. J. Chem. Phys. 2008, 129, 164105. [41] Cioslowski, J.; Pernal, K. Density Matrix Functional Theory of Weak Intermolecular Interactions. J. Chem. Phys. 2002, 116, 4802-4807. [42] Cioslowski, J.; Strasburger, K. (to be published). [43] Pernal, K.; Cioslowski, J. Phase Dilemma in Density Matrix Functional Theory. J. Chem. Phys. 2004, 120, 5987-5992. [44] Zarkadoula, E. N.; Sharma, S.; Dewhurst, J. K.; Gross, E. K. U.; Lathiotakis, N. N. Ionization Potentials and Electron Affinities from ReducedDensity-Matrix Functional Theory. Phys. Rev. A 2012, 85, 032504. Helbig, N.; Theodorakopoulos, G.; Lathiotakis, N. N. Fractional Spin in Reduced Density-Matrix Functional Theory. J. Chem. Phys. 2011, 135, 054109.

28


Page 28 of 30

Page 29 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Lathiotakis, N. N.; Marques, M. A. L. Benchmark Calculations for Reduced Density-Matrix Functional Theory. J. Chem. Phys. 2008, 128, 184103. Lathiotakis, N. N.; Helbig, N.; Gross, E. K. U. Open Shells in ReducedDensity-Matrix-Functional Theory. Phys. Rev. A 2005, 72, 030501. [45] Sharma, S.; Dewhurst, J. K.; Lathiotakis, N. N.; Gross, E. K. U. Reduced Density Matrix Functional for Many-Electron Systems. Phys. Rev. B 2008, 78, 201103. [46] Cioslowski, J.; Pernal, K. Constraints upon Natural Spin Orbital Functionals Imposed by Properties of a Homogeneous Electron Gas. J. Chem. Phys. 1999, 111, 3396-3400. Cioslowski, J.; Pernal, K. Description of a Homogeneous Electron Gas with Simple Functionals of the One-Particle Density Matrix. Phys. Rev. A 2000, 61, 034503. [47] M¨ uller, A. M. K. Explicit Approximate Relation Between Reduced Twoand One-Particle Density Matrices. Phys. Lett. A 1984, 105, 446-452. Buijse, M. A.; Baerends, E. J. An Approximate Exchange-Correlation Hole Density as a Functional of the Natural Orbitals. Mol. Phys. 2002, 100, 401-421. [48] Cioslowski, J.; Pernal, K. Variational Density Matrix Functional Theory Calculations with the Lowest-Order Yasuda Functional. J. Chem. Phys. 2002, 117, 67-71. [49] Gebauer, R.; Cohen, M. H.; Car, R. A Well-Scaling Natural Orbital Theory. Proc. Natl. Acad. Sci. 2016, 113, 12913-12918. ´ [50] Bloch, C. Sur la Théorie des Perturbations des Etats Liés. Nucl. Phys. 1958, 6, 329-347. [51] Mathematica, version 11.3; Wolfram Research: Champaign, IL 2018.

29



TOC 88x35mm (300 x 300 DPI)


Page 30 of 30

Bilinear Constraints upon the Correlation Contribution to the Electron

Recommend Documents