On the Estimation of the Remainder Term in Møller-Plesset MP2 Theory from Limited Configuration Interaction

Abstract

In a previous paper, we have avoided an infinite order perturbation expansion and obtained a closed expression which consists of the second-order Møller-Plesset energy component together with a remainder term. The applicabilty of second-order many-body perturbation theory with a Møller-Plesset reference hamiltonian then rests upon the magnitude of this remainder term rather than the behaviour of the higher order terms on the perturbation series. In the present work, we show how this remainder term can be estimated by limited configuration interaction.

1 Introduction

Many-body perturbation theory with a Møller-Plesset reference hamiltonian is the most widely used approach to the correlation problem in atomic and molecular systems. Second-order theory, which is often designated “MP2” and which incidentally was the order of theory originally pre-sented by Møller and Plesset [1], is computationally efficient and facilitates the use of very large basis sets which allows basis set truncation errors to be reduced to a level where other effects, and in particular relativity, are more significant [2,3,4].

By employing a hybrid partitioning scheme based on both the Rayleigh-Schrödinger and the generalized Brillouin-Wigner approaches, it has been shown [5] that the electron correlation en-ergy expression arising in the second-order many-body perturbation theory with a Møller-Plesset reference hamiltonian can be written as part of a closed expression consisting of the MP2 ap-proximation to the correlation energy together with a remainder term, ℜ. The utility of the MP2 theory therefore rests not on the behaviour of the individual higher order terms in the perturbation expansion but on the magnitude of the remainder term.

In this work we obtain an estimate of the remainder term, ℜ, from limited configuration inter-action and, in particular, from the method most often designated CISD, configuration interaction with single and double excitations. In section 2, we define the remainder term in Møller-Plesset second-order theory for the energy. In section 3, we used Brillouin-Wigner perturbation theory to analyse the limited configuration interaction expansion for a non-degenerate ground state and obtain an estimate for the remainder term ℜ. This is followed by a brief discussion.

2 The remainder term in MP2 theory

Let the time-independent Schrödinger equation be written in the form

$$H|\Psi \rangle =(H_0+\lambda H_1)|\Psi \rangle =\epsilon |\Psi \rangle$$

(1)

with the usual normalization condition $\langle \Psi |\Psi \rangle =1$. Approximations to the eigenvalues of (1) are developed with respect to the solutions of a model eigenproblem

$$H_0|{\Phi }_i\rangle =E_i|{\Phi }_i\rangle$$

(2)

with $\langle {\Phi }_i|{\Phi }_j\rangle ={\mathsf{\delta }}_{ij}$. H₁ is the perturbation operator and λ is the perturbation parameter. Define the projector

$$P=|{\Phi }_0\rangle \langle {\Phi }_0|$$

(3)

and its orthogonal complement Q = I − P. Assuming the intermediate normalization condition, $\langle {\Phi }_0|\Psi \rangle =1$, we have

$$|{\Phi }_0\rangle =P|\Psi \rangle$$

(4)

Following the partitioning approach pioneered by Löwdin [6,7,8,9,10,11,12,13,14,15,16,17,18] , we can use (1) and (3) and write

$$H_{\text{eff}}|{\Phi }_0\rangle =\epsilon |{\Phi }_0\rangle$$

(5)

where the effective hamiltonian is given by

$$H_{\text{eff}}=PHP+PHQ{(\epsilon -QHQ)}^{-1}QHP$$

(6)

The exact energy is then $\epsilon =\langle {\Phi }_0|H_{\text{eff}}|{\Phi }_0\rangle$. Different types of perturbation theory may then be obtained by expanding the denominator in (6) using the operator recursion

$${(X-Y)}^{-1}=X^{-1}+X^{-1}Y{(X-Y)}^{-1}$$

(7)

In our previous work [5], we point out that equation (7) can be written in the more general form

$${(X-Y)}^{-1}=X^{-1}+X^{-1}Y{(X^{\prime }-Y^{\prime })}^{-1}$$

(8)

