Modelling the Pauli Potential in the Pair Density Functional Theory

Modelling the Pauli Potential in the Pair Density Functional Theory C. Amovilliaand ´ A. Nagyb aDipartimento di Chimica e Chimica Industriale, Universit`a di Pisa, Via Risorgimento 35, 56126 Pisa, Italy bDepartment of Theoretical Physics, University of Debrecen, H–4010 Debrecen, Hungary October 22, 2008 Abstract In the ground state the pair density can be determined by solving a single auxiliary equation of a two-particle problem. A novel method for determining the Pauli potential entering this equation is presented and, starting from a reliable description of the pair density, an analytical expression is derived for atomic systems. Test calculations are presented for Be and isoelectronic C2+ and O4+ ions. 1 Introduction Generalized density functional theories have received a growing importance in recent years. For electron systems, the interest has been posed on the pair density as the fundamental variable instead of the one particle density. It turned out that there exist a variational principle for the pair density (analogous to the Hohenberg-Kohn theorems of the density functional theory). It has been shown that - instead of Kohn-Sham equations - in the pair density functional theory [1–4] the ground state 1 problem of an arbitrary system is reduced to a two-particle problem. The twoparticle equation is written [1–4] as −1 2∇2 1−1 2∇2 2+v(r1) + v(r2) + N−1 r+vP(r1,r2)χ(r1,r2) =µχ(r1,r2),(1) where vis the external potential, Nis the number of electrons and the notation r=|r1−r2|is used. The ground-state eigenfunction for this equation, ˜χ0say, corresponds to the pair density amplitude and is related to the pair density nof the real system as n=N(N−1) 2|˜χ0|2.(2) Eq. (1) contains an unknown term, vP, of completely kinetic origin. After a density functional analogy vPis called Pauli potential. Eq. (1) is analogous to the density functional equation for the square root of the density (which dates back to Thomas and Fermi [5] and is analyzed by Levy, Perdew and Sahni [6]). The pair density can be numerically calculated either on the Hartree-Fock level or on highly correlated level. The pair density can also be determined from Eq. (1) in a rather straightforward way if the Pauli potential is known. However, there are no data for the Pauli potential in the literature, yet. Although recently, the electron-electron cusp condition and asymptotic behaviour for the Pauli potential have been derived [7], much work is called for to completely understand how such a potential could be modellized in analogy with the Kohn-Sham potential in ordinary density functional theory (DFT). Considering the present knowledge, this kind of work appears extremely difficult. We believe that an important first step, in order to gain more insight in this direction, is the reconstruction of the Pauli potential 2 from a reliable pair density form for some realistic tractable electron system. For these cases, a six variable functional form for vPshould be, in principle, obtained. This function, at this point, should be viewed as a source of information expecially with the aim of finding those properties that can be transferred to other systems for which the pair density is unknowm. Motivated by the above consideration, in this paper we present model Pauli potential for the Be atom and isoelectronic atomic ions C2+ and O4+. The method, which is intended to capture the main features of vP, is based on an ansatz on the form of the pair density amplitude. We generalize the method of Amovilli et al. [8]. The original method was used to obtain the exact Hamiltonian for an analytic ground-state wave function for He-like ions. Here, a generalization is presented for producing the Pauli potential from a model pair density amplitude. The paper is organized as follows: In section 2 the pair density functional theory is reviewed. In section 3 a model pair density amplitude and the corresponding potential is presented. Section 4 describes numerical examples: the Be and some isoelectronic atomic ions. The last section is devoted to discussion. 2 The pair density functional theory First, the pair-density functional theory [1–3] is summarized. Consider the many electron Hamiltonian H, ˆ H=ˆ T+ˆ Vee + N X i=1 v(ri),(3) 3 where ˆ T= N X i=1 (−1 2∇2 i) (4) is the kinetic energy operator, ˆ Vee = N X i<j 1 |ri−rj|(5) is the electron-electron repulsion energy operator and v(r) is a local external potential. For convenience we consider an even number of particles. The second-ordered reduced density matrix is defined as n2(x1,x2;x0 1,x0 2) = N(N−1) 2ZΨ(x1,x2,x3..., xN)Ψ∗(x0 1,x0 2,x3..., xN)dx3...dxN,(6) where xistands for the spatial and the spin coordinates: ri, σiand the integral symbol when referred to spin denotes summation. The diagonal of the spin-independent second-ordered density matrix n(r1,r2) = X σ1,σ2 n2(r1, σ1,r2, σ2) (7) also called pair density is the key quantity. It is convenient to introduce new position variables qJ= (rj,rj0) = (qJ1, qJ2, qJ3, qJ4, qJ5, qJ6),(8) i.e. the pairs will be denoted by capital indices while the particles in each pair will be identified by the corresponding unprimed and primed letters. With the above notation the number of variables is the same as that of the initial system. Each particle is associated with a single pair, i.e. the number of indices Jis N/2. The ’internal’ potential for the particles in pair Jis given by ˜v(qj) = ˜v(rj,rj0) = 1 |rj−rj0|.(9) 4 while the interaction between pairs Iand Jis WIJ =W(qI,qJ) = W(ri,ri0;rj,rj0) = 1 |ri−rj|+1 |ri−rj0|+1 |ri0−rj|+1 |ri0−rj0|.(10) The energy of the pairs due to the external potential is ˆ U= M X I=1 u(qJ) = N X I=1 (v(rj) + v(rj0)) .(11) Defining the operator ˆ Lrepresenting the internal energy of pairs ˆ L= M X I=1 (−1 2∇2 I+ ˜v(qI)) ,(12) the initial Hamiltonian can be expressed as ˆ H=ˆ L+ˆ W+ˆ U , (13) where ˆ W=1 2 M X I6=J WIJ (14) is the interaction energy between different pairs (M=N/2). ˆ His the same as the initial Hamiltonian, but now it is written in terms of pairs of particles, with ˆ L+ˆ Urepresenting the Hamiltonian of independent (noninteracting with each other) pairs and ˆ Wrepresenting the interpair interaction. The Laplacian in the kinetic energy operator can also be written as ∇2 I=∇2 qI=∇2 i+∇2 i0= 6 X α=1 ∂2 ∂q2 Iα .(15) The energy of the independent pairs has the form Q[n] = min Ψ→nhΨ|ˆ L+ˆ W|Ψi.(16) 5 The search of the minimum is over all antisymmetric wave functions Ψ which yield the given n. Then the ground state energy can be written as E= min n1 N−1Zu(r1,r2)n(r1,r2)dr1dr2+Q[n].(17) The factor 1/(N−1) comes from the normalization of n. The density of pair I n(qI) = n(ri,ri0) = X σi,σi0 n2(ri, σi,ri0, σi0) (18) is the pair density in the original space. The Hohenberg-Kohn theorems [10] have been generalized for the pair density [11,12] of the original space. The ground state inequality is 1 N−1Zn(q)u(q) + Q[n]≥E0.(19) where E0and n0are the ground-state energy and the diagonal of the spin independent second-order density matrix, respectively. In the pair density functional theory the adiabatic connection is defined by the parametrized Hamiltonian ˆ Hα=ˆ L+αˆ W+ˆ Uα,(20) where ˆ Uα=PIuα I(q) is given by the condition that the pair density n(q), of the original space keeps being independent of α. For α= 0 the ’non-interacting Hamiltonian’ ˆ Hα=0 =ˆ L+ˆ Uα=0 =X I=1 hα=0 I(21) is obtained. In this auxiliary system the interaction between the pairs is zero and the auxiliary equations have the form ˆ H0Ψ0=E0Ψ0.(22) 6 The wave function in this auxiliary system can be written as a symmetrized expression of antisymmetric two-particle functions χI: Φ0(x1, ..., xN) = ˆ S(χ1(x1,x2)...χM(xN−1,xN)) .(23) ˆ S=1 N!X P ˆ P(24) is the symmetrizer operator. Pis the permutation operator and the sum is over all permutation of the electron pairs. This wave function is antisymmetric with respect to the exchange of the variables of a single pair and symmetric with respect to the exchange of the pairs. The disadvantage of the present notation is that it does not allow transposition of variables belonging to two different pairs. In the ground state n(q) = NN−1 2X σ |χ0(x1,x2)|2=NN−1 2|˜χ0(q)|2,(25) where the two-particle function ˜χ0satisfy the eigenvalue equation h0(q)˜χ0(q) = −1 2∇2 q+veff (q)˜χ0(q) = ε0˜χ0(q),(26) We mention in passing that it is possible to write the pair density in terms of geminals. The present version of the pair density version of theory has the advantage that the calculation of nis always reduced to the solution of a two-particle equation that is the N-body problem can be reduced to a two-body problem. It has been proved [1] that the auxiliary potential is uniquely determined by the diagonal form of the spin independent second-order density matrix and the effective potential is of the form veff (q) = v(r1) + v(r2) + N−1 r12 +vp,(27) 7 where vp= (N−1)δTP δn (28) and TP=T−T0(29) is the difference of the kinetic energies of the real system (T=hΨ|ˆ T|Ψi) and the auxiliary system T0= M X I=1 Zχ∗ I(x1,x2)[−1 2∇2 q]χI(x1,x2).(30) By a density functional analogy the functional TP[n] is called Pauli energy. The total energy has the form E[n] = T0[n] + TP[n] + Zn(q) r12 dq+1 N−1Zn(q)u(q)dq.(31) The disadvantage of the present treatment is that it is hard to capture the fermionic structure of an electronic system with a single effective potential. However, we have always a two-particle problem to solve independently of the number of electrons. It is worth to make efforts to find adequate approximation for the Pauli potential in order to utilize this benefit. The present study is a step in this direction. The auxiliary equations can also be derived by constrained search [1,13]. The two-particle equation (26) was later derived [14] in a different way which is not restricted to even number of electrons. 8 3 A model pair density amplitude and the corresponding potential As it was shown in the first paper [1] the Pauli potential is uniquely determined by the pair density. That is, from the knowledge of n,vPcan be given by inverting Eq. (1) vP(r1,r2) = −Kloc(r1,r2)−w(r1,r2),(32) where Kloc(r1,r2) = −1 2˜χ0(r1,r2)h∇2 1+∇2 2i˜χ0(r1,r2) (33) and w(r1,r2) = v(r1) + v(r2) + N−1 r−µ. (34) We have recently proved [7] that the Pauli potential asymptotically behaves as vP→(N−2) 1 r1 +1 r2 −1 r(35) when r1→ ∞,r2→ ∞ and r→ ∞. The electron-electron cusp condition has the form: vP=2−N r(36) as r→0. With the aim to reconstruct vPin some functional form, we start out from a model unnormalized pair density amplitude in the form χ=χHF (λr1, λr2)(1 + g(r)) ,(37) 9 a λ E Tp< r−2> < r−1> < r > < r2> < r3> 4 0.9835 –14.660(4) 0.878(3) 9.319 4.279 15.537 54.485 232.96 5 0.9847 –14.669(4) 0.864(3) 9.541 4.320 15.469 54.107 230.88 6 0.9860 –14.669(4) 0.855(3) 9.701 4.349 15.418 53.815 229.23 HF(a)–14.573 1.005 10.536 4.489 15.120 51.956 218.11 corr(b)–14.667 — 9.536 4.337 15.272 52.854 222.48 (a)In the Tpcolumn is reported the difference Tex −T(1) wand moments are from refs. [16,18]. (b)Moments from [18]. Table 2: Total energy (E), Pauli kinetic energy (Tp) and some moments < rk> for Be atom for different choices of the correlation function parameter aand the scaling constant λcalculated in this work and comparison with HF and correlated literature data. Data are in atomic units. a λ E Tp< r−2> < r−1> < r > < r2> < r3> 5 1.00090 –36.538(9) 3.253(8) 25.37 7.540 8.013 14.031 29.466 6 1.00063 –36.538(9) 3.226(8) 25.35 7.577 7.997 13.991 29.369 7 1.00055 –36.539(9) 3.204(8) 25.94 7.604 7.985 13.961 29.296 HF(a)–36.408 3.475 27.06 7.716 7.945 13.863 29.06 corr(b)–36.534 — 25.50 7.548 8.118 14.502 31.14 (a)In the Tpcolumn is reported the difference Tex −T(1) wand moments are from refs. [16,18]. (b)Moments from [16]. Table 3: Total energy (E), Pauli kinetic energy (Tp) and some moments < rk>for C2+ atomic ion for different choices of the correlation function parameter aand the scaling constant λcalculated in this work and comparison with HF and correlated literature data. Data are in atomic units. amplitude as they result from a fitting of the same accurate function. The complete definition of ˜χ0depends at this point by the parameter aentering the correlation function g(r) and the scaling factor λ. We made different choices of such parameters for all the three cases and the final results are collected in Tabs. 2,3,4. The main problem encountered in the calculation of the total energy by Monte Carlo method has been related to the high variance of the function to be averaged 16 a λ E Tp< r−2> < r−1> < r > < r2> < r3> 7 1.00005 -68.40(1) 6.91(1) 49.11 10.712 5.492 6.534 9.280 8 1.00003 -68.40(1) 6.91(1) 49.45 10.741 5.486 6.522 9.260 9 1.00003 -68.41(1) 6.90(1) 49.73 10.763 5.481 6.513 9.244 HF(a)-68.257 7.374 51.8 10.887 5.455 6.469 9.167 corr(b)-68.411 — 49.17 10.694 5.570 6.769 9.843 (a)In the Tpcolumn is reported the difference Tex −T(1) wand moments are from refs. [16,18]. (b)Moments from [16]. Table 4: Total energy (E), Pauli kinetic energy (Tp) and some moments < rk>for O4+ atomic ion for different choices of the correlation function parameter aand the scaling constant λcalculated in this work and comparison with HF and correlated literature data. Data are in atomic units. which is defined in Eq. (59). This requires a long simulation to achieve an energy mean value with an accuracy of the order of some mHartree. The same occurs for Tp. Looking at the results of Tabs. 2,3,4, it is evident that the optimal values of the parameters aand λmust be found by searching for a compromise between the need of getting reliable values of the moments < rk>and the best energy. The results show also that a considerable fraction of correlation energy has been taken into account. It is also important to notice that the variance becomes larger when the nuclear charge increases but also that our approximation, mainly based on short range correlation, should works better in such cases. It is also interesting to look at the values of Tp. From the bounds on the generalized Weizs¨acker-type kinetic energy introduced by Ayers [15] it follows that 0≤Tp≤Tex −T(1) w.(63) Looking at our results, this inequality is satisfied in the range of avalues consid17 ered here. Deviations from this behavior must lead to considerations related to N-representability. Finally, it is worthwhile to look at the plots of the effective potential derived by the approximate pair density amplitude. This has been done for the contributions V1(r1, r2) and V2(r) while V3(r1, r2, r) cannot be easily shown being dependent on three independent variables. For this purpose, and only for Be, we plot V1(r1, r2) in Fig.1 and V2(r) in Fig.2. From Fig.1, it is clear that V1is dominated by the external nuclear potential when r1or r2tends to 0 while is about constant for both large r1 and r2, being −µthe limit in this case. The ripples of the two dimensional surface of Fig.1 are instead a consequence of the exchange interaction and determine the shell structure of the one particle density of Be atom. Finally, V2, shown in Fig.2, is always repulsive. For small r, it behaves as the electron-electron interaction potential while it goes to zero more rapidly for large r. The ripples of V1, the long range behavior of V2and the contribution V3are special features of vP. 5 Conclusions In this work, we have illustrated a method to reconstruct the Pauli potential of pair density functional theory for four electron atomic ions. The potential is derived by inverting the effective two electron equation involving the pair density amplitude assuming that the pair density itself can be written in an analytical tractable form. Cusp and asymptotic conditions have been satisfied and appropriate adjustable parameters have been used in order to reproduce, within a reasonable accuracy, the total energy and some lower moment of the intracule density. Some interesting fea18 tures of the Pauli potential have been found for the systems treated here. These features are contained in the expressions (40), (44) and (45) for V1,V2and V3.V1and V2include also the nuclear and the electron-electron electrostatic potential energies. Some illustrations are given also in Figures 1 and 2. We would like to emphasize that the present method is not restricted to fourelectron sytems. In the the pair density theory one has to solve an effective two electron equation independently on the number of electrons. That is, the novel method introduced here to invert the effective two electron equation can always be applied if the the pair density (or the the pair density amplitude) is available. For the future, it will be interesting to analyze in details each individual term in order to find a generalization of the above expressions for all polyelectronic systems in a form which does not require the inversion of the effective two electron equation worked here. About V2, we would like to refer briefly to the ’average-pair-density theory’ of Gori-Giorgi and Savin. In this theory the spherically and system-averaged pair density f(r) is determined by simple radial equations conjectured by Gori-Giorgi and Savin [21]: h−∇2 r+weff (r)iφi(r) = iφi(r),(64) the solutions of which give f(r) as X i θi|φi(r)|2=f(r),(65) that is, f(r) is given by a weighted sum of the square of some orthogonal ’effective’ geminals φiwith weighting factors of ’occupancy’ θi. The potential weff (r) in Eq. 19 (64) was approximated as weff (r) = w(0) eff (r) + wc eff (r),(66) where w(0) eff (r) = ∇2f1/2 KS f1/2 KS (67) and wc eff (r) = 1 r+r2 2r3 s −3 2rs!θ(rs−r).(68) θ(rs−r) is the Heaviside step function and rs=4π 3%−1/3 ,(69) where %is the average electron density. The correlation potential wc eff (r) originally proposed by Overhauser [22], has been used to solve Eq. (64) for the uniform electron gas [21, 23]. It leads to an accurate description of the short-range part of f. Our potential V2(r1,r2) (44), using the expression (51) for g, has the form V2(r1,r2) = 1 r(1 + ar)2(1 + (a+ 1/2)r).(70) We immediatelly notice that the dominant term in (70) for small ris 1/r. It is the same as the first term in the Overhauser potential, which is also the dominant part of the Overhauser potential for small r. Thus the potential V2(r1,r2) has some resemblance to the Overhauser potential. The 1/r term in the Overhauser potential comes from the cusp condition on f(r) [21]. The dominant term in (70) for small r has the same origin. We also mention in passing that it was derived via a double adiabatic connection by one of the present authors [4] that the square root of the spherically and 20 system averaged pair density is the solution of a simple radial equation, that is, contrary to the theory of Gori-Giorgi and Savin, it is possible to obtain f(r) through a solution of a single equation. If a single geminal is used the fermionic character should be reflected in the potential which is consequently more complicated. If more than one geminals are used the spherically and system averaged pair density has a more complicated form but the potential can be more easily approximated. The number of geminals N(N−1)/2 depends on the number of electrons. Therefore a single geminal approach might gain an important role as the number of electrons increases. In the density functional theory there has been a growing interest in determining the exact exchange, exchange-correlation and Kohn-Sham potentials in the knowledge of the density. Several methods have been worked out [24–29]. The exact potentials are very useful, for example to check the accuracy of approximate methods. An analogous problem in the pair density functional theory is to obtain the Pauli potential in the knowledge of the pair density as here the electron-electron iteraction is exactly treated, but the kinetic energy functional is unknown. The problem here is more complicated in the sense that a two-particle potential vPshould be calculated, instead of a one-body exchange-correlation potential of the density functional theory. On the other hand, it is also simpler as only a single equation has to be inverted instead of several Kohn-Sham eguations in the density functional theory. The accurate form of the Pauli potential obtained by the present method can be used later to find approximate expressions for it. One has to be, however, extremely careful in the construction because of the N-representability problem [11,15,30–43]. Dal Ri et al. [19] derived density matrices from Jastrow-type trial wave func21 tions. The pair density used in this work can be considered as the lowest order approximation to the general, N-representable pair density presented by Dal Ri et al. Consequently, our pair density is, at least approximately , N-representable. References [1] ´ A. Nagy, Phys. Rev. A 66, 022505 (2002). [2] ´ A. Nagy: Pair Density Functional Theory, in The Fundamentals of Electron Density, Density Matrices and Density Functional Theory in Atoms, Molecules and Solid State, Eds. N. I. Gidopoulos and S. Wilson (Kluwer, 2003) p. 79. [3] ´ A. Nagy and C. Amovilli, J. Chem. Phys. 121, 6640 (2004). [4] ´ A. Nagy: J. Chem. Phys. 125, 184104 (2006). [5] L.H. Thomas, Proc. Cambr. Phil. Soc. 23, (1926) 542; E. Fermi, Z. Phys. 48 (1928) 73. [6] M. Levy, J. P. Perdew and V. Sahni, Phys. Rev. A 30, 2745 (1984). [7] ´ A. Nagy and C. Amovilli, J. Chem. Phys. 128, 114115 (2008). [8] C. Amovilli, N. H. March, I. A. Horward and ´ A. Nagy, Phys. Lett. A, 372, 4053 (2008). [9] C. Amovilli and N. H. March, Chem. Phys. Lett. 378, 167 (2003). [10] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). [11] P. Ziesche, Phys. Lett. A 195, 213 (1994). 22 [12] A. Gonis, T. C. Schulthess, J. van Ek and P. E. A. Turchi, Phys. Rev. Lett. 77, 2981 (1996); A. Gonis, T. C. Schulthess, P. E. A. Turchi and J. van Ek, Phys. Rev. B 56, 9335 (1997). [13] M. Levy, Proc. Natl. Acad. Sci. USA 76, 6062 (1979); E. Lieb, Int. J. Quantum. Chem. 24, 243 (1983). [14] F. Furche, Phys. Rev. A 70, 022514 (2004). [15] P. W. Ayers, J. Math. Phys. 46, 062107 (2005). [16] F. J. G´alvez, E. Buend´ıa and A. Sarsa, J. Chem. Phys. 111, 10903 (1999). [17] F. J. G´alvez, E. Buend´ıa and A. Sarsa, Chem. Phys. Lett. 378, 330 (2003). [18] J. Komasa, W. Cencek and J. Rychlewski, Phys. Rev. A 52, 4500 (1995). [19] M. Dal Ri, S. Stringari and O. Bohigas, Nucl. Phys. A376, 81 (1982). [20] M. Higuchi and K. Higuchi, Phys. Rev. A 75, 042510 (2007). [21] P. Gori-Giorgi and A. Savin, Phys. Rev. A 71, 032513 (2005). [22] A. W. Overhauser, Can. J. Phys. 73, 683 (1995). [23] P. Gori-Giorgi and A. Savin, Phys. Rev. A 73,032506 (2006); Phi. Mag. 86, 2643(2006); J. Toulouse, P. Gori-Giorgi and A. Savin, Int. J. Quantum. Chem. 106, 2026 (2006). [24] C.O. Almbladh and A.C. Pedroza, Phys. Rev. A 29, 2322 (1994). [25] F. Aryasetiawan and M.J. Stott, Phys. Rev. B 38, 2974 (1988); J. Chen and M.J. Stott, Phys. Rev. A 44, 2816 (1991); J. Chen, R.O. Esquivel and M.J. Stott, Phil. Mag. B 69, 1001, (1994). 23 [26] Q. Zhao and R.G. Parr, J. Chem. Phys. 98, 543 (1993); R.G. Parr, Phil. Mag. B69, 737 (1994); Q. Zhao, R. C. Morrison and R.G. Parr, Phys. Rev. A 50,2138 (1994). [27] R. van Leeuwen and E. J. Baerends, Phys. Rev. A 49,2421 (1994). [28] ´ A. Nagy, J. Phys. B 26, 43 (1993); Phil. Mag. B 69, 779 (1994). [29] A. G¨orling, Phys. Rev. A 46,3753 (1992);51,4501 (1995). [30] Reduced Density Matrices in Quantum Chemistry (Academic Press, New York, 1976). [31] P. Ziesche, Int. J. Quantum. Chem. 60, 149 (1996). [32] M. Levy and P. Ziesche, J. Chem. Phys.115, 9110 (2001). [33] A. J. Coleman, Rev. Mod. Phys. 35, 668 (1963); A. J. Coleman and V. I. Yukalov, Reduced Density Matrices: Coulson’s Challange (Sringer-Verlag, New York, 2000); J. Cioslowski, Many-electron Densities and Reduced Density Matrices (Kluwer/Plenum, New York, 2000). [34] K. Husimi, Proc. Phys. Math. Soc. Jpn. 22, 264 (1940). [35] P. O. L¨owdin, Phys. Rev. 97, 1474 (1955). [36] E. R. Davidson, Chem. Phys. Lett. 246, 209 (1995); Phys. Rev. A 1, 30 (1970). [37] F. Sasaki, Phys. Rev. 138, B1338 (1965). [38] J. K. Percus, J. Chem. Phys. 122, 234103 (2005); B. Liu and J. K. Percus, Phys. Rev. A 74, 012508 (2006). [39] M. E. Pistol, Chem. Phys. Lett. 400, 548 (2004); Chem. Phys. Lett. 417, 521 (2006). 24 [40] P. W. Ayers and E. R. Davidson, Int. J. Quantum. Chem. 106, 1487 (2006). [41] P. W. Ayers, S. Gordon and M. Levy, J. Chem. Phys. 124, 054101 (2006). [42] P. W. Ayers, and M. Levy, J. Chem. Sci. 117, 507 (2006). [43] C. Garnod and J. K. Percus, J. Math. Phys. 5, 1756 (1964). Acknowledgements This paper was written in the frame of the Bilateral Scientific Cooperation between Italy and Hungary sponsored by Consiglio Nazionale delle Ricerche and the Hungarian Academy of Sciences. Grant OTKA No. T 029469 is gratefully acknowledged. 25