Restoring the Pauli principle in the random phase approximation ground state

TL;DR

This paper introduces a method to correct the violation of the Pauli exclusion principle in the RPA ground state, improving its accuracy in electronic structure calculations.

Contribution

The authors propose a correction to the RPA wavefunction that restores the Pauli principle, enhancing the physical accuracy of RPA-based methods.

Findings

Corrected RPA yields molecular dipole moments closer to experimental values.

Corrected RPA provides electronic kinetic energies comparable to high-level methods.

The method improves the physical consistency of RPA calculations.

Abstract

Random phase approximation ground state contains electronic configurations where two (and more) identical electrons can occupy the same molecular spin-orbital violating the Pauli exclusion principle. This overcounting of electronic configurations happens due to quasiboson approximation in the treatment of electron-hole pair operators. We describe the method to restore the Pauli principle in the RPA wavefunction. The proposed theory is illustrated by the calculations of molecular dipole moments and electronic kinetic energies. The Hartree-Fock based RPA, which is corrected for the Pauli principle, gives the results of comparable accuracy with M{\o}ller-Plesset second order perturbation theory and coupled-cluster singles and doubles method.

Equations

$$|{\Psi }_0⟩\sim e^{\sum T_{p_1{h}_1{p}_2{h}_2}{c}_{p_1}^{†}{c}_{h_1}{c}_{p_2}^{†}{c}_{h_2}}|HF⟩.$$

$$[c_{h_2}^{†}{c}_{p_2},c_{p_1}^{†}{c}_{h_1}]\approx {\delta }_{p_1{p}_2}{\delta }_{h_1{h}_2},$$