where, in general, $X\ne X^{\prime }$ and $Y\ne Y^{\prime }$ but

$$X-Y=X^{\prime }-Y^{\prime }$$

(9)

and in the application of (8) to the effective hamiltonian (6)

$$X-Y=X^{\prime }-Y^{\prime }=\epsilon -QHQ$$

(10)

Putting

$$X=E_0-Q{H}_{0}Q$$

(11)

and

$$Y=Q{H}_{1}Q-\Delta E_0$$

(12)

Where $\Delta E_0=\epsilon -E_0$, yields the Rayleigh-Shcrödinger perturbation series, whilst the choice

$$X^{\prime }=\epsilon$$

(13)

and

$$Y^{\prime }=QHQ$$

(14)

leads to the generalized Brillouin-Wigner expansion, which was introduced by Löwdin in Parts II and XII of his studies of perturbation theory [6,7,8,9,10,11,12,13,14,15,16,17,18] and re-examined recently by the present authors[19].

Employing both the Rayleigh-Schrodinger choice for X and Y, that is (11) and (12), and the generalized Brillouin-Wigner choice for $X^{\prime }$ and $Y^{\prime }$, that is (13) and (14), in the recursion (8) gives

$$\begin{array}{ll}{(\epsilon -QHQ)}^{-1}=& {(E_0-Q{H}_{0}Q)}^{-1}+\\ & {(E_0-Q{H}_{0}Q)}^{-1}(Q{H}_{1}Q-\Delta E_0){(\epsilon -QHQ)}^{-1}.\end{array}$$

(15)

The exact energy eigenvalue is given by

$$\begin{array}{ll}\epsilon & =\langle {\Phi }_0|H_{eff}|{\Phi }_0\rangle \\ & =\langle {\Phi }_0|H|{\Phi }_0\rangle +\langle {\Phi }_0|HQ{(\epsilon -QHQ)}^{-1}QH|{\Phi }_0\rangle \end{array}$$

(16)

which upon the substitution of (15) gives

$$\epsilon =\langle {\Phi }_0|H|{\Phi }_0\rangle +\langle {\Phi }_0|HQ{(E_0-Q{H}_{0}Q)}^{-1}QH|{\Phi }_0\rangle +\Re$$

(17)

where the first term on the right-hand-side is the sum of the zero-order and first order energies

$$\begin{array}{ll}\langle {\Phi }_0|H|{\Phi }_0\rangle & =\langle {\Phi }_0|H_0|{\Phi }_0\rangle +\langle {\Phi }_0|H_1|{\Phi }_0\rangle \\ & =E_0+{\epsilon }_1^{(0)}\end{array}$$

(18)

which in Møller-Plesset theory is just the Hartree-Fock energy. The second term on the right-hand-side of equation (17) is the second-order Møller-Plesset energy

$$\begin{array}{ll}\langle {\Phi }_0|HQ{(E_0-Q{H}_{0}Q)}^{-1}QH|{\Phi }_0\rangle & ={\displaystyle \sum _{k\ne 0}\frac{\langle {\Phi }_0|H_1|{\Phi }_k\rangle \langle {\Phi }_k|H_1|{\Phi }_0\rangle }{E_0-E_k}}\\ & ={\epsilon }_2^{(0)}\end{array}$$

(19)

Unlike the usual perturbation expansions, (17) is not an infinite series, but a closed expression containing a remainder term which we designate ℜ and which has the form

$$\Re =\langle {\Phi }_0|HQ{(E_0-Q{H}_{0}Q)}^{-1}(Q{H}_{1}Q-\Delta E_0){(\epsilon -QHQ)}^{-1}QH|{\Phi }_0\rangle$$

(20)

or

$$\Re =\langle {\Phi }_0|HQ{(E_0-Q{H}_{0}Q)}^{-1}(Q{H}_{1}Q-(\epsilon -E_0)){(\epsilon -QHQ)}^{-1}QH|{\Phi }_0\rangle$$

(21)

