Many-body approximations for atomic binding energies
We benchmark three standard approximations for the many-body problem – the Hartree-Fock, projected Hartree-Fock, and random phase approximations – against full numerical configuration-interaction calculations of the electronic structure of atoms, from Li through to Ne. These configuration-interaction calculations used up to uncoupled basis states, equivalent to coupled basis states (configuration state functions.) Each method uses exactly the same input, i.e., the same single-particle basis and Coulomb matrix elements, so any differences are strictly due to the approximation itself. Although it consistently overestimates the ground state binding energy, the random phase approximation has the smallest overall errors; furthermore, we suggest it may be useful as a method for efficient optimization of single-particle basis functions.
In principle we know how to find the electronic structure and binding energies of atom: solve the many-body time-independent Schrödinger equation for electrons about a nucleus with charge :
or variants thereof including relativistic corrections. In practice one cannot solve this exactly, or even fully numerically much beyond helium.
So one turns to approximations. But for any approximation there are errors intrinsic to the approximation method itself, and errors due to the input, such as truncations in the space or assumptions in the Hamiltonian. Surprisingly, rather little effort has been put into disentangling these two sources of error.
In this paper we benchmark three approximations against an ‘exact’ calculation with a fixed space and fixed nonrelativistic Hamiltonian. For our exact calculation we use full configuration-interaction ( full CI); for a chosen single-particle basis, we can find numerically converged solutions for the ground state (as well as some number of excited states). The full CI spaces are large, up to dimension in the uncoupled (‘M-scheme’) basis, or in the coupled basis of configuration state functions (CSFs); see below for further discussion. We also have a suite of codes that carry out Hartree-Fock, angular-momentum projected Hartree-Fock, and random phase approximations, using the identical input as CI. Thus our goal is to carefully document errors due the approximations themselves (and not due to the single-particle space, relativistic corrections or lack thereof, etc.); in this we are extending prior work benchmarking approximations for the nuclear many-body problem .
In the next section we define our many-body methods and discuss briefly computational issues. In Section 3 we give our results for several different spaces for second row atoms, from Li to Ne, comparing not only the ground state binding energies but also the experimentally more relevant ionization potentials and electron affinities. Because full CI is very time-consuming and our approximations computationally cheap, in Section 4 we consider the possibility of using approximations for fast optimization of the the single-particle basis, and of the considered methods find RPA the most accurate.
2Many-body methods and inputs
We work entirely in occupation or configuration space, meaning we start with some set of orthonormal single-particle states, , with good angular momentum and parity, and introduce in second quantization  fermion creation and annihilation operators that create or remove a particle from the states . For definitiveness we have single-particle states. The Hamiltonian is then 
The one-body and two-body matrix elements, (which includes both kinetic energy and the nuclear central potential) and respectively, are appropriate integrals over the Hamiltonian and the single-particle states. Because the Hamiltonian is an angular momentum scalar, we generate the matrix elements in coupled form. Given the single-particle space, the matrix elements are computed externally to all our codes and read in via a file. Any radial form of the single-particle states and any form of the interaction, including non-local interactions, are allowed and make no practical difference to our many-body codes. Because we work in second quantization, the matrix elements are antisymmeterized and we do not separate direct from exchange terms.
For this paper we restrict the Hamiltonian to the simple nonrelativistic Coulomb interaction (Equation 1). We build our single-particle radial wavefunctions from Slater-type orbitals (STOs), that is, functions of the form , which allows us to compare to other work as well as making the integrals analytic. We orthonormalize our single-particle basis so that the one-body part of (Equation 2) is diagonal. For our specific choices of parameters for our basis, we use the CVB bases sets , described below. Finally, we describe the atom as a two-species systems, with spin-up and spin-down electron occupying distinct orbits. We fix in all calculations, as we have no spin-orbit coupling.
These methods were chosen because we have existing codes for them, and thus leave for future work other methods such as the multi-configurational Hartree-Fock approximation  and many-body perturbation theory 
For full configuration-interaction (full CI) calculations , we use the BIGSTICK code , developed originally for nuclear CI calculations but which can be applied as well to atomic calculations. The only assumption is that the single-particle basis be states of good angular momentum and parity.
BIGSTICK creates a many-body basis of Slater determinants , constructed from the single-particle basis:
where is the number of particles (occupied states). Furthermore, because we have no spin-orbit coupling (a generalization we leave to future work), we fix the number of spin-up and spin-down electrons and thus fix . Because single-particle states each have good , it is almost trivial to generate many-body basis states with fixed total ; in nuclear CI calculations this is called an -scheme basis. And because the Hamiltonian is an angular momentum scalar, the final eigenstates thus will automatically have good total (and ). Although we have the capability for various truncations on the many-body basis, such as allowing only single- or double-electron excitations, we make no such truncations here as they have no correspondence for our approximation methods.
Most atomic CI programs such as CIV3  and GRASP92  use a basis of configuration state functions (CSFs) coupled up to good and (or total , hence in nuclear structure physics this is called a -scheme basis). Such a basis is smaller in dimension, by a factor of 10 to 100, but there are always trade-offs: the basis states are linear combinations of -scheme Slater determinants and hence the many-body Hamiltonian matrix elements are more complicated to calculate. And although the dimensions of an uncoupled -scheme basis are larger, the Hamiltonian matrix is much sparser. Furthermore the simplicity of the uncoupled basis allows us to avoid storing the many-body matrix element, and instead recompute them efficiently on the fly , which further reduces the storage of the Hamiltonian by a factor of 10 to 100. We can go up to uncoupled basis dimensions of on a desktop machine (and similar programs on parallel computers reach up to ), as listed in Table II, while the dimension of the largest equivalent coupled CSF ground state dimension is about , as listed in Table III.
Given these arguments, it is interesting to compare basis sizes. These are not always given in the literature, and many ‘large-scale’ calculations rely upon careful pruning of the configuration basis ; nevertheless, for comparison purposes a large scale relativistic CI calculation  used up to 197,426 basis CSFs, while a more recent nonrelativistic calculation  used up to 901,816 CSFs. Quantum chemists surpassed the one-billion determinant limit some time ago .
In a given uncoupled basis, BIGSTICK computes the many-body Hamiltonian matrix elements and finds the low-lying eigenstates, including the ground state, using the Lanczos algorithm  (the Davidson-Liu algorithm  is efficient for the diagonal-dominated atomic problem but not the more general many-body problem such as nuclear CI; furthermore recomputing the matrix elements on the fly is better suited for the Lanczos algorithm). Convergence of, say, five states take a few hundred Lanczos iterations; on a desktop machine with eight OpenMP threads we can accomplish this for an -scheme basis in a few hours or at most a few days.
With the single-particle space and Hamiltonian matrix elements fixed, the resulting eigenenergies for the full-space CI calculation are numerically exact. All other calculations described in this paper are numerical approximations to these results.
CI calculations have been previously carried out for second row elements, see for example nonrelativistic  and relativistic calculations  , though most of these were not full CI and thus were in much smaller dimensions.
and creates the trial Slater determinant
Then one finds a Slater determinant that minimizes the energy, that is, that minimizes
Because our input matrix elements are antisymmeterized, we fully include both direct and exchange terms. Working in occupation space the exchange term causes us no difficulty (or, to put it another way, any difficulty is off-shored into calculation of the antisymmeterized integrals).
Many Hartree-Fock calculations enforce good angular momentum, for closed-shell systems or closed-shell plus or minus one particle. Our HF code allows the Slater determinant to break rotational invariance, even for closed-shell systems; such Slater determinants are ‘deformed.’ Experience in nuclear physics suggest this to be a fruitful path. In a few cases, when possible, we will compare and contrast spherical and deformed HF solutions.
In order to minimize, we use the standard Hartree-Fock equations, which consist of iteratively solving
where the Hartree-Fock effective one-body Hamiltonian is
Internally we use an rectangular matrix , with the columns representing the transformed occupied states . In that case .
(Because we leave out spin-orbit coupling and conserve and , we use separate Slater determinants for spin-up and spin-down electons; this is a straightforward generalization. )
The Hartree-Fock energy is now
where we have to subtract off the potential to keep from double-counting in the single-particle energies.
2.3Projected Hartree-Fock approximation
If one does not have a closed shell, a single Slater determinant will not in general have good angular momentum. Conversely, an open-shell state with good angular momentum, which can be written as a sum of Slater determinants, will have important correlations that lower the energy.
One can start with single-particle states which have good angular momentum and then construct many-body states with good total angular momentum; these are called configuration state functions (CSFs), and have been applied to second-row elements, e.g. [?].
We take a different approach, however. Our Hartree-Fock code allows for arbitrary deformation: the single-particle states generated by the self-consistent field do not have good angular momentum. The Slater determinant built from these states is called a ‘deformed’ state. In order to project out good angular momentum, we introduce the standard projection operator 
where the Euler angles and , is the Wigner -matrix , and is the rotation operator. To compute and , we use the matrix representation
where and are matrices representing our original and rotated Slater determinants, respectively, and is a square matrix; in our single-particle basis the matrix elements are straightforward:
where is the Wigner -matrix . It is useful to note that the for the rotational matrix are those of the single-particle space, while the for the projection operator (Equation 5) are those of the many-body space.
Matrix elements between two oblique Slater determinants is relatively straightforward . First,
Calculation of the Hamiltonian (and other) expectation value requires the density matrix,
Accounting for two species is straightforward. Note that even if our original Slater determinants were real, by rotating over all angle they can become complex.
We only project out good orbital angular momentum . We could in principle project out exact spin ; we leave such a modification for future work. Because we leave off spin-orbit coupling, however, we automatically get states with good already.
We now can compute the Hamiltonian and “norm” matrix elements,
which are both Hermitian. Then we solve the generalized eigenvalue equation
where we allow to be complex. Although the deformed Hartree-Fock state will have an orientation, the final result will be independent of orientation; we confirmed this by arbitrarily rotating our HF state.
Both the Hartree-Fock and the projected Hartree-Fock energies will be variational upper bounds to the exact CI energy, with the projected Hartree-Fock energy lower than (or equal to, in the HF state has good rotational symmetry) .
2.4Random phase approximation
Once one has a Hartree-Fock calculation, it is straightforward to go on to the random phase approximation, using the matrix formulation :
Here label occupied single-particle levels in the Hartree-Fock state, while label unoccupied levels.
From this we can compute the RPA correlation energy, which is a correction to the HF energy. When working with a rotationally invariant Hamiltonian, as we do here, it is possible to have a HF state which breaks rotational symmetry; this means one could rotate the HF state in any direction and not change the HF energy. In RPA this shows up as a zero-energy mode, which must be dealt with carefully when computing the binding energy . The final result is
The RPA energy, while lower than the HF energy, is not variational. In previous calculations for nuclei in valence spaces  RPA both under- and overestimated the binding energy; in our chosen atomic cases below, RPA consistently overestimates the binding energy. Because we start from a deformed HF state, our RPA calculations do not have good angular momentum. (RPA is often said to ‘restore’ rotational symmetry, but this is not restoration of good angular momentum quantum numbers. Rather, the Hartree-Fock energy is invariant under rotation, and this invariance shows up as a mode that appears in RPA but not in the simpler Tamm-Dancoff approximation or TDA . In this sense RPA ‘restores’ rotational invariance that was broken in TDA.)
2.5Input spaces and interaction
Our codes can use (up to memory limitations, of course) any finite set of single-particle states which have good angular momentum. For convenience we use the CVB1, CVB2, and CVB3 Slater-type orbitals (STOs) and parameter values of Ema et al. . These are radial wavefunctions of the form , with the set of STOs and the values of optimized for Li through Ne. The single-particle spaces ranged from 8 orbits (for CVB1-Li) up to 22 (for CVB3-Ne).
To construct our many-electron basis states, we fixed spin-up electrons and spin-down electrons, with or , which allowed us to regain states of all possible total . We could and did choose in order to compute ionization potentials and electron affinities. We also fixed total , which allowed us to regain states of all possible total . We otherwise had no constraints, e.g., we did not restrict ourselves to single, double, or otherwise excitations; although our CI code can make such restrictions, the HF, PHF, and RPA cannot. Dimensions of the uncoupled many-body CI spaces we used are in Table ?, while the equivalent dimensions for coupled (CSF) bases for the ground states are in Table ?, where we indicated the ground state quantum numbers by .
To construct our interaction file, we first orthornormalized the STO basis. We did this by diagonalizing the one-body part of the Hamiltonian (and so this depended on the charge of the nucleus); then the one-body Hamiltonian is diagonal and expressed in terms of single-particle energies. The two-body matrix elements were first computed with the STOs, which can be done semi-analytically, and then transformed to the orthonormalized basis. Our matrix elements and resulting atomic energies for small cases were validated against an prior, independent code (M. Bromley, private communication).
The following sections give descriptions of the results obtained as well as tables and example figures. The computed ground state energies for HF, PHF, and RPA are compared to CI. The ionization potentials and electron affinities are compared to experimental values .
3.1Ground State Energy
The computed ground state energies for the neutral atoms studied are in Table ?. CI is exact for purposes of comparison. As required by variational theory, Hartree-Fock gives an upper bound to the ground state energy and projected Hartree-Fock improves upon it. The random phase approximation computes lower energies than CI in all of the calculations we performed. There are no CI calculations for several the largest basis sets in combination with the largest atoms because the computational requirements exceeded our available resources; keep in mind that our ‘approximate’ results were computed in full CI spaces, as we allow arbitrary deformation of the single-particle states.
Figure ? shows the results for carbon in graphical form; the results for all other atoms are qualitatively similar, although the error for all approximations increases as we add more electrons. For example, the error in RPA is only a few millihartrees for Li, while it is 100 millihartrees for Ne.
Note there is no difference between the HF and PHF results for Li and for Ne; this is because the HF states for neutral Li and Ne are already in states of good total L; projection cannot change this.
Table ? shows the single ionization potentials for all methods as well as experimental values  for all atoms studied. Figure ? shows the ionization potential for carbon graphically.
Because the ionization potential (and in the next section, the electron affinity) is the difference between two ground state energies, the nice clean trends of the previous section are lost. CI agrees well with experiment, within a few millihartree through C and less than 20 millihartree difference through Ne. Increasing basis size brings better agreement with experiment, though of course we note we have left out relativistic corrections.
HF consistently underestimates the ionization potentials, with PHF closer for Li through N (with HF and PHF identical for Li as the neutral Li and Li I HF states both have good angular momenta) but farther away for O, F, and Ne.
RPA results are the sporadic, sometimes better than HF / PHF and sometimes; again, this is because errors in the approximation do not cancel but exacerbate. Increasing the basis size does not always lead to better estimates; they get worse for O, F, and Ne.
Not all of the second row atoms can bind an additional electron; we only present those cases which do so experimentally: lithium, boron, carbon, oxygen and fluorine. Table ? shows the electron affinities computed by each method and the respective experimental values . Our results are comparable to other work in the field . Figure ? shows the electron affinity of carbon in graphical form.
The quality of agreement is generally poor, even for CI. In retrospect this is not surprising: we have already seen errors for all methods grow rapidly with the number of electrons, and given the smallness of second row electron affinities the results are unsurprising if disappointing. Many of the affinities even have the wrong sign. Adding an electron to the atom increases the runtimes and memory requirements over the neutral atom calculations, therefore, only the smallest basis sets for all but lithium could be used. Without the larger calculations, a clear trend is difficult to establish. It appears that the CVB1 basis set is fairly poor when computing electron affinities, sometimes giving negative results. However, based on the CVB2 calculations for boron and carbon, the CI electron affinities might converge on the experimental values as basis size increases, similar to what is seen with the ground state energy and ionization potential calculations.
All three of our approximations perform poorly, trending neither towards the CI result nor towards experiment. HF and PHF are unable to provide consistent results across either atoms or basis sets and generally tend to be too far from experimental values to be useful; nor do they show consistent trends as the basis size is increased. RPA performs better, but our results do not suggest convergence with basis size. Because it is not clear if this is a fault of the approximation or of the basis for CI (which again is our numerically exact benchmark) it would be good, if and when possible, to continue this investigation with even larger bases.
3.4Optimization of the space
We utilized the CVB STO-basis sets , which had been arrived at through optimization using both Hartree-Fock and CI calculations . But CI calculations, particularly in large dimension spaces, are tedious and time-consuming. Furthermore, a basis optimized for Hartree-Fock calculations may not be optimized for CI calculations. Therefore it might be useful to see how well our approximate methods can work as a proxy for full CI when changing the basis. In other words, if we vary the basis, does a minimum in CI have a corresponding minimum in one or more of our approximate methods?
We carried out an example optimization by altering all of the parameter values in the STO for carbon CVB1 by a scaling factor ranging from to . Table ? gives the results of a quadratic fit to each computational method. The data show that RPA has a similar minimum to CI. Figure ? shows these results graphically.
This optimization suggests that RPA is the most desirable method of those studied for basis set optimization. It is much faster than CI and clearly tracks the ground state energy as basis parameters are altered. PHF is a more useful upper bound to the ground state energy than HF, but cannot track the CI energy well enough to be used as a proxy for CI.
4Conclusions and Future work
This paper is an extension of previous work comparing the random phase approximation against CI calculations in a CI basis , testing how good is an approximation against a numerically exact result. By comparing against a numerically exact result, we eliminate any errors due to choice of model space, interaction, etc. As part of this process we find eigenvalues of very large CI matrices. To simplify matters we used an uncoupled basis, which allowed us to recomputed the many-body matrix elements on the fly in a highly efficient manner.
Not surprisingly, we found in absolute error RPA gave the best approximation to the ground state energy. RPA is not variationally bounded, and in the cases we chose it systematically overestimated the binding energy; this contrasts from previous, similar work for nuclei, where RPA both over- and underestimated the binding energy . Because the overestimate was not uniform, HF and PHF calculations sometimes gave accidentally better estimates of the ionization potentials and electron affinities.
Perhaps most usefully, we found that for an example fixed atomic system, RPA tracks the CI binding as we varied parameters for the basis, suggesting RPA might provide an efficient proxy to speed up optimization of a basis. This will be explored further in future work.
In principle we can also compute and compare transitions and strength functions, as has been done for nuclei , and static and dynamic observables such as static moments, polarizabilities, etc. This we also leave for future work.
We thank Michael Bromley for many valuable conversations and suggestions, and in particular for aid in validating our codes. The U.S. Department of Energy supported this investigation through grants DE-FG02-96ER40985 and DE-FC02-07ER41457.
- I. Stetcu and C.W. Johnson, Phys. Rev. C 66, 034301 (2002).
- I. Stetcu and C.W. Johnson, Phys. Rev. C 67, 044315 (2003); I. Stetcu and C.W. Johnson, Phys. Rev. C 69, 024311 (2004).
- J. T. Staker and C. W. Johnson, to be published.
- I. Lindgren and J. Morrison, Atomic many-body theory, 2nd. ed. (Springer-Verlag, Berlin, 1985).
- P.J. Brussard and P.W.M. Glaudemans, Shell-model applications in nuclear spectroscopy (North-Holland Publishing Company, Amsterdam 1977).
- P. Ring and P. Shuck, The nuclear many-body problem, 1st edition, Springer-Verlag, New York 1980.
- I. Ema et al., J. Comput. Chem. 24, 859 (2003).
- D. B. Cook, Handbook of computational quantum chemistry, (Oxford University Press, Oxford, 1998).
- C. Froese Fischer, T. Brage and P. Jönsson, Computational atomic structure: an MCHF approach (Institute of Physics Publishing, Bristol, 1997).
- F. Jensen, Introduction to computational chemistry, 2nd ed. (John Wiley Sons Tld, West Sussex, 2007).
- P.-O. Löwdin, Phys. Rev. 97, 1474 (1955)
- C.D. Sherrill and H.F. Schaefer III, Adv. Quant. Chem. 34, 143 (1999).
- B. A. Brown and B. H. Wildenthal, Annu. Rev. Nucl. Part. Sci. 38, 29 (1988).
- I. Shavitt, Mol. Phys. 94, 3 (1998).
- A. Hibbert, Comp. Phys. Comm. 3, 141 (1975); A. Hibbert, Comp. Phys. Comm. 35, C-298 (1984).
- F.A. Parpia, C. Froese Fischer, and I. P. Grant, Comp. Phys. Comm. 94, 249 (1996); F.A. Parpia, C. Froese Fischer, and I. P. Grant, Comp. Phys. Comm. 175, 745 (2006).
- W. E. Ormand and C. W. Johnson, to be published.
- E. Caurier and F. Nowacki, Acta Phys. Pol. B 30 (1999) 705
- N. C. Deb and A. Hibbert, At. Data and Nucl. Data Tables, 95, 184 (2009).
- G. Corrégé and A. Hibbert, At. Data and Nucl. Data Tables, 86, 19 (2004).
- G. P. Gupta and A. Z. Msezane, At. Data and Nucl. Data Tables, 89, 1 (2005).
- M. H. Chen, K. T. Cheng, and W. R. Johnson, Phys. Rev. A 64, 042507 (2001).
- M. W. J. Bromley and J. Mitroy, Phys. Rev. A 81, 052708 (2010).
- J. Olsen, P. Jørgensen, and J. Simons, Chem. Phys.Lett. 169, 463 (1990).
- R. R. Whitehead, A. Watt, B. J. Cole, and I. Morrison, Advances in Nuclear Physics (Plenum, New York, 1977), Vol. 9, p. 123.
- E. R. Davidson, Comput. Phys. 7, 519 (1993).
- M. L. Leininger, C. D. Sherrill, W. D. Allen, and H. F. Schaefer III, J. Comp. Chem. 22, 1574 (2001).
- A. W. Weiss, Phys. Rev. 122 (6), 1826-1836 (1961).
- C. F. Bunge, Phys. Rev. 168, 92 (1968)
- F. Saskai and M. Yoshimine, Phys. Rev. A 9 (1), 17-25 (1974).
- C. F. Bunge, Phys. Rev. A 14 , 1965 (1976); At. Data Nucl. Data Tables 18, 126 (1976).
- D. Silver, S. Wilson, and C. F. Bunge, Phys. Rev. A 19 (4), 1375-1382 (1975).
- O. Jitrik and C. F. Bunge, Phys. Rev. A 56 (4), 2614-2623 (1997).
- C. F. Bunge, J. Chem. Phys. 125, 014107 (2006).
- C. X. Almora-Diaz, and C. F. Bunge, Int. J. Quantum Chem. 110 (15), 2982-2988 (2010).
- P.-O. Löwdin, Phys. Rev. 97, 1490 (1955)
- S. Lunell, Phys. Rev. 173, 86 (1968).
- B. Laskowski and S. Lunell, Int. J. Quantum Chem. 9 (S9), 175-182 (1975).
- R. D. Lawson, Theory of the nuclear shell model (Clarendon Press, Oxford, 1980).
- A. R. Edmonds, Angular momentum in quantum mechanics (Princeton University Press, New Jersey, 1996).
- G. H. Lang, C. W. Johnson, S. E. Koonin, and W. E. Ormand, Phys. Rev. C 48, 1518 (1993); E. Y. Loh, Jr. and J. E. Gubernatis, in Electronic Phase Transitions, edited by W. Hanke and Yu. V. Kopaev (Elsevier Science Publishers B.V., New York, 1992).
- David R. Lide, CRC Handbook of Chemistry and Physics 85th edition, (CRC Press Inc, Florida 2005).
- I. Stetcu, PhD. thesis, Louisiana State University 2003 (unpublished).
- H. Hotop and W. C. Lineberger, J. Phys. Chem. Ref. Data 4 (3), 539-576 (1975).
- H. Jiang and E. Engel, J. Chem. Phys. 127, 184108 (2007).
- A. Grüneis et al., J. Chem. Phys. 131, 154115 (2009).
- L. A. Cole and J. F. Perdew, Phys Rev A 25, 1265 (1982).