Proposal for a modified Møller-Plesset perturbation theory

Proposal for a modified Møller-Plesset perturbation theory Alejandro Cabo,1Francisco Claro,2Eduardo Menéndez-Proupin,3,*Norge Cruz-Hernández,4and Javier Fernández-Sanz4 1Grupo de Física Teórica, Instituto de Cibernética, Matematica y Física, Calle E, No. 309, Vedado, La Habana, Cuba 2Facultad de Física, Pontificia Universidad Católica de Chile, Vicuña Mackenna 4860 Santiago, Chile 3Departamento de Física, Facultad de Ciencias, Universidad de Chile, Las Palmeras 3425, 780-0024 Ñuñoa, Santiago, Chile 4Facultad de Química, Departamento de Química Física, Universidad de Sevilla, E-41012 Sevilla, Spain 共Received 2 June 2005; published 23 January 2006兲 A modified version of the Møller-Plesset approach for obtaining the correlation energy associated with a Hartree-Fock ground state is proposed. The method is tested in a model of interacting fermions that allows for an exact solution. Using up to third order terms improved results are obtained, even in the limit of loosely bound particles. Tested in molecules as well, the modified method appears to give improved results in symmetric systems. DOI: 10.1103/PhysRevA.73.012510 PACS number共s兲: 31.15.Md, 03.65.Ge, 31.25.⫺v, 02.70.⫺c I. INTRODUCTION The study of interacting many-particle systems is seriously constrained by the large dimension of the Hilbert space. Several approximation schemes have been devised over the years, among which the Hartree-Fock 共HF兲method is one of the oldest and most fruitful, notably in atomic and molecular physics. Because it treats interactions in a mean field way particle correlations are left out, however, a shortcoming that can limit severely the validity of its results. One may improve over HF by treating correlations as a perturbation. In the so-called Møller-Plesset method, a Rayleigh- Schrödinger perturbative expansion that is naturally suggested by the same structure of the HF solution 关1–4兴is adopted. This formalism reduces the correlation energy to an infinite series in the perturbation, of which only the first few terms need to be computed in practice. This scheme has been used for a long time as a good starting point to study correlation effects in molecular systems. However, this method is useful only if the perturbation series is rapidly convergent, which is not always true 关5兴. Failures of the Møller-Plesset method have been documented even for small molecules 关5–7兴. In this paper we present a variation of the Møller- Plesset approach that in some cases appears to give accurate results in low-order perturbation schemes. Consider the Hamiltonian of a general system of identical particles interacting through a pair potential H=兺 i h共i兲+兺 i,j⬎i v共i,j兲,共1兲 where h共i兲=h共ri兲is the sum of the one-particle kinetic energy plus external potential energy, and v共i,j兲=v共ri,rj兲is the pair-interaction potential. The HF approximation to the ground state of a system of identical fermions leads to a variational wave function ⌽in the form of a Slater determinant, constructed with one-particle orbitals that satisfy the self-consistent field 共SCF兲equations 再 h共i兲+兺 b 关Fb共i兲−Kb共i兲兴 冎 ␾ n共i兲= ⑀ n ␾ n共i兲,共2兲 where Fband Kbare the Coulomb and exchange operators, respectively 关4兴. Here and in what follows the summation over indices a,b,crun over all occupied states, while nruns over all possible states. The eigenvalues satisfy the relation ⑀ n=具n兩h兩n典+兺 b 具nb储nb典,共3兲 where 具n兩h兩n典=兰 ␾ n共1兲*h共1兲 ␾ n共1兲d1 are the matrix element of one-particle operator hin the basis 兵 ␾ n其, and 具mn储mn典=具mn兩mn典−具mn兩nm典,共4兲 具mn兩pq典= 冕 ␾ m共1兲* ␾ n共2兲*v共1,2兲 ␾ p共1兲 ␾ q共2兲d1d2共5兲 are the matrix elements of the pair-interaction operator. In the above expression 1 and 2 represent the one-electron variables of coordinate and spin. The energy of the HF state is EHF =兺 a 具a兩h兩a典+1 2兺 a,b 具ab储ab典=兺 a ⑀ a−1 2兺 a,b 具ab储ab典, 共6兲 where to obtain 共6兲Eq. 共3兲is used 关4兴. The many-particles operator HHF =兺 n ⑀ nc ˆn †c ˆn−1 2兺 a,b 具ab储ab典共7兲 is a natural choice for an approximate independent-particles description of the interacting system. In Eq. 共7兲c ˆn †共c ˆn兲is the creation 共annihilation兲operator of a particle in the state ␾ n. HHF, known as Hartree-Fock Hamiltonian, is diagonal in the basis of all the Slater determinants composed with the HF orbitals, and its ground state energy is the HF energy EHF. This Hamiltonian operator is the starting point to construct a perturbation expansion for the correlation energy using as perturbation VHF=H−HHF 关4兴. Hence, the energy of the interacting electron system is obtained as a series *Electronic address: [email protected] PHYSICAL REVIEW A 73, 012510 共2006兲 1050-2947/2006/73共1兲/012510共5兲/$23.00 ©2006 The American Physical Society012510-1 E=E共0兲+E共1兲+E共2兲+E共3兲+¯.共8兲 With the choice 共7兲for the zero-order Hamiltonian, we have that E共0兲=EHF, and E共1兲=0. A well known feature of the HF theory is that the many body ground state energy EHF is not the bare sum of energies of filled single-particle orbitals. Interactions are counted twice and this overestimation is corrected by subtracting the constant explicit in Eq. 共6兲. A second choice for a Hamiltonian is also possible, however. Equation 共6兲may be written in the alternative form, EHF =兺 a 冉 ⑀ a−1 2兺 b 具ab储ab典 冊 ,共9兲 suggesting as an alternative Hamiltonian also diagonal in the HF orbitals HMHF =兺 n ⑀ ˜ nc ˆn †c ˆn,共10兲 where ⑀ ˜ n= ⑀ n−1 2兺 b 具nb储nb典=1 2共 ⑀ n+具n兩h兩n典兲.共11兲 In 共10兲the operators c ˆn †and c ˆnare the same as in 共7兲. Hence, HMHF and HHF have the same eigenfunctions, and the same ground state energy, but different excitation energies. The correlation energy can again be described in terms of a perturbation, this time of the form VMHF=H−HMHF. For instance, the second-order correction to the HF ground state energy has the form 关4兴 E共2兲=兺 ⌽⬘ 円具⌽兩H兩⌽⬘典円2 EHF −E⌽⬘ =1 4兺 abrs 兩具ab储rs典兩2 EHF −E⌽⬘ .共12兲 This is the first finite correction in the Rayleigh-Schrödinger perturbation expansion since for either choice of the Hamiltonian the first order term vanishes. Only doubly excited states ⌽⬘=⌽ab rs 共 ␾ areplaced by ␾ rand ␾ breplaced by ␾ sin the Slater determinant ⌽, with ⑀ r, ⑀ sabove the Fermi energy兲 contribute 关4兴. The choice of the zero-order Hamiltonian affects the excitation energies in the denominator, which have the form EHF −E⌽ab rs = 再 ⑀ a+ ⑀ b− ⑀ r− ⑀ s,if H0=HHF, ⑀ ˜ a+ ⑀ ˜ b− ⑀ ˜ r− ⑀ ˜ s,if H0=HMHF. 冎 共13兲 Replacing in Eq. 共12兲both forms clearly lead to different numerical values. The analysis of the third- and higher-order terms in the perturbation expansion yields additional modifications. In the third-order correction E共3兲, besides the replacements of ⑀ m→ ⑀ ˜ min the standard Møller-Plesset expression 关4兴, an additional term ⌬E共3兲appears ⌬E共3兲=−E共2兲−1 4兺 abrs 共hrr +hss −haa −hbb兲兩具ab储rs典兩2 共 ⑀ ˜ r+ ⑀ ˜ s− ⑀ ˜ a− ⑀ ˜ b兲2, 共14兲 where hnn=具n兩h兩n典. Curiously, the first element of Eq. 共14兲 cancels the second-order energy correction. Yet, the sum that follows contains the second-order energy as well, with a plus sign, so that ⌬E共3兲is third order in the interaction. TABLE I. Ground state energies of a system of two harmonically confined spin-1 2particles interacting through a harmonic potential of strength k. Several approximations are included: Hartree-Fock 共HF兲, Møller- Plesset perturbation theory of orders n=2, 3 共MPn兲, modified MPn共MMPn兲, as well as the exact values. Repulsive interactions are represented by negative values of the elastic constant k. kHF MP2 MP3 MMP2 MMP3 Exact −0.25 1.732 1.655 1.836 1.702 1.710 1.707 −0.24 1.744 1.681 1.803 1.717 1.724 1.721 −0.22 1.766 1.725 1.784 1.745 1.750 1.748 −0.20 1.789 1.760 1.791 1.772 1.776 1.775 −0.18 1.811 1.791 1.808 1.798 1.801 1.800 −0.16 1.833 1.819 1.828 1.823 1.825 1.825 −0.09 1.908 1.905 1.906 1.905 1.906 1.906 −0.04 1.960 1.959 1.959 1.959 1.959 1.959 −0.01 1.990 1.990 1.990 1.990 1.990 1.990 0.00 2.000 2.000 2.000 2.000 2.000 2.000 0.04 2.040 2.039 2.039 2.039 2.039 2.039 0.16 2.154 2.150 2.149 2.149 2.149 2.149 0.36 2.332 2.319 2.314 2.316 2.313 2.311 0.64 2.561 2.534 2.522 2.525 2.516 2.510 1.00 2.829 2.784 2.762 2.767 2.749 2.732 CABO et al. PHYSICAL REVIEW A 73, 012510 共2006兲 012510-2 In order to assess the convenience of either formulation for obtaining the correlation energy we have applied them to two systems of known exact solutions, and for a small number of molecules. The following sections are devoted to this analysis. II. EXACT MODELS Our first model system is a two spin-1 2particle in a bidimensional harmonic potential, interacting through a harmonic force. This is an interacting system involving identical fermions that is solved exactly 关8兴. The Hamiltonian is given by 共1兲with h共i兲=1 2共−ⵜi 2+ri 2兲,v共i,j兲=1 2k共ri−rj兲2.共15兲 Here, the coordinates and the energy are given in the oscillator units of the confinement potential. The exact ground state energy is E0=1+冑1+2k.共16兲 The model contains the parameter kthat allows the study of attractive 共k⬎0兲as well as repulsive 共k⬍0兲interactions. Equation 共16兲shows that k=−0.5 is the lowest value for which a bound state exists. To obtain a numerical solution at the Hartree-Fock level, we use a basis of noninteracting harmonic oscillators eigenfunctions ␾ nx,ny共x,y兲= ␸ nx共x兲 ␸ ny共y兲,0艋nx+ny艋5, 共17兲 where the ␸ 共x兲are the usual one-dimensional harmonic oscillator eigenfunctions. The set 共17兲is an exact solution when k=0. The HF solution is obtained by solving the SCF equations 共2兲in a spin restricted configuration, i.e., constraining the spatial part of the wave function to be equal for the two occupied states with opposite spins 关4兴. Table I shows the ground state energy for the Hamiltonian 共15兲calculated using the standard Møller-Plesset perturbation theory of order n=2, 3 共MPn兲, and our modified form, MMPn. Figure 1 shows the exact correlation energy Ecorr =EHF−Efor different values of k, together with results obtained for the two choices of the HF Hamiltonian, in the second and third order of the perturbation theory. Notice that MMPnyields better results throughout. Notice also that the usual MPnfails poorly in both orders of approximation when the system becomes more loosely bound, as kapproaches the critical value −0.5. In fact, convergence problems prevent solving the HF SCF equations for kbeyond −0.25. By contrast, the MMPncontinue to be good approximations even in this range of k. The simplest molecular system is the hydrogen molecule, for which we present the potential energy surface 共PES兲in both spin restricted and unrestricted configurations 关8兴. The TABLE II. Bondlengths 共in pm兲of diatomic molecules and the error made in different approximations. The basis used is cc-pVTZ. Exp HF MP2 MMP2 B3LYP H274.1 −0.7 −0.4 −0.2 +0.2 HF 91.7 −1.9 +0.03 −1.6 +0.5 OH+102.9 −2.2 −0.5 −1.6 +1.0 NH 103.6 −1.9 −0.9 −1.7 +0.3 NO+106.3 −3.6 +1.5 +2.7 −0.6 BH 123.2 −1.0 −1.5 −0.8 +0.0 FIG. 2. Potential energy surface for the H2molecule as a function of the distance H-H. The basis set used is 6-31G**. Top: spinrestricted configuration. Bottom: spin-unrestricted configuration. FIG. 1. Correlation energy for a system of two harmonically confined spin-particles interacting through a harmonic potential of strength k. The exact result is included, as well as corrections up to second and third order for the Møller-Plesset 共MP兲and modified Møller-Plesset 共MMP兲choices of zeroth-order Hamiltonian. PROPOSAL FOR A MODIFIED MØLLER-PLESSET …PHYSICAL REVIEW A 73, 012510 共2006兲 012510-3 molecular electronic Hamiltonian in atomic units is given by 共1兲with h共i兲=−1 2ⵜi 2,v共i,j兲=1 兩ri−rj兩.共18兲 The HF wave functions are obtained by the expansion in standard basis sets of Gaussian functions 关4兴共e.g., 6-31G**, DZP, etc.兲and applying the SCF method as implemented in the code HONDO 关9兴. The MMP2 and MMP3 energies are obtained by a minor modification of the MP2 and MP3 routines in the same code. The PES of H2has been calculated by Kolos and Wolniewicz 关10兴with extreme precision and can be regarded as exact, thus allowing a reference for approximative methods. It can be seen in Fig. 2 that the MMP2 method improves the correlation energy over the standard MP2 in the whole range of distances. A well known failure of the spin restricted HF and derived methods is the inability to describe the dissociation of the molecule, as the two occupied states must have the same orbital wave function. This problem is solved by the spin-unrestricted HF method, where the restriction that opposite spin orbitals share the orbital part is lifted. Again, our modified scheme provides a solution closer to the exact one to second order. However, it can be appreciated that for unrestricted calculations both the MMP2 and MP2 are worse than the restricted case in the range 1.2–1.6 Å. This shows that the problem of spin contamination affects both MP2 and MMP2. III. MOLECULAR SYSTEMS In this section, we evaluate the MMP approximation for several simple molecules, for which no exact wave functions are known. Table II shows the error made at the second order of approximation in the bondlength of a set of diatomic molecules, as compared with the experimental length. It can be appreciated that the bondlength is improved by the perturbation schemes saved for BH, where the MP2 increases the error. MMP2 gives better results for H2and BH, while MP2 is more accurate in the other cases. For reference, the results with the popular density functional method B3LYP are included in the table. Next, we turn our attention to the value of the correlation energy at the experimental geometries. Tables III and IV show a comparison of the standard MPnenergies with the present MMPnenergies. For reference and comparison, the values with the most accurate configuration interaction 共CI兲 quantum chemistry method are included. For these molecules it turns out that the correlation energy is, in general, better in MPnthan in MMPnmethods. Only for the symmetric pairs the methods are comparable, MMPngiving a slightly better result for H2. The reason for the rather low correlation energy in the MMP2 can be explained in terms of Eqs. 共12兲and 共13兲.It corresponds to an increase of the energy gap between the occupied and unoccupied states. This can be appreciated in Fig. 3, where the Fock eigenvalues 共3兲and the modified eigenvalues 共11兲for the HF are shown. As can be noticed, the correction to the Fock eigenvalues is greater for the 共dou- TABLE III. Correlation energy 共in atomic units兲for several diatomic molecules at the experimental bondlengths. The basis used is cc-pVTZ. EHF MP2 MP3 MMP2 MMP3 CI H2−1.133 −0.032 −0.037 −0.034 −0.038 −0.039 HF −100.058 −0.290 −0.290 −0.228 −0.268 −0.282 C2−75.402 −0.385 −0.347 −0.328 −0.346 −0.407 NH −54.875 −0.172 −0.191 −0.149 −0.179 −0.207 NO+−128.966 −0.439 −0.422 −0.369 −0.408 −0.445 OH+−74.866 −0.182 −0.202 −0.156 −0.189 −0.222 BH −25.130 −0.084 −0.101 −0.077 −0.095 −0.110 TABLE IV. Calculated energies 共in atomic units兲for several diatomic and triatomic molecules at the experimental geometries. Geometries and basis functions are specified as follows: BH: r=2.329 a0, DZP. HF: r=1.733 a0,DZP.CH 2:r=2.11 a0, ␪ =102.4°, DZP. H2O: r=1.88973 a0, ␪ =104.5°, DZP. NH3: r=1.911 65 a0, ␪ =106.7°, DZV. a0: Bohr radius, DZV: Dunning’s 共9s,5p兲/共3s,2p兲basis set, DZP: DZV plus polarization. EHF MP2 MP3 MMP2 MMP3 CI BH −25.124 −25.183 −25.199 −25.187 −25.198 −25.209 HF −100.048 −100.232 −100.233 −100.186 −100.220 −100.251 CH2−38.885 −38.998 −39.017 −38.992 −39.013 −39.027 H2O −76.041 −76.242 −76.247 −76.200 −76.234 −76.254 NH3−56.176 −56.290 −56.297 −56.268 −56.289 −56.304 CABO et al. PHYSICAL REVIEW A 73, 012510 共2006兲 012510-4 bly兲occupied levels 共level艋5兲than for the unoccupied ones. This is related to the values of the interaction integrals 具nb储nb典present in Eq. 共11兲. These integrals are generally positive and have larger values for the occupied levels because the corresponding orbitals ␾ nhave more overlap with the occupied orbitals 共 ␾ b兲. For example, in hydrogen flouride the occupied orbitals are concentrated around the F atom, while the lowest unoccupied molecular orbital and higher molecular orbitals are more dense around the H atom. For the C2molecule the effect is less pronounced though qualitatively similar. The case of H2and the two-fermion oscillator is different. In this case, there are only two occupied states with the same energy, and the sum over the bstates in Eq. 共11兲contains only one term, as 具nn储nn典=0. On the other hand, for the unoccupied states the bsum contains two terms and the negative shift of these eigenvalues is larger, thus reducing the gap. This effect is also present in larger molecules, but it is overcome by the differences in the spatial distribution of the wave functions, as discussed above. IV. CONCLUSIONS We have proposed and analyzed critically an alternative formulation of the Møller-Plesset perturbation series for the ground state of many-particle fermion systems. The new zeroth-order Hamiltonian has the same wave functions as the Hartree-Fock Hamiltonian, but with a different spectrum. Here, the correction to obtain the correct Hartree-Fock energy is included in the single-particle energies, instead of through a global energy shift. This opens a possibility of including part of the correlation energy in the zero-order Hamiltonian, correcting at the same time for the well known double count of the interaction energy present in the ordinary Hartree-Fock formulation. We have applied the new method to a number of different systems and have found that the new method is as good or better than the Møller-Plesset scheme in symmetric systems. However, for molecular systems with more than two electrons, the standard MP series seems to provide better accuracy. It turns out that the correction to the one-particle energy levels is larger for the occupied levels than for the unoccupied ones, increasing the gap between occupied and unoccupied states. This reduces the size of the lowest-order contributions to the correlation energy, decreasing the rate of convergence of the perturbation series. However, we note that in order to obtain the HF energy in the zero-order approximation to the many-particles problem, Eq. 共11兲needs to be applied only to the occupied HF states. We stress that this choice is motivated by the attempt to eliminate the double counting of the mean field interaction energies in the HF energy. There is no reason other than simplicity to apply the same scheme to the empty levels. This degree of freedom could be exploited to improve the convergence of the many-particles perturbation series. ACKNOWLEDGMENTS Support from FONDECYT Grant Nos. 1020829, 1050293, 7020829, and the Third World Academy of Sciences, is gratefully acknowledged. N.C-H. thanks the Spanish Ministerio de Ciencia y Tecnología–Ramón y Cajal program. We also thank M. Dupuis for providing his code HONDO, and P. Fuentealba and J. Garza for helpful discussions. 关1兴C. Møller and M. S. Plesset, Phys. Rev. 46, 618 共1934兲. 关2兴J. M. Eisenberg and W. Greiner, Microscopic Theory of the Nucleus 共North-Holland, Amsterdam, 1972兲. 关3兴N. H. March, Self-Consistent Fields in Atoms 共Pergamon Press, Oxford, 1975兲. 关4兴A. Szabo and N. S. Ostlung, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory 共Dover Publications Inc., New York, 1996兲. 关5兴R. H. Nobes, J. A. Pople, L. Radom, N. C. Handy, and P. J. Knowles, Chem. Phys. Lett. 138, 481 共1987兲. 关6兴R. H. Nobes, D. Moncrieff, M. W. Wong, L. Radom, P. M. W. Gill, and J. A. Pople, Chem. Phys. Lett. 182, 216 共1991兲. 关7兴N. L. Ma, B. J. Smith, and L. Radom, Chem. Phys. Lett. 193, 386 共1992兲. 关8兴M. Moshinsky, O. Novaro, and A. Calles, J. Phys. 共Paris兲31, C4 共1970兲. 关9兴M. Dupuis, HONDO-2004, based on HONDO-95, available from the Quantum Chemistry Program Exchange, University of Indiana, Bloomington, IN 47405;M. Dupuis, J. D. Watts, H. O. Villar, and G. J. B. Hurst, Comput. Phys. Commun. 52, 415 共1989兲;134, 150 共2001兲;M. Dupuis and A. Marquez, J. Chem. Phys. 114, 2067 共2001兲. 关10兴W. Kolos and L. Wolnievicz, J. Chem. Phys. 49, 404 共1968兲. FIG. 3. Fock eigenvalues 共upper group兲and modified eigenvalues 共lowest group兲of the molecule HF, calculated with different basis sets: STO-3G, 6-31G, and cc-pVTZ. The lines are guides to the eyes. The dashed vertical line indicates the highest occupied molecular orbital. PROPOSAL FOR A MODIFIED MØLLER-PLESSET …PHYSICAL REVIEW A 73, 012510 共2006兲 012510-5