The approximation to the total energy of an atomic or molecular system given by MP2 theory can thus be written in the form

$$\epsilon =E_0+{\epsilon }_1^{(0)}+{\epsilon }_2^{(0)}+\Re$$

(22)

The utility of the method rests on the requirement that the remainder term satisfies

$$\Re <\tau$$

(23)

for some arbitrarily chosen small τ. Specifically, the use of the MP2 method does not depend on the behaviour of the individual higher order terms in the Møller-Plesset perturbation expansion.

3 Estimation of the remainder term from limited configuration interaction

In this paper, we consider the use of limited configuration interaction in obtaining an estimate for $\Re$. We restrict our attention to a non-degenerate ground state and employ the Brillouin-Wigner perturbation expansion [20,21,22,23] in our analysis. The energy of such a system is given by

$$\epsilon =\langle {\Phi }_0|H_0|{\Phi }_0\rangle +\langle {\Phi }_0|H_1+H_1{ℬ}_0{H}_1+\dots |{\Phi }_0\rangle$$

(24)

where $ℬ$₀ is the Brillouin-Wigner resolvent

$${ℬ}_0={\displaystyle \sum _{k\ne 0}\frac{|{\Phi }_k\rangle \langle {\Phi }_k|}{\epsilon -E_k}}$$

(25)

Introducing the reaction operator, ${\nu }_B$, as

$${\nu }_B=H_1+H_1{ℬ}_0{H}_1+\dots$$

(26)

we can write (24) as

$$\epsilon =\langle {\Phi }_0|H_0|{\Phi }_0\rangle +\langle {\Phi }_0|{\nu }_B|{\Phi }_0\rangle$$

(27)

${\nu }_B$ can also be written as

$${\nu }_B=H_1+H_1{ℬ}_0{\nu }_B$$

(28)

which is a Lippmann-Schwinger-like equation [24]. The wave operator [25,26,27] is defined by

$$|\Psi \rangle =\mathsf{\Omega }|{\Phi }_0\rangle$$

(29)

so that, since from (1) and (4) we have

$$\begin{array}{ll}\epsilon & =\langle {\Phi }_0|H|\Psi \rangle \\ & =\langle {\Phi }_0|H_0|{\Phi }_0\rangle +\langle {\Phi }_0|H_1|\Psi \rangle ,\end{array}$$

(30)

we can compare with equation (27) and write

$${\nu }_B=H_1\mathsf{\Omega }.$$

(31)

Furthermore, using (28), we can write

$$\mathsf{\Omega }=1+{ℬ}_0{H}_1\mathsf{\Omega }.$$

(32)

Substituting (31) into (27) gives

$$\epsilon =\langle {\Phi }_0|H_0|{\Phi }_0\rangle +\langle {\Phi }_0|H_1\mathsf{\Omega }|{\Phi }_0\rangle$$

(33)

and then, by introducing (32) into (33) we get

$$\epsilon =\langle {\Phi }_0|H_0|{\Phi }_0\rangle +\langle {\Phi }_0|H_1|{\Phi }_0\rangle +\langle {\Phi }_0|H_1{ℬ}_0{H}_1\mathsf{\Omega }|{\Phi }_0\rangle$$

(34)

Equations (32) and (34) together define the infinite order Brillouin-Wigner perturbation expansion. In principle, the wave operator, Ω, is an n-particle operator and so we can write

$$\mathsf{\Omega }={\mathsf{\Omega }}^{(0)}+{\mathsf{\Omega }}^{(1)}+{\mathsf{\Omega }}^{(2)}+\dots +{\mathsf{\Omega }}^{(i)}+\dots ++{\mathsf{\Omega }}^{(n)}$$

(35)

where ${\mathsf{\Omega }}^{(n)}$ involves n particles only. Without loss of generality, we can always put ${\mathsf{\Omega }}^{(0)}=1$. Substituting (35) into (34) gives

(36)

or

$$\epsilon =\langle {\Phi }_0|H|{\Phi }_0\rangle +{\epsilon }_2^{BW}+{\Re }^{\prime }$$

