Atoms in Generalized Orbital Configurations: Towards Atom-Dedicated Density Functionals

Abstract

Deriving a practical formula for the atomic body with generalized shell occupations, we perform a detective analysis of the radial distribution in the exchange energy, hinting at ideas about new types of density functionals, dedicated to the specifics of the electronic structure of atoms, exploiting the intrinsic spherical symmetry.

1. Introduction

The electronic structure of the atoms is tacitly perceived nowadays as something definitively resolved. Aiming at immediate applications, from the start of the era of low power computers, the quantum chemists leaned toward rather drastic approximations and numerical compromises. This trend is continued even nowadays, in the quest for computing larger structures, challenging the domain of nano-chemistry. In this race, feedback with a sound account of the atomic structure is somewhat ignored.

A clear example is the extensive and intensive use of Gaussian type orbitals (GTOs) [1,2]. The users are aware that the r^k·exp(−ζ·r²) primitives are not the best choices, since the known exact solutions for the hydrogen atom suggests that the Slater type orbitals (STOs) [3], with r^k·exp(−ζ·r) components, are more natural options. The GTO-based codes [4,5] are more massively employed than the rather rare STO-type ones [6]. In addition, as we observed recently [7], the GTO computational establishment carries, since early editions of computer codes, a hidden drawback, overlooked by most users. Namely, the GTO codes and repositories are built on drastically limited radial power-patterns of r^k·exp(−ζ·r²) primitives. Thus, the s-type orbitals are always without radial prefactors (i.e., k = 0), the p-type ones get only k = 1, while the d functions correspond to k = 2, and so on. An inspiration borrowed from hydrogen prototypic exact solution will suggest that a larger variety of k powers for each type of shell is conceptually welcomed and technically beneficial. Another collateral loss induced by the success of GTO technology comes from the evaluation of all the integrals, including the atomic ones, by products dichotomized with respect to x, y, and z Cartesian coordinates, because this is technically convenient [8].

However, in this way, a whole wisdom of factorizing the atomic problems in radial and angular parts is lost. The black-box treatment of the atomic integral ignores the advantage of thinking of the inter-electronic effects (Coulomb, exchange, and correlation) with the help of the so-called Slater–Condon integrals [9]:

$$F_{ab}^k\equiv F_{n_a{l}_a,n_b{l}_b}^k={{\displaystyle \int }}_{r_1=0}^{\infty }{{\displaystyle \int }}_{r_2=0}^{\infty }{R}_{n_a{l}_a}{\left(r_1\right)}^2{R}_{n_b{l}_b}{\left(r_2\right)}^2\frac{min{(r_1,r_2)}^k}{max{(r_1,r_2)}^{k+1}}{r}_1^2{r}_2^{2}d{r}_{1}d{r}_2$$

(1)

$$G_{ab}^k\equiv G_{n_a{l}_a,n_b{l}_b}^k={{\displaystyle \int }}_{r_1=0}^{\infty }{{\displaystyle \int }}_{r_2=0}^{\infty }{R}_{n_a{l}_a}\left(r_1\right)R_{n_b{l}_b}\left(r_1\right)R_{n_a{l}_a}\left(r_2\right)R_{n_b{l}_b}\left(r_2\right)\frac{min{(r_1,r_2)}^k}{max{(r_1,r_2)}^{k+1}}{r}_1^2{r}_2^{2}d{r}_{1}d{r}_2$$

(2)

defined as a function of the radial components R(r) of atomic functions belonging to a given couple of shells, a and b, each identified with their main and secondary quantum numbers, n and l. In the body of this article we will perform analyses of the exchange density in atoms, rationalizing the energy with the help of Slater–Condon formulas.

