Atomic Hartree–Fock limit calculations using Lambda functions

Nonrelativistic Hartree–Fock limit calculations are performed for the group 18 atoms using Lambda functions, which are Laguerre-type basis functions (LTFs). Since Lambda functions form a complete orthonormal set for bound states, the total energy approaches the Hartree–Fock limit monotonically upon increasing the number of expansion terms, as determined by the maximum number (N) of the principal quantum number (n) in a set. The convergence behavior of the total energy in relation to the number of expansion terms is investigated. For Rn, N = 116 is required to satisfy the convergence criterion ∣ Δ E / E ∣ < 10 − 15 . Here E is the total energy and ΔE is (E(N−1) − E(N)). Total energies are obtained with 30 significant digits for He, with 14 digits for Ne, Ar, Kr, Xe, and Rn, and with 13 digits for Og. This approach gives −2.86167999561223887877554374002 au as the Hartree–Fock energy of He.

We have developed a computer program to evaluate molecular integrals over Lambda functions. Using this integral program, we have performed nonrelativistic Hartree-Fock (HF) calculations for single atoms.
Calculation of the Hartree-Fock (HF) energy is a starting point for determining the electron correlation energy. Hylleraas [7,8] used LTFs for the HF, and correlation calculations for He. To date, LTFs have been used mainly in studying light atoms. In the present study, the nonrelativistic Hartree-Fock limit is calculated using Lambda functions for atoms, including heavier atoms. Seven atoms make up the group 18 elements: namely, He (atomic number Z=2), Ne (Z=10), Ar (Z=18), Kr (Z=36), Xe (Z=54), Rn (Z=86) and Og (Z=118). These have electronic configurations of a closed-shell, 1 S.
In heavy atoms, relativistic effects become non-negligible. We are also interested in relativistic calculations such as Dirac-Hartree-Fock (DHF) calculations, however it is well known that DHF calculations in a basis set expansion may cause some troubles such as variation collapse [24,25]. Therefore, to verify the effectiveness of Lambda functions, it would be preferable to perform nonrelativistic calculations at the first stage.
We use here all functions specified by principal quantum number (n) and orbital angular quantum number (l), and we write N for the maximum n in a set. Then N determines the number (N dim ) of entire expansion terms, i.e.the dimension of the relevant Fock matrix. We investigate the convergence behavior of the total energy (E) as N increases. We determine the value of N for which the convergence criterion (| | ) The HF limit values have to date been calculated by the numerical Hartree-Fock (NHF) method. In 1972, Froese-Fischer [26,27] performed NHF calculations for He-Ar and obtained the HF energy to 7 significant figures. The HF energy is vital as a target value in developing basis sets for atoms and molecules. Tatewaki et al [28] extended the calculation for He-Ar to 10 significant figures in 1994. Shortly after, Koga et al [29] attained similar accuracy to Xe. In 2002, Koga et al [30] reached as far as Lr (Z=103). The number of significant figures remains at 10.
To our knowledge, there have been no further cyclopedic NHF calculations for heavier atoms than Lr after Koga et al [30]. However, Saito [31] performed HF calculations for He-Og by applying the B-spline expansion method. The resulting energy values run to 12 significant figures. For He-Lr, the values match those of Koga et al [30] up to 10 significant digits. Thus Saito's results are the most accurate at present. For Rf (Z=104) to Og, no data exist comparable to Saito's. In the relativistic domain, however, numerical Dirac-Hartree-Fock values are available for He-Og [25,[32][33][34].
The present paper determines HF energy values for all atoms of group 18 by applying Lambda functions. We use Lambda functions because they form a complete orthonormal set, so that the total energy that results is guaranteed to converge to the true value upon increasing the number of expansion terms. We are unaware of any Hartree-Fock level calculations using Lambda functions for heavy atoms.

Computational details
The Lambda function used here is written as a product of a radial component and an angular component. A spherical coordinate system is used: Here, Y lm is a spherical harmonic, --+ L n l l 1 2 2 is an associated Laguerre polynomial, and z l is called the exponent. A 2e-integral is written as follows:  Here, r > and r < respectively denote the smaller and the larger of r 1 and r 2 . The integral is separated into angular and radial parts. The radial integral F l in equation (5) corresponds to the formula (1.15) in [35] By the transformation of variables, The function F l can be found simply by multiplying ζ with the other part, i.e. [( ) N N N N 2 1 2 3 4 × integrals]. The latter part is independent of ζ. In the present project we aim to build universal basis sets, and an identical set is therefore applied to all elements. Only ζ is optimized for each element, however.
To evaluate integrals, we used the algorithm developed by Freund and Hill [35]. Their formulation is based on integer arithmetic so as not to lose precision in the integrals. In fact, modern CPUs are designed such that operations for double precision floating-point numbers are processed most efficiently. We therefore wrote a computer program using the multiple-precision floating-point arithmetic package MPfun90 [36], instead of integer arithmetic.
Here, one-electron nucleus-electron attraction integrals are also evaluated using the algorithm of Freund and Hill [35]. Kinetic-energy integrals are calculated according to the formula given by Hagstrom and Shull [37].
In the present study, we use all of the basis functions of the maximum l under n. Thus the dimensions of Fock matrices, N dim , are N for s-electron atom (He), 2N−1 for sp-electron atoms (Ne, Ar), 3N−3 for spd-electron atoms (Kr, Xe), and 4N−6 for spdf-electron atoms (Rn, Og). Here, N is the maximum of n in a set.
For SCF, we used two programs. In the early stage of this project we used the ATOMCI program [38,39], and later on the atmscf program [40]. Both programs have integral generation modules for GTFs/STFs but not for LTFs. Thus the SCF modules of these programs were interfaced to our LTF integral program. We used the ATOMCI program for small N. For large N, we used the atmscf program. However, since the atmscf program requires an initial vector which is quite close to the final solution vector, we at first obtained a solution using the ATOMCI program and transferred this solution vector to the atmscf program.
The original atmscf program is of double precision; we enhanced it to quadruple precision. In the quadruple precision version, the threshold for diagonalization of Fock matrices is reduced to 10 −30 , and the SCF convergence criterion (the relative minimum energy difference) is set at 10 −21 . Since the present basis set is orthonormal, we use unit matrices for overlap matrices. Orthonormality means that HF calculations are free from the problems of linear dependence. The total energy is expected to approach the HF limit monotonically as increasing numbers of basis functions are included.
Expectation values of r n (n=2, 1, −1, −2) for each atomic orbital are evaluated analytically by the algorithm based on the recursion relations (see equations (3.3) and (3.4) of [35]). When integrals of double precision are used, accuracy of 10 significant digits is ensured for these expectation values.

Accuracy of 2e-integrals
The algorithm for evaluating 2e-integrals proposed by Freund and Hills [35] requires many bits of mantissa to avoid significant loss of accuracy when floating-point arithmetic is used. The double precision and quadruple precision use 52 bits and 112 bits for the mantissas, respectively. They correspond to about 16 and 34 decimal digits. In MPfun90, the significant digits is controlled by specifying the value of init. For example, init=200 corresponds to about 200 decimal digits. We first investigated the accuracy of 2e-integrals by varying the init. The values of F l are listed for some values of l in table 1, fixing n at 110 and ζ as 1.
When we specify a sufficiently large init value, we can consider the calculated value as the true value. In table 1, the first example gives the same value until 30 significant digits for init=220 and init=240. Thus this value can be considered as the true value. However, for init=180, significant digits are lost severely and the value diverges from the true value. Figures differing from the true value are underlined. It is also seen from table 1 that significant loss of accuracy occurs readily for (ssss)-type integrals in particular. The reason is that the number of product-sum operations and the repetition of recursion are (n−l−1). In the case of init=220, 2e-integrals have at least 31 significant digits for n=110.

Atomic HF energies (i) He
As the He atom has only 2 electrons, its electronic structure has been thoroughly investigated by many researchers. In 1977, Gázquez and Silverstone [41] performed Hartree-Fock calculations for He using spline functions, and obtained −2.8616799956122 au for the HF energy. Recently, King et al [21,22] performed fully correlated (FC) calculations for He using Lambda functions in the frame of the perimetric coordinate system. They obtained −2.861679995612239 au as the HF energy, using 20 expansion terms.
We performed Hartree-Fock calculations for He with Lambda functions, using integrals of quadruple precision in SCF. By optimizing ζ, we obtained HF energies as shown in table 2. Figures up to 16th digit after the decimal point coincide between energies at N=20 and at N=22. We therefore consider that the number of significant digits at N=20 is 17, and the HF energy value at N=20 is therefore −2.8616799956122388 au. This is same as the value given by King et al [21] up to 16 digits.    We next performed HF calculations by increasing N from 1 to 100, but with ζ fixed as 2. When N is increased sufficiently, the resultant value is expected to approach the true value even without optimizing the ζ value because the Lambda functions form a complete set. Table 3 shows the convergence behavior with increasing N. By comparing tables 2 and 3, it is apparent that the lowest energy values in both tables are independent of ζ. This means that each functional space is expanded using enough terms and the optimization over ζ is saturated. Finally, we present −2.86167999561223887877554374002 au (30 significant digits) as the HF-limit energy of He. Other properties at the N-value for the saturated energy are set out in Supplementary material (stacks.iop. org/JPCO/4/085006/mmedia).
The above results show that HF energies can be calculated to any number of significant figures by increasing the number of expansion terms available online at stacks.iop.org/JPCO/4/085006/mmedia.
(ii) Ne-Rn We continue HF calculations for the remaining heavier atoms. Since these atoms have more electrons than He, the computation is more costly. In this case, we use integrals of double precision in SCF and aim to end up with 14 significant digits for the resulting HF energies. We impose two criteria for the total energy and define N d12 and N d14 : Usually | | DE is appropriate for the convergence criterion. However, since we are concerned with the number of significant digits here, we adopt | | DE E instead. By applying criterion (a) or (b), accuracy of 12 or 14 significant digits is ensured for the total energy, respectively. We refer to the rounded energy value at N d14 as -E HF limit , and present this as a reference value for the Hartree-Fock limit energy. Total energies satisfying these criteria are 14 is the electronic configuration of the ground state of the Og atom. This configuration is consistent with the results of the Dirac-Hartree-Fock [42] and relativistic coupled-cluster [43] calculations.
Since Og has many electrons, a large number of expansion terms are needed to ensure accuracy. This increases the computational cost. As a result, we abandon accuracy of 14 significant figures and seek accuracy only to 13 significant figures. It is appropriate to determine the value of N for which a total energy accurate to 13 significant digits is ensured, and stop the computation at this N-value. For this purpose we give two criteria: We searched the best criterion to obtain 13 significant figures with the lowest computational cost and as a result we obtained the criterion (b) above. Figure 1 shows that the total energy decreases monotonically as N increases. The value of E d12 is −46324.35581505695 au at N d12 =115, and the value of ΔE is 4.15×10 −9 au. Further, E d13 is equal to −46324.35581507534 au at N d13 =125, and -E HF limit is equal to −46324.35581508 au to 13 significant digits. This value matches the HF energy obtained by Saito [31] to 12 significant digits. We confirmed that -E HF limit coincides with the value extrapolated to N=200 by fitting 11 data points between (N d13 −10) and N d13 to the following model function f (t): Here, a, k, and c are fitting parameters. The function f (t) is plotted in figure 1. This model function works well, but it has no physical grounds. For the fitting and plotting, we used Wolfram Mathematica [44].

Conclusions
We have performed nonrelativistic Hartree-Fock calculations for the single atoms of group 18 using Lambda functions. Lambda functions are a type of Laguerre-type basis functions (LTFs), and form a complete orthonormal set for bound states in spherical polar coordinates. The constructed basis sets are universal. Based on the algorithm developed by Freund and Hill [35], one-and two-electron integrals are generated, with their accuracy maintained using multiple-precision arithmetic package MPfun90. These integrals are saved in files in double/quadruple precision and transferred to SCF programs.
By increasing the number of expansion terms (N), we investigated the convergence behavior of the total energy to the Hartree-Fock limit. The total energy decreases monotonically as N increases. Finally, HF energies are obtained to 30 significant figures for He, and to 14 significant figures for Ne-Rn. For Og, the HF energy to 13 significant figures was obtained at N=125. The completeness and orthonormality properties of Lambda functions facilitate precise HF calculation.
The integral program developed in the present project is universal. It is applicable not only to atoms having electron configurations of a closed-shell but also to atoms having an open-shell configuration. In the present article, we showed applications to Hartree-Fock calculations, but it can be used for post-HF calculations such as CI or perturbation calculations. We have also performed preliminary HF calculations for several atoms (Zn, No) rather than the group 18 atoms and obtained promising results.