(37)

Where ${\epsilon }_2^{BW}$ is the second order Brillouin-Wigner energy component

$${\epsilon }_2^{BW}=\langle {\Phi }_0|H_1{ℬ}_0{H}_1|{\Phi }_0\rangle$$

(38)

and the remainder term, ${\Re }^{\prime }$, is given by

$${\Re }^{\prime }={\displaystyle \sum _{i=1}^n\langle {\Phi }_0|H_1{ℬ}_0{H}_1{\mathsf{\Omega }}^{(i)}|{\Phi }_0\rangle }$$

(39)

Introducing the sum-over-states form of the Brillouin-Wigner resolvent (25), we have for the second order energy

$${\epsilon }_2^{BW}={\displaystyle \sum _{k\ne 0}\frac{\langle {\Phi }_0|H_1|{\Phi }_k\rangle \langle {\Phi }_k|H_1|{\Phi }_0\rangle }{\epsilon -E_k}}$$

(40)

and for the remainder

$${\epsilon }_2^{BW}={\displaystyle \sum _{i=1}^n}\left[{\displaystyle \sum _{k\ne 0}\frac{\langle {\Phi }_0|H_1|{\Phi }_k\rangle \langle {\Phi }_k|H_1{\mathsf{\Omega }}^{(i)}|{\Phi }_0\rangle }{\epsilon -E_k}}\right]$$

(41)

The matrix elements $\langle {\Phi }_k|H_1{\mathsf{\Omega }}^{(i)}|{\Phi }_0\rangle$ are obtained by using (32)

(42)

Equation (42) is the working equation for Brillouin-Wigner based configuration interaction.

In order to recover an approximation for the total energy in the form given in (22), we use the identity

$${(\epsilon -E_k)}^{-1}={(E_0-E_k)}^{-1}-{(E_0-E_k)}^{-1}\Delta E_0{(\epsilon -E_k)}^{-1}$$

(43)

to write (37) as

(44)

where ${\Re }_0$ is the Rayleigh-Schrödinger resolvent

$${\Re }_0={\displaystyle \sum _{k\ne 0}\frac{|{\Phi }_k\rangle \langle {\Phi }_k|}{E_0-E_k}}$$

(45)

The first and second terms on the left-hand-side of (44) are just the Møller-Plesset series through second order in the energy and so we have

$$\epsilon =E_0+{\epsilon }_1^{(0)}+{\epsilon }_2^{(0)}+{\Re }^{″}$$

(46)

where the remainder term ${\Re }^{″}$ is

(47)

Now we expect one- and two-body effects to be most significant and we, therefore, make the approximation

$$\mathsf{\Omega }\approx 1+{\mathsf{\Omega }}^{(1)}+{\mathsf{\Omega }}^{(2)}$$

(48)

Furthermore, if we restrict our attention to single and double replacements with respect to Φ₀ in equation (42), then equations (40), (41) and (42) provide a computational scheme which realizes the limited configuration interaction method in its CISD form. In this approximation the evalua-tion of the remainder term ${\Re }^{″}$ has roughly the same computational complexity as the MP2 energy itself. It provides information about the the validity of MP2 in cases where single and double substitutions play a dominant role.

4 Discussion

In previous work [28,29] we have shown how Brillouin-Wigner perturbation theory can be used to obtain a posteriori corrections to the method of limited configuration interaction which allow for the non-linear scaling with electron number. A posteriori corrections were obtained for both the single reference [28] and the multireference [29] case. We have also discussed the use of multireference Brillouin-Wigner methods for many-body systems [30].

In the present work, an analysis of the limited configuration interaction CISD method using the Brillouin-Wigner expansion has been used to obtain an estimate of the remainder term in MP2. The computational complexity associated with this estimate of the remainder term is roughly comparable with that of the MP2 energy expression itself.

We note that the Brillouin-Wigner coupled cluster expansion could also be used to obtain an estimate of the remainder term ${\Re }^{″}$.

References