The exchange energy in atoms deserves continuous attention. The importance of which is weighted by the high popularity gained by the application of density functional theory (DFT) in practical calculations [10]. The sensible issue is that practically the whole applied DFT stays on extrapolations originating from the uniform electron gas [11]. The many alleviations of this rather limiting hypothesis, belonging to the classes of generalized gradient approximation (GGA) [12,13] or hybrid functionals [14,15], are tacitly an offshoot of the local density approximation. In the atom, the variation of radial density is very steep, decaying exponentially from the height of the nuclear cusp. In such circumstances, the foundation based on electron gas can be questioned, aiming for new paradigms, and exploiting the regularities due to the spherical symmetry. We will attempt this idea in the actual work.

2. Results and Discussion

2.1. The Energy of the Atom with a General Shell Occupation Scheme

As a basic tool of investigation in the aimed area, we need an explicit formula, with the function of arbitrary occupations of atomic shells, compatible with calculations performed by state-averaged complete active space (CAS), or, in certain cases, in the restricted Hartree–Fock (RHF) frame. Basic textbooks [16] offer many details of the atomic structure theory, which, however, are not handy and general enough to be applied as ancillary tools, compliant with the modern numeric computational procedures. Then, we will tailor a closed formula for the total atomic energy as a function of shell occupations.

The derivation can be started with master formulas met in the so-called DFT + U methods [17,18], where the two-electron energy of a given configuration, defined by its α and β occupation numbers, is written as:

$$V_{ee}[{\left\{n\right\}}_{\alpha },{\left\{n\right\}}_{\beta }]=\frac{1}{2}{\displaystyle \sum }_a^{a\in \alpha }{\displaystyle \sum }_{a\prime \ne a}^{a\in \alpha }{n}_a{n}_{a\prime }\left(U_{aa\prime } - J_{aa\prime }\right) + \frac{1}{2}{\displaystyle \sum }_b^{b\in \beta }{\displaystyle \sum }_{b\prime \ne b}^{b\prime \in \beta }{n}_b{n}_{b\prime }\left(U_{bb\prime } - J_{bb″}\right) + {\displaystyle \sum }_a^{a\in \alpha }{\displaystyle \sum }_b^{b\in \beta }{n}_a{n}_b{U}_{ab}$$

(3)

where U refers to Coulomb-type integrals, actually equivalent to the F_ab⁰ Slater–Condon parameters, and J to the averaged exchange integrals, taken as combinations of F^k or G^k integrals.

The above formulation, ascribed in the spin-unrestricted spirit, does not offer the final answer, being valid for a given S_z = S spin projection resulted from the balance of α and β populations, but not for the S spin itself. Rethinking the situation as a spin-restricted case, the electron–electron potential for a definite S spin state can be obtained as the difference between the above generic formula taken at S_z = S and those for the S_z = S − 1 projection, averaging over the whole set of configurations that span these S_z quantum numbers.

The obtained general formula is given in Equation (4), admitting general occupations, $p_{l_i}$, over the whole set of shells, l denoting the secondary quantum number, while i counts the repetition of the given shell type:

$$\begin{gathered} E={\displaystyle \sum }_l{\displaystyle \sum }_i{p}_{l_i}{h}_{l_i}+ {\displaystyle \sum }_l{\displaystyle \sum }_i\left\{\frac{1}{2}{p}_{l_i}\left(p_{l_i} - 1\right)F_{l_i{l}_i}^0 + \left[{\alpha }_l\cdot p_{l_i} + {\beta }_{ l}\cdot p_{l_i}^2 + {\gamma }_l\cdot {\sigma }_{l_i}\left({\sigma }_{l_i} + 1\right)\right]J_{l_i{l}_i }\right\} \\ + {\displaystyle \sum }_{l\le l\prime }{\displaystyle \sum }_{l_i\ne l{\prime }_{i\prime }}\left\{p_{l_i}{p}_{l{\prime }_{i\prime }}\left(F_{l_{i}l{\prime }_{i\prime }}^0 - \frac{1}{2}{J}_{l_{i}l{\prime }_{i\prime }}\right) - 2{\sigma }_{l_i}{\sigma }_{l{\prime }_{i^{\prime }}}\cdot J_{l_{i}l{\prime }_{i\prime }}\right\} \end{gathered}$$

