- Email: [email protected]
Adiabatic connections and properties of coupling-integrated exchange–correlation holes and pair densities in density functional theory
Chemical Physics Letters 480 (2009) 136–139
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Adiabatic connections and properties of coupling-integrated exchange–correlation holes and pair densities in density functional theory M.-E. Pistol a,*, C.-O. Almbladh b a b
Division of Solid State Physics/The Nanometer Structure Consortium, Lund University, P.O. Box 118, SE-221 00 Lund, Sweden Theoretical Physics, Lund University, P.O. Box 118, SE-221 00 Lund, Sweden
a r t i c l e
i n f o
Article history: Received 16 June 2009 In final form 20 August 2009 Available online 23 August 2009
a b s t r a c t We demonstrate a method that likely will improve the accuracy of functionals used in density functional theory and pair density functional theory. In deriving the method we show that the coupling-integrated exchange–correlation hole used in density functional theory obeys all the N-representability conditions as non-integrated holes and we give two such conditions, not normally used by the density functional theory community. We will also show that in addition to the normal adiabatic connection there exist a dual adiabatic connection where the kinetic energy is integrated over the inverse electron mass. Ó 2009 Elsevier B.V. All rights reserved.
1. Introduction It is fair to say that density functional theory (DFT) now plays a central role in computational many-body physics. DFT is related to pair density functional theory [1] and pair density matrix methods [2]. The basic variable in DFT is the electron density and a theorem of Hohenberg and Kohn says that the ground state energy of a many-body system is a unique functional of the electron density [3]. The energy due to the interaction with an external potential, e. g. a set of atomic nucleii is easy to find exactly. However, the part of the functional which is due to the kinetic energy of the electrons as well as the energy due to the interactions between the electrons is unknown. Kohn and Sham [4] showed that the DFT problem can be treated as a Hartree–Fock problem but where the kinetic energy was still unknown. Later it was shown, by using the adiabatic connection [5], that the kinetic energy can be treated ‘exactly’ and what is unknown can be lumped together into a functional of the exchange–correlation hole. The whole problem is thus reduced to finding out what the exchange–correlation hole looks like. The exchange–correlation hole is known to obey a set of exact conditions, known as the sum rule, the positivity rule and a ‘cusp condition’ relating the derivative of the on-top hole with the magnitude of the on-top hole [6,7]. Essentially all implementations of DFT use the adiabatic connection and it is very important to have as many constraints as possible on the exchange–correlation hole. We will in this Letter give new constraints on the exchange– correlation hole and we will also give two algorithms that will correct any exchange–correlation hole that violates the new con-
* Corresponding author. E-mail addresses: [email protected] (M.-E. Pistol), carl-olof.almblad@ teorfys.lu.se (C.-O. Almbladh). 0009-2614/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2009.08.047
straints. This will thus give a general method to ‘purify’ existing exchange–correlation functionals in DFT. The same method will also improve pair density functional theory, if the adiabatic connection is used. However pair density matrix methods are not addressed by our method since no adiabatic connection is used in this context. First, we will show that there exists a dual adiabatic connection in which the kinetic energy term (instead of the interaction energy) is integrated. We will then show that the exchange– correlation hole, after performing the normal adiabatic connection, obeys all the conditions obeyed by the exchange–correlation hole for an interacting system. These, so called N-representability conditions, are in principle known, but computationally hard. An N-representable exchange–correlation hole is a hole function that can be obtained from a wave-function and similarly for an N-representable pair density. Finally we will present our algorithms. Our results apply to bosons as well as to fermions. We use atomic units throughout the Letter and ‘N-representable’ always means ‘pure-state N-representable’. In DFT the energy of an interacting electron system is given by:
Z Z 1 ; x2 Þ 1 nðx1 Þnðx2 Þ 1 nðx1 Þhðx E½n; V ¼ T 0 ½n þ dx1 dx2 þ dx1 dx2 2 jx1 x2 j 2 jx1 x2 j Z þ nðx1 ÞVðx1 Þ dx1 ; ð1Þ where T 0 ½n is the kinetic energy of a set of non-interacting elecR 1 Þnðx2 Þ dx1 dx2 is the major trons giving the electron density nðx1 Þ; nðx jx1 x2 j R part of the Coulomb energy between the electrons, nðx1 ÞVðx1 Þ dx1 is the interaction energy of the electrons with an external potential and
Exc ¼
1 2
Z
1 ; x2 Þ nðx1 Þhðx dx1 dx2 jx1 x2 j
ð2Þ
137
M.-E. Pistol, C.-O. Almbladh / Chemical Physics Letters 480 (2009) 136–139
is the unknown exchange–correlation energy. Vðx1 Þ is the external 1 ; x2 Þ is the ‘coupling strength integrated’ expotential, and hðx change–correlation hole, which we will return to below. The exchange–correlation hole is defined as:
Cðx1 ; x2 Þ hðx1 ; x2 Þ ¼ nðx2 Þ 1 ; nðx1 Þnðx2 Þ
ð3Þ
where Cðx1 ; x2 Þ is the pair density, which in turn is defined by
Cðx1 ; x2 Þ ¼ NðN 1Þ
Z
jWðx1 ; x2 ; . . . ; xN Þj2 dx3 dxN
ð4Þ
in terms of the many-body wavefunction, Wðx1 ; x2 ; . . . ; xN Þ, and the number of electrons, N. A pair density, Cðx1 ; x2 Þ, is N-representable if there is a wavefunction such that Eq. (4) is true. It is easy to see that C is non-negative, symmetric and equal to zero if the arguments (where both space and spin are included) are equal. However, N-representability is a far more stringent condition. Most functions satisfying the above conditions are not N-representable. In order to make the Letter self-contained we will now describe the adiabatic connection (which has been described several times before, for instance in Ref. [5]). When deriving Eq. (1) one imagines a set of partially interacting systems, parametrised by a constant k, such that we have the following Hamiltonian, Hk ¼ T þ V k þ P k 1 i – j jxi xj j. The external potential V k is chosen such that the elec2 tron density does not depend on k. The energy of a partially interacting system is then given by:
Ek ¼ T k ½n þ
Z
nðxÞV k ðxÞ dx þ
k 2
Z
Ck ðx1 ; x2 Þ jx1 x2 j
dx1 dx2 ;
ð5Þ
Z Z k nðx1 Þnðx2 Þ Ek ¼ T k ½n þ nðxÞV k ðxÞ dx þ dx1 dx2 2 jx1 x2 j Z k nðx1 Þhk ðx1 ; x2 Þ þ dx1 dx2 ; 2 jx1 x2 j
ð6Þ
Z 1 Ck ðx1 ; x2 Þ dx1 dx2 2 jx1 x2 j Z Z 1 nðx1 Þnðx2 Þ ¼ nðxÞ dV k ðxÞ dx þ dk dx1 dx2 2 jx1 x2 j Z 1 nðx1 Þhk ðx1 ; x2 Þ þ dk dx1 dx2 ; 2 jx1 x2 j nðxÞ dV k ðxÞ dx þ dk
ð7Þ
where the kinetic term is zero due to the Hellman–Feynman theorem. Integrating over k we get:
E1 E0 ¼
Z
nðxÞðV 1 ðxÞ V 0 ðxÞÞ dx ZZ 1 Ck ðx1 ; x2 Þ þ dx1 dx2 dk 2 jx1 x2 j Z Z 1 nðx1 Þnðx2 Þ dx1 dx2 ¼ nðxÞðV 1 ðxÞ V 0 ðxÞÞ dx þ 2 jx1 x2 j ZZ 1 nðx1 Þhk ðx1 ; x2 Þ þ dx1 dx2 dk: 2 jx1 x2 j
Since E0 ¼ T 0 ½n þ
E1 ¼ T 0 ½n þ or
Z
1 2
Z
R
ð8Þ
ð9Þ
1 nðxÞV 1 ðxÞ dx þ 2
Cðx1 ; x2 Þ
El ¼ lT l ½n þ
jx1 x2 j
Z
dEl ¼ dlT l ½n þ
ð11Þ
nðxÞV l ðxÞ dx þ
Z
1 2
Cl ðx1 ; x2 Þ jx1 x2 j
dx1 dx2 :
ð12Þ
Z
nðxÞ dV l ðxÞ dx;
ð13Þ
where now the exchange–correlation term is zero due to the Hellman–Feynman theorem. Integrating over l we get:
E1 E0 ¼
Z
T l ½n dl þ
R
Z
nðxÞðV 1 ðxÞ V 0 ðxÞÞ dx:
nðxÞV 0 ðxÞ dx þ 12
Z
R
nðxÞV 1 ðxÞ dx þ
C0 ðx1 ;x2 Þ jx1 x2 j
1 2
Z
ð14Þ
dx1 dx2 we have:
C0 ðx1 ; x2 Þ jx1 x2 j
dx1 dx2 ;
ð15Þ
R where T½n ¼ T l ½n dl is the kinetic energy integrated over the inverse electron mass. It may be of concern that for l ¼ 0 we do not have a differential equation and the system is classical which has been investigated by Seidl et al. [9]. However we can always choose the starting point for the integration at l ¼ with being arbitrarily close to zero. Eq. (15) can of course also be written in terms of an exchange–correlation hole, instead of a pair density. Eq. (15) is quite similar to Eq. (5) with k ¼ 1. The difference is that in Eq. (5) we have little knowledge of C1 whereas C0 in Eq. (15) is known, even if the underlying wavefunction is not specified. By performing the dual adiabatic connection we can perform DFT with exact exchange–correlation, although it is necessary to have good approximations to T½n. It may now be contemplated that it may be possible to obtain two more dual pairs of adiabatic connections – between the external potential and the kinetic energy and between the external potential and the interaction energy, where the remaining term keeps the electron density constant. However since the energy due to the external potential only depends on the electron density, no extra information is obtained by such adiabatic connections. If we summarise, we have the following three expressions for the total energy of an electronic system:
E1 ¼ T 1 ½n þ
E1 ¼ T 0 ½n þ
ð10Þ
nðx1 Þnðx2 Þ dx1 dx2 jx1 x2 j
Again T l ½n means that the kinetic energy depends on l via the underlying wavefunction. We have the following relation:
E1 ¼ T½n þ dx1 dx2 ;
1 2
R ¼ R hk dk are the coupling-integrated pair where C ¼ Ck dk and h density and exchange–correlation hole, respectively. There are variations of this adiabatic connection where the interpolation between the non-interacting and the interacting system is not linear in k [8]. We will now show that there is a dual adiabatic connection, be and C. Again we imagine a fore we investigate the properties of h set of interacting systems, parametrised by a constant l, such that P 1 .l we have the following Hamiltonian, Hl ¼ lT þ V l þ 12 i – j jx x i jj is thus a parameter which turns on the kinetic energy and which can be seen as the inverse of the electron mass. The external potential V l is chosen such that the electron density does not depend on l. The energy of a system, as a function of l, is then given by:
nðxÞV 0 ðxÞ dx we have:
Z
nðxÞV 1 ðxÞ dx þ
Z
1 ; x2 Þ nðx1 Þhðx dx1 dx2 ; jx1 x2 j
E1 ¼ T½n þ
because of Eq. (3). Although the electron density does not depend on k we use the notation T k ½n to emphasize that the kinetic energy depends on k since the underlying wavefunction depends on k. We have the following two relations:
Z
þ
Since E0 ¼
or, equivalently:
dEk ¼
E1 ¼ T 0 ½n þ
Z
Z
Z
Z
nðxÞV 1 ðxÞ dx þ
1 2
nðxÞV 1 ðxÞ dx þ
1 2
nðxÞV 1 ðxÞ dx þ
1 2
Z
C1 ðx1 ; x2 Þ jx1 x2 j
Z
Z
Cðx1 ; x2 Þ jx1 x2 j
C0 ðx1 ; x2 Þ jx1 x2 j
dx1 dx2 ;
ð16Þ
dx1 dx2 ;
ð17Þ
dx1 dx2 :
ð18Þ
The question is now, what can we say about T and C? About T we are strongly constrained. cannot say much [10] but C, and thus h,
138
M.-E. Pistol, C.-O. Almbladh / Chemical Physics Letters 480 (2009) 136–139
It is known that pure-state N-representable pair densities belongs to a convex set, P½N [11–16]. P½N has been characterised completely, in terms of its boundary, which consists precisely of the classical pair configurations [11,12]. With a trivial modification of the method used in Refs. [11,12] is it possible to show that P½N; n, the set of N-representable pair densities that give a density nðxÞ, is also convex and characterised by its classical boundary. Unfortunately it is an NP-hard problem to determine membership in P½N, as shown by Deza and Laurent [17] and the problem is also quantum Merlin–Arthur complete [18]. The convexity of these sets means that any convex combination of N-representable pair densities is also N-representable. In other words,
Ci 2 P½N )
X
ai Ci 2 P½N;
ð19Þ
i
P
where ai P 0 and i ai ¼ 1. By inspecting Eq. (3) we see that if the electron density is fixed then the exchange–correlation hole belongs to a convex set as well which we denote P h ½N; n, which is independent of the coupling strength, k. When deriving the adiabatic connection we do indeed have a fixed electron density. Since can be the coupling-integrated exchange–correlation hole, h, approximated arbitrarily closely by a Riemann sum,
h
X
hki ðki ki1 Þ;
fully, imposing stronger N-representability conditions can improve the accuracy of DFT over LDA, for instance by decreasing the self-interaction energy [22] and by making overbinding harder to occur. When using the pair density as the basic variable in a DFT-like formalism [1,23,24] it is common to use the interacting system, i.e. Eq. (16), and for instance use the von Weisäcker approximation for the kinetic energy term [25]. One N-representability condition that is fairly strong is:
Z
e ðx1 ; x2 Þf ðx2 Þ dx1 P 0 f ðx1 Þ C
ð21Þ
for arbitrary f ðxÞ, where
e ðx1 ; x2 Þ ¼ Cðx1 ; x2 Þ þ dðx1 x2 Þnðx1 Þnðx2 Þ; C
ð22Þ
e is a positive definite operator, is a ‘modified’ pair density. That is, C which was shown by Garrod and Percus [26]. Another, related condition, is given by the fact that the Fourier transform of the distance density of a quantum many-body system is positive [27]. If we define the distance density by:
dðxÞ ¼
Z
Cðx1 ; x1 þ xÞ dx1 ;
ð23Þ
ð20Þ
i
2 P ½N; n. The which is a convex combination, we conclude that h h same argument applied to the pair density shows that C 2 P½N. This means that every property that must be true for hk ðCk Þ due to N ðCÞ. Some of these properties of representability is also true for h h have been derived previously [19] but not in this generality. The argument is unchanged for bosons, and also in spin-resolved DFT ðCÞ as for we have the same N-representability conditions for h hk ðCk Þ, independently of spin configuration. It is necessary to be careful when minimising the energy using variational calculus. The ground state energy for an interacting system can be found by minimizing the right hand side of Eqs. (5) or (16) over the set of C1 2 P½N or by minimising the right hand side of Eq. (6) over both nðxÞ and h1 2 P h ½N; n. The N-representability conditions on nðxÞ are very simple [20] and can be taken as n being integrable, positive and continuously differentiable except in a finite set of points (corresponding to the position of the atomic nucleii). The ‘only’ problem is that the kinetic energy term is not known but we emphasize that all other terms are known. By using Eqs. (10), (11) or (17) we can find the kinetic energy ‘exactly’ but we are no longer allowed to minimise the expressions over 2 P ½N; n. The problem is not variational in this case C 2 P½N or h h belong to P½N and P ½N; n, respecbut can be made so. C and h h tively. However there may be elements in P½N that give a lower energy than the correct C. By deleting these elements we get a variational problem again, assuming a non-degenerate ground state. Unfortunately C then becomes a member of an unknown subset of P½N which we may call S½N. The same reasoning can which then belongs to an unknown set be applied to h Sh ½N; n Ph ½N; n. Thus by performing the adiabatic connection we simplify the calculation of the kinetic energy, but this comes at the expense of getting unknown conditions on the pair density or the exchange–correlation hole. Similar problems arise if we use the dual adiabatic connection. It can be argued that nothing is gained by performing the adiabatic connection, either we have an unknown kinetic term where we know that the exchange–correlation hole belongs to Ph ½N; n or we have a known kinetic term where the exchange–correlation hole belongs to an unknown Sh ½N; n. However it appears that nature is kind, since the local density approximation (LDA) [21], using the adiabatic connection and very weak N-representability conditions on the exchange–correlation hole gives very good results in practical calculations. Hope-
then F½dðxÞ P 0 where F is the Fourier transform operator. Both of these conditions are true also in the spin-resolved case, i.e. for C assuming that the electrons have either parallel or antiparallel spin. They are true for bosons as well. An efficient characterisation of P½N is impossible at this time due to its NP-hardness [17]. However for some NP-hard problems there are efficient approximate algorithms [28]. It is also known that the N-representability conditions become more stringent the more particles we have, since P½M P½N if M > N [12,29], which is not reflected in the conditions given above. Studies of random pair densities have shown that there is structure in the spectra of certain matrices related to the pair densities and that this structure is N-dependent [30]. We find that the exchange–correlation hole introduced by Ernzerhof and Perdew [31] does not manifestly satisfy Eq. (22), and we suspect this can be true for many phenomenological exchange–correlation holes. e ic ðx1 ; x2 Þ, violates Eq. (21) there are If a modified pair density, C e ic ðx1 ; x2 Þ can be seen simple ways to correct it. First we note that C as the kernel of an integral operator. Thus it has eigenfunctions and eigenvalues, where some eigenvalues very likely are negative. However a positive definite operator has only non-negative eigene ic in a finite values which can be easily corrected. First we expand C set of its eigenfunctions:
e ic ðx1 ; x2 Þ ¼ C
X
ki fi ðx1 Þfi ðx2 Þ;
ð24Þ
i
where ki are the eigenvalues. We then replace negative eigenvalues e ic ðx1 ; x2 Þ to a positive definite operator which by zero, converting C then satisfies Eq. (21). e ic to a matrix, C e M;ic , by It is also possible to convert an incorrect C discretising the problem. Then diagonalise the matrix, forming e M;ic;diag will contain negative entries on the e M;ic;diag ¼ S1 C e M;ic S. C C diagonal which should be put to zero, making it a positive definite e M;c;diag . An inverse similarity transform will then create a operator, C e M;c;diag S1 which is close e M;c ¼ S C corrected modified pair density, C in Frobenius norm to the initial operator and which is positive definite. Naturally it will be necessary to correct the normalisation of e M;c . C We are presently implementing these algorithms in a DFT-program in order to check whether the algorithms will improve the accuracy of a given exchange–correlation potential. Since our method to improve DFT is so general, we expect that, at least for some functionals, this will be the case.
M.-E. Pistol, C.-O. Almbladh / Chemical Physics Letters 480 (2009) 136–139
Acknowledgements We thank C. Verdozzi for a critical reading of the manuscript. This work was performed within the nanometer structure consortium in Lund and supported by the Swedish Foundation for Strategic Research (SSF) and the Swedish Research Council (VR).
References [1] P. Ziesche, Phys. Lett. A 195 (3–4) (1994) 213. [2] E.R. Davidson, Reduced Density Matrices in Quantum Chemistry, Academic Press, London, 1976. [3] P. Hohenberg, W. Kohn, Phys. Rev. 136 (3B) (1964) B864. [4] W. Kohn, L.J. Sham, Phys. Rev. 140 (4A) (1965) 1133. [5] A.D. Becke, J. Chem. Phys. 88 (2) (1988) 1053, doi:10.1063/1.454274. [6] R.T. Pack, W.B. Brown, J. Chem. Phys. 45 (2) (1966) 556, doi:10.1063/ 1.1727605. [7] A.K. Rajagopal, J.C. Kimball, M. Banerjee, Phys. Rev. B 18 (5) (1978) 2339, doi:10.1103/PhysRevB.18.2339. [8] W. Yang, J. Chem. Phys. 109 (1) (1998) 10107. [9] M. Seidl, P. Gori-Giorgi, A. Savin, Phys. Rev. A 75 (4) (2007) 042511.
139
[10] L. Wittgenstein, Tractatus Logico-Philosophicus, Proposition 7, Routledge and Kegan Paul, 1922. [11] M.-E. Pistol, Chem. Phys. Lett. 400 (2004) 548. [12] M.-E. Pistol, Chem. Phys. Lett. 422 (2006) 363. [13] M.L. Yoseloff, Trans. Am. Math. Soc. 190 (1974) 1. [14] W.B. McRae, E.R. Davidson, J. Math. Phys. 13 (10) (1972) 1527. [15] P.W. Ayers, S.B. Liu, Phys. Rev. A 75 (2007) 022514. [16] P.W. Ayers, Phys. Rev. A 12 (2006) 042502. [17] M. Deza, M. Laurent, J. Comput. Appl. Math. 55 (1994) 217. [18] Y.-K. Liu, M. Christandl, F. Verstraete, Phys. Rev. Lett. 98 (11) (2007) 110503. [19] O. Gunnarsson, B.I. Lundqvist, Phys. Rev. B 13 (10) (1976) 4274. [20] T.L. Gilbert, Phys. Rev. B 12 (1975) 2111. [21] J. Perdew, Y. Wang, Phys. Rev. B 45 (23) (1992) 13244. [22] P. Mori-Sanchez, A.J. Cohen, W. Yang, J. Chem. Phys. 125 (2006) 201102. [23] M. Higuchi, K. Higuchi, Phys. Rev. A 75 (2007) 042510. [24] A. Nagy, J. Chem. Phys. 125 (18) (2006) 184104. [25] P. Ayers, J. Math. Phys. 46 (6) (2005) 062107. [26] C. Garrod, J.K. Percus, J. Math. Phys. 5 (1964) 1756. [27] M.-E. Pistol, Chem. Phys. Lett. 431 (2006) 216. [28] C.H. Papadimitriou, Computational Complexity, Addison-Wesley, Reading, Mill Valley, CA, 1994. [29] E. Davidson, Chem. Phys. Lett. 246 (3) (1995) 209. [30] M.-E. Pistol, Chem. Phys. Lett. 449 (1) (2007) 208. [31] M. Enzerhof, J.P. Perdew, J. Chem. Phys. 109 (9) (1998) 3313.
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Adiabatic connections and properties of coupling-integrated exchange–correlation holes and pair densities in density functional theory M.-E. Pistol a,*, C.-O. Almbladh b a b
Division of Solid State Physics/The Nanometer Structure Consortium, Lund University, P.O. Box 118, SE-221 00 Lund, Sweden Theoretical Physics, Lund University, P.O. Box 118, SE-221 00 Lund, Sweden
a r t i c l e
i n f o
Article history: Received 16 June 2009 In final form 20 August 2009 Available online 23 August 2009
a b s t r a c t We demonstrate a method that likely will improve the accuracy of functionals used in density functional theory and pair density functional theory. In deriving the method we show that the coupling-integrated exchange–correlation hole used in density functional theory obeys all the N-representability conditions as non-integrated holes and we give two such conditions, not normally used by the density functional theory community. We will also show that in addition to the normal adiabatic connection there exist a dual adiabatic connection where the kinetic energy is integrated over the inverse electron mass. Ó 2009 Elsevier B.V. All rights reserved.
1. Introduction It is fair to say that density functional theory (DFT) now plays a central role in computational many-body physics. DFT is related to pair density functional theory [1] and pair density matrix methods [2]. The basic variable in DFT is the electron density and a theorem of Hohenberg and Kohn says that the ground state energy of a many-body system is a unique functional of the electron density [3]. The energy due to the interaction with an external potential, e. g. a set of atomic nucleii is easy to find exactly. However, the part of the functional which is due to the kinetic energy of the electrons as well as the energy due to the interactions between the electrons is unknown. Kohn and Sham [4] showed that the DFT problem can be treated as a Hartree–Fock problem but where the kinetic energy was still unknown. Later it was shown, by using the adiabatic connection [5], that the kinetic energy can be treated ‘exactly’ and what is unknown can be lumped together into a functional of the exchange–correlation hole. The whole problem is thus reduced to finding out what the exchange–correlation hole looks like. The exchange–correlation hole is known to obey a set of exact conditions, known as the sum rule, the positivity rule and a ‘cusp condition’ relating the derivative of the on-top hole with the magnitude of the on-top hole [6,7]. Essentially all implementations of DFT use the adiabatic connection and it is very important to have as many constraints as possible on the exchange–correlation hole. We will in this Letter give new constraints on the exchange– correlation hole and we will also give two algorithms that will correct any exchange–correlation hole that violates the new con-
* Corresponding author. E-mail addresses: [email protected] (M.-E. Pistol), carl-olof.almblad@ teorfys.lu.se (C.-O. Almbladh). 0009-2614/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2009.08.047
straints. This will thus give a general method to ‘purify’ existing exchange–correlation functionals in DFT. The same method will also improve pair density functional theory, if the adiabatic connection is used. However pair density matrix methods are not addressed by our method since no adiabatic connection is used in this context. First, we will show that there exists a dual adiabatic connection in which the kinetic energy term (instead of the interaction energy) is integrated. We will then show that the exchange– correlation hole, after performing the normal adiabatic connection, obeys all the conditions obeyed by the exchange–correlation hole for an interacting system. These, so called N-representability conditions, are in principle known, but computationally hard. An N-representable exchange–correlation hole is a hole function that can be obtained from a wave-function and similarly for an N-representable pair density. Finally we will present our algorithms. Our results apply to bosons as well as to fermions. We use atomic units throughout the Letter and ‘N-representable’ always means ‘pure-state N-representable’. In DFT the energy of an interacting electron system is given by:
Z Z 1 ; x2 Þ 1 nðx1 Þnðx2 Þ 1 nðx1 Þhðx E½n; V ¼ T 0 ½n þ dx1 dx2 þ dx1 dx2 2 jx1 x2 j 2 jx1 x2 j Z þ nðx1 ÞVðx1 Þ dx1 ; ð1Þ where T 0 ½n is the kinetic energy of a set of non-interacting elecR 1 Þnðx2 Þ dx1 dx2 is the major trons giving the electron density nðx1 Þ; nðx jx1 x2 j R part of the Coulomb energy between the electrons, nðx1 ÞVðx1 Þ dx1 is the interaction energy of the electrons with an external potential and
Exc ¼
1 2
Z
1 ; x2 Þ nðx1 Þhðx dx1 dx2 jx1 x2 j
ð2Þ
137
M.-E. Pistol, C.-O. Almbladh / Chemical Physics Letters 480 (2009) 136–139
is the unknown exchange–correlation energy. Vðx1 Þ is the external 1 ; x2 Þ is the ‘coupling strength integrated’ expotential, and hðx change–correlation hole, which we will return to below. The exchange–correlation hole is defined as:
Cðx1 ; x2 Þ hðx1 ; x2 Þ ¼ nðx2 Þ 1 ; nðx1 Þnðx2 Þ
ð3Þ
where Cðx1 ; x2 Þ is the pair density, which in turn is defined by
Cðx1 ; x2 Þ ¼ NðN 1Þ
Z
jWðx1 ; x2 ; . . . ; xN Þj2 dx3 dxN
ð4Þ
in terms of the many-body wavefunction, Wðx1 ; x2 ; . . . ; xN Þ, and the number of electrons, N. A pair density, Cðx1 ; x2 Þ, is N-representable if there is a wavefunction such that Eq. (4) is true. It is easy to see that C is non-negative, symmetric and equal to zero if the arguments (where both space and spin are included) are equal. However, N-representability is a far more stringent condition. Most functions satisfying the above conditions are not N-representable. In order to make the Letter self-contained we will now describe the adiabatic connection (which has been described several times before, for instance in Ref. [5]). When deriving Eq. (1) one imagines a set of partially interacting systems, parametrised by a constant k, such that we have the following Hamiltonian, Hk ¼ T þ V k þ P k 1 i – j jxi xj j. The external potential V k is chosen such that the elec2 tron density does not depend on k. The energy of a partially interacting system is then given by:
Ek ¼ T k ½n þ
Z
nðxÞV k ðxÞ dx þ
k 2
Z
Ck ðx1 ; x2 Þ jx1 x2 j
dx1 dx2 ;
ð5Þ
Z Z k nðx1 Þnðx2 Þ Ek ¼ T k ½n þ nðxÞV k ðxÞ dx þ dx1 dx2 2 jx1 x2 j Z k nðx1 Þhk ðx1 ; x2 Þ þ dx1 dx2 ; 2 jx1 x2 j
ð6Þ
Z 1 Ck ðx1 ; x2 Þ dx1 dx2 2 jx1 x2 j Z Z 1 nðx1 Þnðx2 Þ ¼ nðxÞ dV k ðxÞ dx þ dk dx1 dx2 2 jx1 x2 j Z 1 nðx1 Þhk ðx1 ; x2 Þ þ dk dx1 dx2 ; 2 jx1 x2 j nðxÞ dV k ðxÞ dx þ dk
ð7Þ
where the kinetic term is zero due to the Hellman–Feynman theorem. Integrating over k we get:
E1 E0 ¼
Z
nðxÞðV 1 ðxÞ V 0 ðxÞÞ dx ZZ 1 Ck ðx1 ; x2 Þ þ dx1 dx2 dk 2 jx1 x2 j Z Z 1 nðx1 Þnðx2 Þ dx1 dx2 ¼ nðxÞðV 1 ðxÞ V 0 ðxÞÞ dx þ 2 jx1 x2 j ZZ 1 nðx1 Þhk ðx1 ; x2 Þ þ dx1 dx2 dk: 2 jx1 x2 j
Since E0 ¼ T 0 ½n þ
E1 ¼ T 0 ½n þ or
Z
1 2
Z
R
ð8Þ
ð9Þ
1 nðxÞV 1 ðxÞ dx þ 2
Cðx1 ; x2 Þ
El ¼ lT l ½n þ
jx1 x2 j
Z
dEl ¼ dlT l ½n þ
ð11Þ
nðxÞV l ðxÞ dx þ
Z
1 2
Cl ðx1 ; x2 Þ jx1 x2 j
dx1 dx2 :
ð12Þ
Z
nðxÞ dV l ðxÞ dx;
ð13Þ
where now the exchange–correlation term is zero due to the Hellman–Feynman theorem. Integrating over l we get:
E1 E0 ¼
Z
T l ½n dl þ
R
Z
nðxÞðV 1 ðxÞ V 0 ðxÞÞ dx:
nðxÞV 0 ðxÞ dx þ 12
Z
R
nðxÞV 1 ðxÞ dx þ
C0 ðx1 ;x2 Þ jx1 x2 j
1 2
Z
ð14Þ
dx1 dx2 we have:
C0 ðx1 ; x2 Þ jx1 x2 j
dx1 dx2 ;
ð15Þ
R where T½n ¼ T l ½n dl is the kinetic energy integrated over the inverse electron mass. It may be of concern that for l ¼ 0 we do not have a differential equation and the system is classical which has been investigated by Seidl et al. [9]. However we can always choose the starting point for the integration at l ¼ with being arbitrarily close to zero. Eq. (15) can of course also be written in terms of an exchange–correlation hole, instead of a pair density. Eq. (15) is quite similar to Eq. (5) with k ¼ 1. The difference is that in Eq. (5) we have little knowledge of C1 whereas C0 in Eq. (15) is known, even if the underlying wavefunction is not specified. By performing the dual adiabatic connection we can perform DFT with exact exchange–correlation, although it is necessary to have good approximations to T½n. It may now be contemplated that it may be possible to obtain two more dual pairs of adiabatic connections – between the external potential and the kinetic energy and between the external potential and the interaction energy, where the remaining term keeps the electron density constant. However since the energy due to the external potential only depends on the electron density, no extra information is obtained by such adiabatic connections. If we summarise, we have the following three expressions for the total energy of an electronic system:
E1 ¼ T 1 ½n þ
E1 ¼ T 0 ½n þ
ð10Þ
nðx1 Þnðx2 Þ dx1 dx2 jx1 x2 j
Again T l ½n means that the kinetic energy depends on l via the underlying wavefunction. We have the following relation:
E1 ¼ T½n þ dx1 dx2 ;
1 2
R ¼ R hk dk are the coupling-integrated pair where C ¼ Ck dk and h density and exchange–correlation hole, respectively. There are variations of this adiabatic connection where the interpolation between the non-interacting and the interacting system is not linear in k [8]. We will now show that there is a dual adiabatic connection, be and C. Again we imagine a fore we investigate the properties of h set of interacting systems, parametrised by a constant l, such that P 1 .l we have the following Hamiltonian, Hl ¼ lT þ V l þ 12 i – j jx x i jj is thus a parameter which turns on the kinetic energy and which can be seen as the inverse of the electron mass. The external potential V l is chosen such that the electron density does not depend on l. The energy of a system, as a function of l, is then given by:
nðxÞV 0 ðxÞ dx we have:
Z
nðxÞV 1 ðxÞ dx þ
Z
1 ; x2 Þ nðx1 Þhðx dx1 dx2 ; jx1 x2 j
E1 ¼ T½n þ
because of Eq. (3). Although the electron density does not depend on k we use the notation T k ½n to emphasize that the kinetic energy depends on k since the underlying wavefunction depends on k. We have the following two relations:
Z
þ
Since E0 ¼
or, equivalently:
dEk ¼
E1 ¼ T 0 ½n þ
Z
Z
Z
Z
nðxÞV 1 ðxÞ dx þ
1 2
nðxÞV 1 ðxÞ dx þ
1 2
nðxÞV 1 ðxÞ dx þ
1 2
Z
C1 ðx1 ; x2 Þ jx1 x2 j
Z
Z
Cðx1 ; x2 Þ jx1 x2 j
C0 ðx1 ; x2 Þ jx1 x2 j
dx1 dx2 ;
ð16Þ
dx1 dx2 ;
ð17Þ
dx1 dx2 :
ð18Þ
The question is now, what can we say about T and C? About T we are strongly constrained. cannot say much [10] but C, and thus h,
138
M.-E. Pistol, C.-O. Almbladh / Chemical Physics Letters 480 (2009) 136–139
It is known that pure-state N-representable pair densities belongs to a convex set, P½N [11–16]. P½N has been characterised completely, in terms of its boundary, which consists precisely of the classical pair configurations [11,12]. With a trivial modification of the method used in Refs. [11,12] is it possible to show that P½N; n, the set of N-representable pair densities that give a density nðxÞ, is also convex and characterised by its classical boundary. Unfortunately it is an NP-hard problem to determine membership in P½N, as shown by Deza and Laurent [17] and the problem is also quantum Merlin–Arthur complete [18]. The convexity of these sets means that any convex combination of N-representable pair densities is also N-representable. In other words,
Ci 2 P½N )
X
ai Ci 2 P½N;
ð19Þ
i
P
where ai P 0 and i ai ¼ 1. By inspecting Eq. (3) we see that if the electron density is fixed then the exchange–correlation hole belongs to a convex set as well which we denote P h ½N; n, which is independent of the coupling strength, k. When deriving the adiabatic connection we do indeed have a fixed electron density. Since can be the coupling-integrated exchange–correlation hole, h, approximated arbitrarily closely by a Riemann sum,
h
X
hki ðki ki1 Þ;
fully, imposing stronger N-representability conditions can improve the accuracy of DFT over LDA, for instance by decreasing the self-interaction energy [22] and by making overbinding harder to occur. When using the pair density as the basic variable in a DFT-like formalism [1,23,24] it is common to use the interacting system, i.e. Eq. (16), and for instance use the von Weisäcker approximation for the kinetic energy term [25]. One N-representability condition that is fairly strong is:
Z
e ðx1 ; x2 Þf ðx2 Þ dx1 P 0 f ðx1 Þ C
ð21Þ
for arbitrary f ðxÞ, where
e ðx1 ; x2 Þ ¼ Cðx1 ; x2 Þ þ dðx1 x2 Þnðx1 Þnðx2 Þ; C
ð22Þ
e is a positive definite operator, is a ‘modified’ pair density. That is, C which was shown by Garrod and Percus [26]. Another, related condition, is given by the fact that the Fourier transform of the distance density of a quantum many-body system is positive [27]. If we define the distance density by:
dðxÞ ¼
Z
Cðx1 ; x1 þ xÞ dx1 ;
ð23Þ
ð20Þ
i
2 P ½N; n. The which is a convex combination, we conclude that h h same argument applied to the pair density shows that C 2 P½N. This means that every property that must be true for hk ðCk Þ due to N ðCÞ. Some of these properties of representability is also true for h h have been derived previously [19] but not in this generality. The argument is unchanged for bosons, and also in spin-resolved DFT ðCÞ as for we have the same N-representability conditions for h hk ðCk Þ, independently of spin configuration. It is necessary to be careful when minimising the energy using variational calculus. The ground state energy for an interacting system can be found by minimizing the right hand side of Eqs. (5) or (16) over the set of C1 2 P½N or by minimising the right hand side of Eq. (6) over both nðxÞ and h1 2 P h ½N; n. The N-representability conditions on nðxÞ are very simple [20] and can be taken as n being integrable, positive and continuously differentiable except in a finite set of points (corresponding to the position of the atomic nucleii). The ‘only’ problem is that the kinetic energy term is not known but we emphasize that all other terms are known. By using Eqs. (10), (11) or (17) we can find the kinetic energy ‘exactly’ but we are no longer allowed to minimise the expressions over 2 P ½N; n. The problem is not variational in this case C 2 P½N or h h belong to P½N and P ½N; n, respecbut can be made so. C and h h tively. However there may be elements in P½N that give a lower energy than the correct C. By deleting these elements we get a variational problem again, assuming a non-degenerate ground state. Unfortunately C then becomes a member of an unknown subset of P½N which we may call S½N. The same reasoning can which then belongs to an unknown set be applied to h Sh ½N; n Ph ½N; n. Thus by performing the adiabatic connection we simplify the calculation of the kinetic energy, but this comes at the expense of getting unknown conditions on the pair density or the exchange–correlation hole. Similar problems arise if we use the dual adiabatic connection. It can be argued that nothing is gained by performing the adiabatic connection, either we have an unknown kinetic term where we know that the exchange–correlation hole belongs to Ph ½N; n or we have a known kinetic term where the exchange–correlation hole belongs to an unknown Sh ½N; n. However it appears that nature is kind, since the local density approximation (LDA) [21], using the adiabatic connection and very weak N-representability conditions on the exchange–correlation hole gives very good results in practical calculations. Hope-
then F½dðxÞ P 0 where F is the Fourier transform operator. Both of these conditions are true also in the spin-resolved case, i.e. for C assuming that the electrons have either parallel or antiparallel spin. They are true for bosons as well. An efficient characterisation of P½N is impossible at this time due to its NP-hardness [17]. However for some NP-hard problems there are efficient approximate algorithms [28]. It is also known that the N-representability conditions become more stringent the more particles we have, since P½M P½N if M > N [12,29], which is not reflected in the conditions given above. Studies of random pair densities have shown that there is structure in the spectra of certain matrices related to the pair densities and that this structure is N-dependent [30]. We find that the exchange–correlation hole introduced by Ernzerhof and Perdew [31] does not manifestly satisfy Eq. (22), and we suspect this can be true for many phenomenological exchange–correlation holes. e ic ðx1 ; x2 Þ, violates Eq. (21) there are If a modified pair density, C e ic ðx1 ; x2 Þ can be seen simple ways to correct it. First we note that C as the kernel of an integral operator. Thus it has eigenfunctions and eigenvalues, where some eigenvalues very likely are negative. However a positive definite operator has only non-negative eigene ic in a finite values which can be easily corrected. First we expand C set of its eigenfunctions:
e ic ðx1 ; x2 Þ ¼ C
X
ki fi ðx1 Þfi ðx2 Þ;
ð24Þ
i
where ki are the eigenvalues. We then replace negative eigenvalues e ic ðx1 ; x2 Þ to a positive definite operator which by zero, converting C then satisfies Eq. (21). e ic to a matrix, C e M;ic , by It is also possible to convert an incorrect C discretising the problem. Then diagonalise the matrix, forming e M;ic;diag will contain negative entries on the e M;ic;diag ¼ S1 C e M;ic S. C C diagonal which should be put to zero, making it a positive definite e M;c;diag . An inverse similarity transform will then create a operator, C e M;c;diag S1 which is close e M;c ¼ S C corrected modified pair density, C in Frobenius norm to the initial operator and which is positive definite. Naturally it will be necessary to correct the normalisation of e M;c . C We are presently implementing these algorithms in a DFT-program in order to check whether the algorithms will improve the accuracy of a given exchange–correlation potential. Since our method to improve DFT is so general, we expect that, at least for some functionals, this will be the case.
M.-E. Pistol, C.-O. Almbladh / Chemical Physics Letters 480 (2009) 136–139
Acknowledgements We thank C. Verdozzi for a critical reading of the manuscript. This work was performed within the nanometer structure consortium in Lund and supported by the Swedish Foundation for Strategic Research (SSF) and the Swedish Research Council (VR).
References [1] P. Ziesche, Phys. Lett. A 195 (3–4) (1994) 213. [2] E.R. Davidson, Reduced Density Matrices in Quantum Chemistry, Academic Press, London, 1976. [3] P. Hohenberg, W. Kohn, Phys. Rev. 136 (3B) (1964) B864. [4] W. Kohn, L.J. Sham, Phys. Rev. 140 (4A) (1965) 1133. [5] A.D. Becke, J. Chem. Phys. 88 (2) (1988) 1053, doi:10.1063/1.454274. [6] R.T. Pack, W.B. Brown, J. Chem. Phys. 45 (2) (1966) 556, doi:10.1063/ 1.1727605. [7] A.K. Rajagopal, J.C. Kimball, M. Banerjee, Phys. Rev. B 18 (5) (1978) 2339, doi:10.1103/PhysRevB.18.2339. [8] W. Yang, J. Chem. Phys. 109 (1) (1998) 10107. [9] M. Seidl, P. Gori-Giorgi, A. Savin, Phys. Rev. A 75 (4) (2007) 042511.
139
[10] L. Wittgenstein, Tractatus Logico-Philosophicus, Proposition 7, Routledge and Kegan Paul, 1922. [11] M.-E. Pistol, Chem. Phys. Lett. 400 (2004) 548. [12] M.-E. Pistol, Chem. Phys. Lett. 422 (2006) 363. [13] M.L. Yoseloff, Trans. Am. Math. Soc. 190 (1974) 1. [14] W.B. McRae, E.R. Davidson, J. Math. Phys. 13 (10) (1972) 1527. [15] P.W. Ayers, S.B. Liu, Phys. Rev. A 75 (2007) 022514. [16] P.W. Ayers, Phys. Rev. A 12 (2006) 042502. [17] M. Deza, M. Laurent, J. Comput. Appl. Math. 55 (1994) 217. [18] Y.-K. Liu, M. Christandl, F. Verstraete, Phys. Rev. Lett. 98 (11) (2007) 110503. [19] O. Gunnarsson, B.I. Lundqvist, Phys. Rev. B 13 (10) (1976) 4274. [20] T.L. Gilbert, Phys. Rev. B 12 (1975) 2111. [21] J. Perdew, Y. Wang, Phys. Rev. B 45 (23) (1992) 13244. [22] P. Mori-Sanchez, A.J. Cohen, W. Yang, J. Chem. Phys. 125 (2006) 201102. [23] M. Higuchi, K. Higuchi, Phys. Rev. A 75 (2007) 042510. [24] A. Nagy, J. Chem. Phys. 125 (18) (2006) 184104. [25] P. Ayers, J. Math. Phys. 46 (6) (2005) 062107. [26] C. Garrod, J.K. Percus, J. Math. Phys. 5 (1964) 1756. [27] M.-E. Pistol, Chem. Phys. Lett. 431 (2006) 216. [28] C.H. Papadimitriou, Computational Complexity, Addison-Wesley, Reading, Mill Valley, CA, 1994. [29] E. Davidson, Chem. Phys. Lett. 246 (3) (1995) 209. [30] M.-E. Pistol, Chem. Phys. Lett. 449 (1) (2007) 208. [31] M. Enzerhof, J.P. Perdew, J. Chem. Phys. 109 (9) (1998) 3313.