$$\displaystyle a^{\dagger}_{k\sigma}=\left\{\begin{array}[]{l}c^{\dagger}_{k\sigma},\text{ if }k\text{ is virtual orbital},\\ c_{k\sigma},\text{ if }k\text{ is occupied orbital},\end{array}\right.$$

$$\displaystyle a_{k\sigma}=\left\{\begin{array}[]{l}c_{k\sigma},\text{ if }k\text{ is virtual orbital},\\ c^{\dagger}_{k\sigma},\text{ if }k\text{ is occupied orbital}.\end{array}\right.$$

$$a_{k\sigma }|HF⟩=0\text{ for all }k,\sigma .$$

$$C_{ph}^{†}(JM)=\sum _{\sigma {\sigma }^{\prime }}⟨\frac{1}{2}\sigma \frac{1}{2}{\sigma }^{\prime }|JM⟩a_{p\sigma }^{†}{a}_{h\overline{{\sigma }^{\prime }}}^{†}.$$

$$|{\Psi }_0⟩=N{e}^S|HF⟩,$$

$$S=\frac{1}{2}\sum _{JM}\sum _{ph{p}^{\prime }{h}^{\prime }}{T}_{ph{p}^{\prime }{h}^{\prime }}^J{C}_{ph}^{†}(JM)C_{p^{\prime }{h}^{\prime }}^{†}(\overline{JM}).$$

$$\sum _{JM}=\sum _{J=0}^1\sum _{M=-J}^J.$$

$$Q_{JMi}|{\Psi }_0⟩=0,\text{ for all }J,M,\text{ and }i,$$

$$Q_{JMi}^{†}=\sum _{ph}{X}_{ph}^{Ji}{C}_{ph}^{†}(JM)-Y_{ph}^{Ji}{C}_{ph}(\overline{JM})$$

$$\sum _{p^{\prime }{h}^{\prime }}{A}_{ph,p^{\prime }{h}^{\prime }}^J{X}_{p^{\prime }{h}^{\prime }}^{Ji}+B_{ph,p^{\prime }{h}^{\prime }}^J{Y}_{p^{\prime }{h}^{\prime }}^{Ji}={\omega }_{Ji}{X}_{ph}^{Ji},$$

$$\sum _{p^{\prime }{h}^{\prime }}{B}_{ph,p^{\prime }{h}^{\prime }}^J{X}_{p^{\prime }{h}^{\prime }}^{Ji}+A_{ph,p^{\prime }{h}^{\prime }}^J{Y}_{p^{\prime }{h}^{\prime }}^{Ji}=-{\omega }_{Ji}{Y}_{ph}^{Ji}.$$

$$A_{ph,p^{\prime }{h}^{\prime }}^J=({ϵ}_p-{ϵ}_h){\delta }_{p{p}^{\prime }}{\delta }_{h{h}^{\prime }}+\left[1+(-1{)}^J\right](ph|p^{\prime }{h}^{\prime })-(h{h}^{\prime }|p^{\prime }p),$$

$$B_{ph,p^{\prime }{h}^{\prime }}^J=\left[1+(-1{)}^J\right](ph|p^{\prime }{h}^{\prime })-(p^{\prime }h|p{h}^{\prime }),$$

$$Q_{JMi}{e}^S|HF⟩=0$$

$$e^S(Q_{JMi}+[Q_{JMi},S]+\frac{1}{2!}[[Q_{JMi},S],S]+...)|HF⟩=0.$$

$$(Q_{JMi}+[Q_{JMi},S])|HF⟩=0.$$

$$\frac{1}{2}\sum _{JM}\sum _{p_1{h}_1{p}_2{h}_2}\sum _{p_3{h}_3}{T}_{p_1{h}_1{p}_2{h}_2}^J{X}_{p_3{h}_3}^{J_{3}i}\left[C_{p_3{h}_3}(J_3{M}_3),C_{p_1{h}_1}^{†}(JM)C_{p_2{h}_2}^{†}(\overline{JM})\right]|HF⟩$$

$$=\sum _{p_3{h}_3}{Y}_{p_3{h}_3}^{J_{3}i}{C}_{p_3{h}_3}^{†}(\overline{J_3{M}_3})|HF⟩.$$

$$\displaystyle\sum_{JM}\sum_{p_{1}h_{1}p_{2}h_{2}}\sum_{p_{3}h_{3}}T^{J}_{p_{1}h_{1}p_{2}h_{2}}X^{J_{3}i}_{p_{3}h_{3}}\Big{\{}\;\delta_{JJ_{3}}\delta_{MM_{3}}\delta_{p_{1}p_{3}}\delta_{h_{1}h_{3}}$$

$$\displaystyle-\frac{1}{2}\sum_{\sigma_{3}\sigma_{3}^{\prime}\sigma_{1}\sigma_{1}^{\prime}}\langle\frac{1}{2}\sigma_{3}\frac{1}{2}\sigma_{3}|J_{3}M_{3}\rangle\langle\frac{1}{2}\sigma_{1}\frac{1}{2}\sigma_{1}^{\prime}|JM\rangle\Delta(p_{3}\sigma_{3},p_{1}\sigma_{1},h_{3}\sigma_{3}^{\prime},h_{1}\sigma_{1}^{\prime})\Big{\}}\;$$

$$\times C_{p_2{h}_2}^{†}(\overline{JM})|HF⟩=\sum _{p_3{h}_3}{Y}_{p_3{h}_3}^{J_{3}i}{C}_{p_3{h}_3}^{†}(\overline{J_3{M}_3})|HF⟩,$$

$$\Delta (p_1{\sigma }_1,p_2{\sigma }_2,h_1{\sigma }_1^{\prime },h_2{\sigma }_2^{\prime })={\delta }_{p_1{p}_2}{\delta }_{{\sigma }_1{\sigma }_2}{a}_{h_2\overline{{\sigma }_2^{\prime }}}^{†}{a}_{h_1\overline{{\sigma }_1^{\prime }}}+{\delta }_{h_1{h}_2}{\delta }_{{\sigma }_1^{\prime }{\sigma }_2^{\prime }}{a}_{p_2{\sigma }_2}^{†}{a}_{p_1{\sigma }_1}.$$

$$\left[C_{p_1{h}_1}(J_1{M}_1),C_{p_2{h}_2}^{†}(J_2{M}_2)\right]={\delta }_{p_1{p}_2}{\delta }_{h_1{h}_2}{\delta }_{J_1{J}_2}{\delta }_{M_1{M}_2}$$

$$-\underset{\text{Pauli correction}}{\underbrace{\sum _{{\sigma }_1{\sigma }_1^{\prime }{\sigma }_2{\sigma }_2^{\prime }}⟨\frac{1}{2}{\sigma }_1\frac{1}{2}{\sigma }_1^{\prime }|J_1{M}_1⟩⟨\frac{1}{2}{\sigma }_2\frac{1}{2}{\sigma }_2^{\prime }|J_2{M}_2⟩\Delta (p_1{\sigma }_1,p_2{\sigma }_2,h_1{\sigma }_1^{\prime },h_2{\sigma }_2^{\prime })}}.$$

$$\sum _{p^{\prime }{h}^{\prime }}{T}_{ph{p}^{\prime }{h}^{\prime }}^J{X}_{p^{\prime }{h}^{\prime }}^{Ji}-\underset{\text{Pauli correction}}{\underbrace{\sum _{p^{\prime }{h}^{\prime }{J}^{\prime }}{T}_{p{h}^{\prime }{p}^{\prime }h}^{J^{\prime }}W(J^{\prime },J)X_{p^{\prime }{h}^{\prime }}^{Ji}}}=Y_{ph}^{Ji}.$$

$$W(J^{\prime },J)=\sum _{M^{\prime }=-J^{\prime }}^{J^{\prime }}\sum _{{\sigma }_1{\sigma }_1^{\prime }{\sigma }_2{\sigma }_2^{\prime }}(-1{)}^{J+J^{\prime }+M+M^{\prime }}$$

$$\times ⟨\frac{1}{2}{\sigma }_2\frac{1}{2}{\sigma }_1^{\prime }|JM⟩⟨\frac{1}{2}{\sigma }_1\frac{1}{2}{\sigma }_1^{\prime }|J^{\prime }{M}^{\prime }⟩⟨\frac{1}{2}{\sigma }_2\frac{1}{2}{\sigma }_2^{\prime }|J^{\prime }-M^{\prime }⟩⟨\frac{1}{2}{\sigma }_1\frac{1}{2}{\sigma }_2^{\prime }|J-M⟩.$$

$$W(0,0)=\frac{1}{2},W(0,1)=-\frac{1}{2},$$

$$W(1,0)=-\frac{3}{2},W(1,1)=-\frac{1}{2}.$$

$${\rho }_{p^{\prime }p}=\sum _{\sigma }{a}_{p\sigma }^{†}{a}_{p^{\prime }\sigma },{\rho }_{h^{\prime }h}=\sum _{\sigma }{a}_{p\sigma }^{†}{a}_{p^{\prime }\sigma }.$$

$$[{\rho }_{p^{\prime }p},C_{p_1{h}_1}^{†}(JM)]={\delta }_{p^{\prime }{p}_1}{C}_{p{h}_1}^{†}(JM),$$

$$[{\rho }_{p^{\prime }p},S]=\sum _{JM}\sum _{h_1{p}_2{h}_2}{T}_{p^{\prime }{h}_1{p}_2{h}_2}^J{C}_{p{h}_1}^{†}(JM)C_{p_2{h}_2}^{†}(\overline{JM}).$$

Full text

1 Introduction

Random phase approximation (RPA) has recently emerged as a very promising method to study correlated ground states in molecules [1, 2, 3, 4, 5, 6, 7, 8]. RPA calculations, especially when they are performed based on the density functional theory, give accurate total molecular correlation energies, reaction barriers, electron affinities, ionisation potentials [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13] and provide a correct description of van der Waals intermolecular interactions [14, 15, 16]. The connection of RPA with the coupled-cluster doubles approach has been recently established [2, 17], which puts the method into the context of traditional quantum chemistry practices.

In this paper we would like to draw attention to the problem of violation of the Pauli exclusion principle in the RPA ground state and discuss the way of correcting it. This problem has been known in nuclear physics since 1968 [18] and since then, has been a subject of intensive studies there [19, 20, 21]. In quantum chemistry, despite early pioneering works of Simons [22], Linderberg and Ohrn [23], the issues associated with the antisymmetry of RPA wavefunction have been largely overlooked. To understand the problem, let us look at the (schematically written) RPA ground state wavefunction [24, 25]

[TABLE]

Here and throughout the paper, indices $h$ and $p$ refer to Hartree-Fock occupied and virtual molecular orbitals, respectively, and $c^{\dagger}_{k}(c_{k})$ creates (annihilates) an electron in the Hartree-Fock molecular orbital $k$. $|HF\rangle$ is the Hartree-Fock ground state vector. The spin index is suppressed in the introduction for the simplicity. RPA treats particle-hole pair creation operators $c^{\dagger}_{p}c_{h}$ as bosons, approximating the commutator (so-called quasiboson approximation[24, 25]) as

[TABLE]

and therefore neglecting the internal fermionic structure of these operators. This leads to the over-counting of configurations in the RPA ground state wavefunction. The configurations with $p_{1}=p_{2}$ and $h_{1}=h_{2}$ in the exponent (1), which otherwise would be forbidden by the Pauli principle, do contribute to the RPA ground state. The violation of the Pauli principle is amplified in the RPA due to the exponential structure of the ground state – expanding the exponent we see that terms containing two, four, six, eight and more identical electrons are allowed to occupy the same single-particle state. Therefore, the RPA ground state wavefunction (if it is treated in quasiboson approximation) over-correlates the electronic systems by including unphysical configurations and by letting the electrons to come closer in the violation of the Pauli principle.

How can we remove these unphysical states from the RPA ground state wavefunction and restore the Pauli principle? How can we correct the final expression for computed observables for the Pauli principle and at the same time preserve the original RPA equations? These are the questions which interest us in this paper.

2 Theory

2.1 RPA ground state vector with the Pauli principle corrections

We will work with the Hartree-Fock quasiparticle creation and annihilation operators, $a^{\dagger}_{k\sigma}$ and $a_{k\sigma}$, which are related to the original electron creation and annihilation operators in the Hartree-Fock spin-orbital, $c^{\dagger}_{k\sigma}$ and $c_{k\sigma}$, via the transformations [26]:

[TABLE]

and

[TABLE]

Here $\sigma=-1/2,1/2$ are two projections of electron spin. The quasiparticle operators defined above have Hartree-Fock ground state $|HF\rangle$ as a vacuum

[TABLE]

The building blocks for the RPA wavefunction are the particle-hole pair creation operators $C^{\dagger}_{ph}(JM)$ defined as

[TABLE]

It creates the particle-hole excited pair with spin $J$ and spin projection $M$. Over-line over spin indices means angular momentum time-reversal state $C_{ph}(\overline{JM})=(-1)^{J+M}C_{ph}({J-M})$ and $a^{\dagger}_{h\overline{\sigma}}=(-1)^{1/2+\sigma}a^{\dagger}_{h-\sigma}$. The Clebsch-Gordan coefficient $\langle\frac{1}{2}\sigma\frac{1}{2}\sigma^{\prime}|JM\rangle$ couples electronic spins of occupied and virtual molecular orbitals into singlet $J=0,M=0$ or triplet $J=1,M=-1,0,1$ states.

The RPA ground state wavefunction has the following form[24, 25]

[TABLE]

where

[TABLE]

Here and everywhere in the paper

[TABLE]

The correlated RPA ground state $|\Psi_{0}\rangle$ is a vacuum for the RPA excitation annihilation operators

[TABLE]

where the excitation annihilation operator is obtained from the corresponding creation operator

[TABLE]

via Hermitian conjugation.

The structure of operator $Q^{\dagger}_{JMi}$ and related RPA excitation energy $\omega_{Ji}$ are determined from the standard system of RPA equations, which we write separately for singlet ($J=0$) and triplet ($J=1$) electronic excitations

[TABLE]

Matrices $A$ and $B$ are given by the following expressions

[TABLE]

[TABLE]

where $\epsilon_{k}$ is Hartree-Fock energy of $k$ molecular orbital and $(kl|mn)$ is a two-electron integral in the Mulliken notations.

Suppose that we know $X^{Ji}_{ph}$ and $Y^{Ji}_{ph}$ in (15) from the solution of the RPA equations. Our goal in this section is to solve (14) directly and determine coefficient $T$ without the use of the quasiboson approximation (2). We rewrite (14) as

[TABLE]

Using the standard operator identity we get

[TABLE]

We neglect higher order terms $O(S^{2})$ in this expansion and the equation for operator $S$ becomes

[TABLE]

Substituting explicit expression for $Q_{JMi}$ into (22) and rearranging the terms, we get

[TABLE]

The commutator in the l.h.s. of (23) can be exactly computed and (23) becomes

[TABLE]

where

[TABLE]

To perform this derivation we use exact commutator between particle-hole pair creation and annihilation operators

[TABLE]

The commutator (26) deserves special discussion. If we kept only the first term in r.h.s (26), this would correspond to standard RPA quasiboson approximation. The Pauli correction part is roughly proportional to the total number of particle and hole excitations in the RPA state, that means the more RPA ground state deviates from Hartree-Fock Stater determinant, the larger the Pauli corrections contribution to the r.h.s of the commutator (26) become and the more pronounced deviation from the bosonic behaviour of operators $C^{\dagger}(C)$ we see.

Coming back to (24) and taking into account that $|HF\rangle$ is vacuum for the Hartree-Fock quasiparticle operators, after some straightforward algebra we transform (24) to the following form:

[TABLE]

Here matrix $W$ is

[TABLE]

Matrix $W$ (28) is calculated numerically by summing up the corresponding Clebsch-Gordan coefficients and it is

[TABLE]

The similar to (27) equation was obtained to describe RPA ground state correlations in spherical nuclei by Lenske and Wambach [20].

2.2 Calculation single-particle density matrices

Many molecular observables can be expressed in terms of single-particle density matrix. We define single-particle density matrix operators as

[TABLE]

We show detailed derivations of the expectation value for operator $\rho_{p^{\prime}p}$ in the RPA ground state. The calculations for $\rho_{h^{\prime}h}$ are almost identical, we simply write the final expression for its expectation value in the end of the derivations.

We begin with the calculations of two commutators

[TABLE]

[TABLE]

The expectation value is

[TABLE]

Taking into account that $\rho_{p^{\prime}p}|HF\rangle=0$ (since $|HF\rangle$ is vacuum for Hartree-Fock annihilation operators) and that $[[\rho_{p^{\prime}p},S],S]=0$ (since $[\rho_{p^{\prime}p},S]\sim C^{\dagger}$) so that it and all higher terms vanish, we get

[TABLE]

Then, again, since the commutator $[\rho_{p^{\prime}p},S]$ commutes with $S$, the expectation value becomes

[TABLE]

Substituting the commutator (33) into (35) we get

[TABLE]

We express $C^{\dagger}_{p_{2}h_{2}}(\overline{JM})$ in terms of RPA excitation and annihilation operators

[TABLE]

and substitute it back into equation (36):

[TABLE]

Contracting the amplitude $X$ and coefficient $T$ with the use of equation (27), we get

[TABLE]

The calculation of the first term in (38) is straightforward

[TABLE]

The calculation of the second term in (38) is more involved and is given in some details in the appendix. The result is very interesting: The second term in (38), that is the Pauli principle corrections to the expectation value of single particle density matrix $\rho_{p^{\prime}p}$, is the expectation value of the density matrix itself taken with negative sign:

[TABLE]

Combining these two terms (39,40) together into (2.2) we get

[TABLE]

which means

[TABLE]

The similar calculations are performed for single-particle density matrix for occupied molecular orbitals with the result

[TABLE]

3 Test results

We illustrate the role of the Pauli principle violation in the RPA ground state by the calculation of molecular dipole moments and electronic kinetic energies. A dipole moment is one the best observable for our purposes since it is directly related to the single-particle density matrices and it has readily available accurate experimental values for many molecules. We also computed electronic kinetic energy to assess the role of Pauli principle in the calculations of energy related quantities. The kinetic energy is chosen since it is the only part of electronic energy expressed in terms of expectation values of a single-particle operator.

Let us consider dipole moment operator written in the second quantised form

[TABLE]

where $\langle k|d|l\rangle$ is the matrix element of dipole moment between Hartree-Fock orbitals. We perform the derivations for the dipole moment, the calculations of electronic kinetic energies can be performed along the exactly the same lines. Transforming it to Hartree-Fock quasiparticle basis (5,8), we get

[TABLE]

Here the first term is the Hartree-Fock value for the molecular dipole moment. We omit in (45) terms proportional to $a^{\dagger}_{p\sigma}a^{\dagger}_{h\sigma}$ and $a_{h\sigma}a_{p\sigma}$. These terms are linearly proportional to $C^{\dagger}_{ph}(JM)$ and $C_{ph}(JM)$, respectively, therefore they are also linear in excitation creation $Q^{\dagger}_{JMi}$ and annihilation operators $Q_{JMi}$ and do not contribute to RPA ground state expectation value due to (14).

Using the expressions obtained in the previous section for single particle density matrices in the RPA ground state, we obtain the following equation for the molecular dipole moment

[TABLE]

The extra factor $1/2$ originates from the restoration of the Pauli principle in the RPA ground state.

We compute molecular dipole moments and electronic kinetic energies for a representative set of molecules (lithium hydride, hydrogen fluoride, hydrogen chloride, water, hydrogen sulfide, ammonia, methyl fluoride, and methanol). All calculations have been performed with a development version of Mendeleev computer program for ab initio quantum chemical calculations [27]. The current RPA implementation uses Hartree-Fock reference molecular orbitals and the details of the numerical implementation of RPA method is described in [28]. RPA with the Pauli principle correction calculations are compared with the standard RPA approach as well as with the Hartree-Fock (HF), Møller-Plesset second order perturbation theory (MP2) and coupled-cluster singles and doubles method (CCSD). We used 6-311++G(3df,3pd) basis set (cartesian form) in all calculations. The molecular geometries are taken from NIST Computational Chemistry Comparison and Benchmark Database [29]. The computed values of molecular dipole moments are also assessed against the experimental results.

The results of the calculations are summarised in Fig.1 and Fig.2. As seen in Fig.1, Hartree-Fock theory systematically overestimates the molecular dipole moments. The RPA (without correction for violation of the Pauli principle) produces too strongly correlated electronic ground states lowering the value of the dipole moment well below the experimental results. Restoring the Pauli principle in the RPA ground state (RPA+Pauli curve in Fig.1) brings the computed values close to the experimental data and gives the accuracy on a par with MP2 and CCSD methods. The similar behaviour is observed for electronic kinetic energies (Fig.2). The standard RPA significantly overestimates the kinetic energy while RPA with Pauli corrections yields the results close to the MP2 and CCSD values.

4 Discussions and conclusions

One of the main results of the paper are equations (27) and (41). Eq.(27) determines the operator $S$ in the RPA ground state $e^{S}|HF\rangle$ considering the exact commutation relations between particle-hole pair excitation operators (26), therefore it takes into account of the Pauli principle in construction of the ground state wavefunction. The additional Pauli correction term in (27) leads to the extra term in the single particle density matrix. The Pauli principle correction to the single particle density matrix is the same density matrix itself taken with negative sign, see (41). The proper consideration of the Pauli principle leads to extra $\frac{1}{2}$ factor in the expressions for molecular observables (42,43,46).

To summarise, we have developed theoretical method to restore the Pauli exclusion principle in the RPA correlated wavefunction. We obtained new, corrected for the Pauli principle expression for the RPA ground state correlated wavefunction. The method is a posteriori approach: The RPA eigenvalue problem remains intact and only the final expressions for molecular observables are corrected by $\frac{1}{2}$ factor. The role of the proposed Pauli principle corrections was illustrated by the calculations of molecular dipole moments and electronic kinetic energies. The corrected for the Pauli principle RPA ground gives results of comparable accuracy with Møller-Plesset second order perturbation theory and coupled-cluster singles and doubles method, while standard RPA significantly overestimates the role of electronic correlations.

Acknowledgement

We thank Alan Dzhioev for many valuable discussions.

Appendix A Derivation of eq.(40)

[TABLE]