(4)

The $h_{l_i}$ elements are the orbital one-electron energies (kinetic and electron-nuclear parts), the $F_{l_i{l}_i}^0$ and $F_{l_{i}l{\prime }_{i\prime }}^0$ are the respective intra- and inter-shell Coulomb integrals, while the $J_{l_i{l}_i}$ and $J_{l_{i}l{\prime }_{i\prime }}$ are the corresponding averaged exchange integrals. For a lⁿ sub-configuration, having $n \equiv p_{l_i}$ electrons in the l_i shell, more exactly with n_α spin-up and n_β spin-down particles, the ${\sigma }_{l_i}$ = (n_α − n_β)/2 denotes the net spin of the shell. The shell-distributed spin quantities are summed to the total spin of the atom $S = {{\displaystyle \sum }}_l{{\displaystyle \sum }}_i{\sigma }_{l_i}$. Note that, according to the previous discussion, the obtained energy refers to S as a good quantum number. The tacit assumption is that the spin from possible multiple open shells of the general configuration is coupled all parallel. This condition induces a certain limitation, but yet the situation is sufficiently general.

The coefficients ascribed to the intra-shell exchange were obtained by induction, after analyzing the results of spherical averaging on the all possible sets of integer occupations of the shells. Their expressions were found as follows:

$${\alpha }_l=\frac{4l+3}{4\left(l+1\right)} ,\hspace{1em}{\beta }_l=-\frac{2l+3}{8\left(l+1\right)} ,\hspace{1em}{\gamma }_l=\frac{2l+1}{2\left(l+1\right)} ,$$

(5)

if l > 0. For the s-type orbitals, there is no intra-shell exchange (all the above factors being null).

Table 1. Generic formulas for inter-shell averaged exchange integrals.

Shell Couples	s	p	d	f
s	$J_{s_1{s}_2}=G_{n_{1}s{n}_2}^0$
p	$J_{ps}=J_{sp}=\frac{1}{3}{G}_{sp}^1$	$J_{p_1{p}_2}=\frac{1}{3}{G}_{n_{1}p{n}_{2}p}^0 + \frac{2}{15}{G}_{n_{1}p{n}_{2}p}^2$
d	$J_{ds}=J_{sd}=\frac{1}{5}{G}_{sd}^2$	$J_{dp}=J_{pd}=\frac{2}{15}{G}_{dp}^1 + \frac{3}{35}{G}_{dp}^3$	$J_{d_1{d}_2}=\frac{1}{5}{G}_{n_{1}d{n}_{2}d}^0 + \frac{2}{35}{G}_{n_{1}d{n}_{2}d}^2 + \frac{2}{35}{G}_{n_{1}d{n}_{2}d}^4$
f	$J_{fs}=J_{sf}=\frac{1}{7}{G}_{sf}^3$	$J_{fp}=J_{pf}=\frac{3}{35}{G}_{fp}^2 + \frac{4}{63}{G}_{fp}^4$	$J_{fd}=J_{df}=\frac{3}{35}{G}_{fd}^1 + \frac{4}{105}{G}_{fd}^3 + \frac{10}{231}{G}_{fd}^5$	$J_{f_1{f}_2}=\frac{1}{5}{G}_{n_{1}f{n}_{2}f}^0 + \frac{2}{35}{G}_{n_{1}f{n}_{2}f}^2 + \frac{2}{77}{G}_{n_{1}f{n}_{2}f}^4 + \frac{100}{3003}{G}_{n_{1}f{n}_{2}f}^6$

shows the formulas for the inter-shell exchange integrals. The diagonal contains the same type of shells, but spanning different functions.

The J_ll are the average values of the exchange integrals over all the orbital couples within a given shell. For the respective p, d, and f cases, these quantities are

$$J_{pp}=\frac{1}{5}{F}_{pp}^2$$

(6a)

$$J_{dd}=\frac{1}{14}{F}_{dd}^2 + \frac{1}{14}{F}_{dd}^4$$

(6b)

$$J_{ff}=\frac{2}{45}{F}_{ff}^2 + \frac{1}{33}{F}_{ff}^4 + \frac{50}{1287}{F}_{ff}^6$$

(6c)

as function of Slater–Condon radial parameters.

To the best of our knowledge, a closed formula for the atomic energy of the atomic body, as a function of general shell occupation numbers, is not presented in the specialized literature, particularly in the concern of the intra-shell exchange terms fulfilling the meaning of state-averaged CAS over multiple open-shell configurations.

We realized our own code for atomic calculations based on GTO primitives and explicit use of Slater–Condon integrals. The Formula (4) was checked to yield the same results with state-averaged or single-configuration ab initio calculations. It must be emphasized that retrieving the ab initio result with the above formula is not trivial. First of all, this validates the treatment. Then, recall the above sentence, that the actual quantum chemistry codes are not working with Slater–Condon parameterization, because of the lost interest in radial-angular factorization. However, equating to spherically adapted proper parameters brings transparency in the causal engine and clears the way for new rationales.

2.2. A Methodology for the Decomposition of Two-Electron Terms

Before going to the analysis step, to describe the exchange energy, further adequate tools must be prepared. In this view, we will proceed to the decomposition of two-electron terms, conducting their evaluation on numeric grids. Note that we are able to, and we do the full analytic calculation of two-electron integrals, repeating the procedure in the numeric mode only for detective purposes. The numeric integration is better done with a grid with an exponential spacing of the points. Such an integration scheme, due to Weber et al. [19], uses the following definition for the m-th point:

$$r_m=\delta r_0\frac{exp\left(mh\right) - 1}{exp\left(h\right) - 1}$$

(7)

depending on the spacing of the first point with respect of the origin, δr₀, and a scale parameter, h. The same grid is used in several programs conducting numeric calculations of the atoms, to prepare pseudo-potentials for plane-wave calculations, as is the case of the ATOMPAW module [20] from the ABINIT suite [21].

Deciding a maximal radial extension of the set, r_max, the number of points results as

$$m_{max}=\mathrm{int} \left[\frac{1}{h}\ln \left(1 + r_{max}\frac{exp\left(h\right) - 1}{\delta r_0}\right)\right]$$

(8)

The square brackets meaning the integer value. Conversely, if a certain number of points over a given radial interval is desired, the h parameter must be fitted, correspondingly. The points are associated with a set of weights:

$$w_m=h\cdot \delta r_0\frac{exp\left(kh\right)}{exp\left(h\right) - 1}$$

(9)

so that the numeric integration of a given function f(r) can be formulated as a weighted sum over the grid:

$${{\displaystyle \int }}_{r=0}^{\infty }f\left(r\right)dr\approx {\displaystyle \sum }_{m=1}^{m_{max}}{w}_m·f\left(r_m\right)$$

(10)

A double sum of this sort being performed for two-dimensional integration.

Then, to account for the given exchange element, J_ab, the entities to be integrated over the 2D grid are

$$X_{mn}^{ab}=r_m^2{r}_n^2{R}_a\left(r_m\right)R_b\left(r_m\right)R_a\left(r_n\right)R_b\left(r_n\right) {\displaystyle \sum }_{m=1}^{m_{max}}{\chi }_{ab}^k\cdot \frac{min{(r_m,r_n)}^k}{max{(r_m,r_n)}^{k + 1}}$$

(11)

The ${\chi }_{ab}^k$ coefficients formalize the definition of J_ab integrals for a given shell couple, ab, as a combination of Slater–Condon parameters with k superscript. The ${\chi }_{ab}^k$ values can be picked from Equation (6a,b) and from Table 1. Conventionally, the m and n indices run on the electrons labeled 1 and 2, respectively, in Equations (1) and (2). One may produce a partial summation, emulating the integration over the electron #2, which corresponds to the one-electron operator in the mean-field treatment of an exchange element:

$$Y_m^{ab}\approx {\displaystyle \sum }_{n=1}^{n_{max}}{X}_{mn}^{ab}{w}_n \to \widehat{ J_{ab}}\left(r_1\right)$$

(12)

The numeric estimation of the whole exchange coupling integrals can be termed as

$$J_{ab}\approx {\displaystyle \sum }_{m=1}^{m_{max}}{w}_m\cdot Y_m^{ab} ={\displaystyle \sum }_{n=1}^{n_{max}}{\displaystyle \sum }_{m=1}^{m_{max}}{w}_m\cdot w_n\cdot X_{mn}^{ab}$$

(13)

We probed that, in the calculations described below, the numerically estimated integrals resemble the fourth digit, or better (in atomic units), the analytic ones, working over a 300 point grid, established with the following parameters: δr₀ = 0.001 Bohr, r_max = 20 Bohr, and h = 0.02.

2.3. The Radial Distribution of the Exchange Energy in Atoms

Now we can attempt the aimed for insight into the atomic structure. For the first check, we will confine the calculations to the B-Ne series, first submitted to CAS or RHF calculations, and afterward reloaded in our codes of analytic atomic calculations, followed by numeric decomposition of the exchange integrals. The RHF refers to the closed-shell case of neon, and to the single-configuration (non-degenerate ground state) of the half-filled p shell in nitrogen. All the other situations imply a degenerate orbital term (P), treated as CAS average over the corresponding three states.

Let us start with the simple task of drawing the radial density profiles, as shown in Figure 1. One may note that going from B to Ne, the area under the curves increased, directly corresponding to the total number of electrons in the atoms. The main changes occurred in the part corresponding to the progressive population of the p shell.

In the following, we produce the density distribution of the exchange. With the leverage described in Equation (12), one may select from the sum of total energy; in Equation (4), the exchange components integrated only on one electron:

$$V_x^{l{l}^{\prime }}\left(r_1\right)={\displaystyle \sum }_{l\le l\prime }{\displaystyle \sum }_{l_i\ne l{\prime }_{i\prime }}\left\{ - p_{l_i}{p}_{l{\prime }_{i\prime }}\left(\frac{1}{2} + 2{\sigma }_{l_i}{\sigma }_{l{\prime }_{i^{\prime }}}\right)\cdot {\widehat{ J }}_{l_{i}l{\prime }_{i\prime }}\left(r_1\right)\right\}$$

(14a)

$$V_x^l\left(r_1\right)={\displaystyle \sum }_{l_i}\left[{\alpha }_l\cdot p_{l_i} + {\beta }_{ l}\cdot p_{l_i}^2 + {\gamma }_l\cdot {\sigma }_{l_i}\left({\sigma }_{l_i} + 1\right)\right]{\widehat{ J }}_{l_i{l}_i}\left(r_1\right)$$

(14b)

Then, summing over all the intra- and inter-shell contributions, we obtained the energy profiles of total exchange density, V_x, shown in Figure 2.

Contemplating Figure 2, one may see that it looks like a “lake” reflection of the “hills” portrayed in Figure 1. This prompts the idea that the exchange V_x(r) ≡ V_x[ρ(r)] in an atom is better related to the radial density distribution, r²ρ(r), than with the ρ(r)^4/3 celebrated formula originating from the theory of uniform electron gas.

Figure 3 details the exchange energy in the X_mn components derived from Equation (11) for the case of a neon atom. The maps for the other atoms look qualitatively similar. Figure 3a sums all the X_mn contributions at each (m,n) point of the 2D radial grid. The other panels detail the distinct shell contributions. Figure 3b shows the interaction between the two occupied s shells, which can be ascribed as $X_{mn}^{1s,2s}$. This shows positive zones, due to the negative 1s(r₁)2s(r₂) areas, whose action was amended with the negative sign of the exchange in the total energy, −2J_1s,2s. However, the total balance of the exchange was negative due to the net positiveness of the exchange parameters. Figure 3c shows the s–p inter-shell exchange, each point containing the $X_{mn}^{1s,2p} + X_{mn}^{2s,2p}$ terms. Figure 3d shows the p intra-shell exchange, small in relative value, as compared to the other parts. For comparability, all the 3D maps from Figure 3 (exchange energy vs. r₁ and r₂ grid points) are drawn at the same vertical scale.

In the following, we will extract a function useful for the discussion of the Fermi hole, h(r₁,r₂). Namely, considering the generic formulation of the exchange energy as integration over ρ(r₁)h(r₁,r₂)/r₁₂, we emulate h(r₁,r₂) multiplying the exchange potential by r₁₂ and dividing by ρ(r₁). The transformation V_xr₁₂/ρ(r₁) is applied to all the exchange components. Thus, with this reasoning, we propose as interesting quantities the following transformations of Slater–Condon primitives involved in the exchange integrals:

$$u_{ab}\left(r_1,r_2\right)=\frac{1}{\mathsf{\rho }\left(r_1\right)}{R}_{n_a{l}_a}\left(r_1\right)R_{n_b{l}_b}\left(r_1\right)R_{n_a{l}_a}\left(r_2\right)R_{n_b{l}_b}\left(r_2\right)\frac{min{(r_1,r_2)}^k}{max{(r_1,r_2)}^k}{r}_2^2$$

(15)

The above expression is the result after converting the r₁₂ factor in terms of minimum and maximum from individual electron coordinates, r₁ and r₂, via multipolar expansion, implicit in the definition of Slater–Condon parameters. Performing such a mutation in the body of each exchange term and summing them in the same way as used for the obtaining of the total energy, one draws the map seen in Figure 4. The standard definition of the Fermi hole follows the integration over the elements of density matrices [22,23]. The above procedure can be taken as an alternative way, useful in the further quest for new empirical recipes for the Fermi hole shape. Obeying the demands of spherical symmetry, we propose here the idea of redrawing the Fermi zone in atoms as a spherical crust, i.e., a radial profile, instead of the actual image, as a local hole in the homogenous or non-homogenous electronic density. Analyzing sections along the r₂ coordinate in Figure 4, one may believe that a sharp Gaussian profile of the Fermi radial distribution may be a reasonable approximation. The subject demands further investigation and come-back studies on different classes of atoms and bases.

3. Methods

The primary calculations were done with GAMESS (General Atomic and Molecular Electronic Structure System) code [24], using the 6-31+G* basis set [25]. The developed analyses were done with original codes written in the MATLAB–Octave environment [26,27].

4. Conclusions

This work takes a constructive critical contribution the actual state-of-the-art in quantum description of the electronic structure in atoms. The actual GTO-based technologies removed the conceptually valuable idea of radial-spherical factorization and the explicit account of the Coulomb, exchange, and correlation effects in the atom by the Slater–Condon parameters. To restore transparency, we derived a handy general formula for the atom in arbitrary shell occupations, obeying the total spin resulted from open shells, as good quantum number. This formula reproduces results from CAS and RHF calculations with dedicated ab initio codes, being implemented with the analytical expansion of all the integrals, in original MATLAB–Octave scripts. In addition, the Slater–Condon parameters were evaluated by numeric quadrature, this way enabling the insight, by analyzing the elementary contributions. In this manner, we drew the radial distribution in the exchange energy density for a prototypic series, including the B-Ne atoms. We found a clear qualitative correlation of the exchange with the radial electronic density distribution. Further numeric handling and corresponding mapping suggested the idea of new density functionals dedicated to the specific of spherical atomic symmetry. Namely, instead of considering the Fermi hole as a local void, it is more appropriate to regard it as a depletion on a spherical crust (i.e., as a radial profile), obeying the spherical symmetry intrinsic to the atom. We will develop this paradigm in future investigations.