Numerical Procedures for Relativistic Atomic Structure Calculations

: Variational methods are used extensively in the calculation of transition rates for numerous lines in a spectrum. In the GRASP code, solutions of the multiconﬁguration Dirac–Hartree–Fock (MCDHF) equations that optimize the orbitals are represented by numerical values on a grid using ﬁnite differences for integration and differentiation. The numerical accuracy and efﬁciency of existing procedures are evaluated and some modiﬁcations proposed with heavy elements in mind.


Introduction
Benchmark calculations have shown what can be achieved using the relativistic multiconfiguration Hartree-Fock method [1] as implemented in the General Relativistic Atomic Structure Program (GRASP), the most recent version being GRASP2018 [2]. An example is the paper by Li et al. [3]. Reported were 110, 70, 114, 53 levels of C I-IV, respectively with excitation levels in very good agreement with National Institute of Standards and Technology (NIST) database [4]. Line strengths, weighted oscillator strengths and transition probabilities, were given for radiative electric dipole (E1) transitions involving levels up to 1s 2 2s 2 2p6s for C I, up to 1s 2 2s 2 7 f for C II, up to 1s 2 2s7 f for C III and up to 1s 2 8g for C IV, respectively. The data for each of the transitions were accompanied by an accuracy indicator. The average uncertainty was 8.05%, 7.20%, 1.77% and 0.28% respectively for C I-IV.
In another example, Si et al. [5] report results from extended calculations of spectroscopic accuracy for energy levels and radiative transition rates of Fluorine-like ions for elements with atomic number Z = 24 − 30. Many electric dipole (E1, E2, E3) transition rates are reported along with magnetic dipole (M1) and quadrupole (M2) rates. They believe that the energies and wavelengths are accurate enough to aid line identifications involving the n = 3 and n = 4 configurations, for which observations are largely missing. The calculated wavelengths and transition data will be useful in the modeling and diagnostics of astrophysical and fusion plasmas. Both these examples are for relatively light electrons.
For some heavy, but highly ionized systems, Zhang et al. [6] report benchmark calculations of spectroscopic accuracy for excitation energies and wavelengths in sulfur-like tungsten. Energy levels, wavelengths, and transition parameters involving the 88 lowest levels of W 58+ (W LIX) were calculated. The relative importance of the valence-and core-valence electron correlation effects, the Breit interaction, the higher-order retardation correction beyond the Breit interaction through the transverse photon interaction, and the quantum electrodynamical corrections were shown. The energy levels were highly accurate, with uncertainties close to what can be achieved from spectroscopy.
Benchmark calculations for opacity applications where extensive data are needed, not necessarily of high accuracy, have been reported by Gaigalas et al. [7]. Multiconfiguration Dirac-Hartree-Fock calculations were performed in connection with the analysis of data from the coalescence of a binary neutron star that gave rise to electromagnetic emission powered by radioactive decays of r-process nuclei. These observations provide a unique opportunity to study heavy element synthesis in the universe. However, the atomic data of r-process elements are not yet complete enough to decipher the light curves and spectral features of such kilonovae. In their paper, they performed extensive atomic calculations of neodymium (Nd II-IV, Z = 60 and ground states 4 f 4 6s, 4 f 4 , 4 f 3 , respectively) to study the impact of the accuracy of atomic calculations on astrophysical opacities. They found that wavelength dependent factors of opacity are affected by the accuracy of the atomic calculations.
GRASP calculations played a significant role in a recent paper promoting the complex anion, Th − (Z = 90) , as a candidate for laser cooling. The latter is a well-established technique for the creation of ensembles of ultracold neutral atoms or positive ions. This capability has opened up many new research fields over the past 40 years. However, no negatively charged ions have been directly laser cooled, because a cycling transition is rare in atomic anions. For more than a decade, La − was the most promising candidate, but Tang et al. [8] have evaluated Th − for this purpose. Electron affinities were measured and computed and found to be almost twice as large as previously predicted. In particular, they found that the ground state was 6d 3 7s 2 4 F 3/2 rather than 6d 2 7s 2 7p 4 G 0 5/2 . Consequently, there are several strong electric dipole transitions between 6d 3 7s 2 and 6d 2 7s 2 7p in Th − , supporting anion as a new promising candidate for laser cooling.
Generally, for light elements or highly ionized systems of heavier elements, good accuracy is achieved with GRASP. However, for heavy, more neutral atoms or ions, problems have been encountered that are not yet resolved. For the electron affinity of At (Z = 85), the computed results were not as accurate as had been expected [9,10], at least in part because core correlation had not been accounted for sufficient accuracy. In the lanthanides, the energies of the thirteen levels of the 4 f 2 configuration could not be resolved in Pr +3 [11], where there is some uncertainty about the observed 1 S 0 state and could not be resolved in Ce +2 [12] where all levels have been classified [4]. Difficulties have also occurred in the study of Th II (Z= 90) with a 6d 2 7s ground state [13]. A list of 91 strong Th II spectral lines was established on the basis of a pseudo-relativistic Hartree-Fock model and some radiative decay parameters determined yielding results in reasonably good agreement with the experiment. Attempts were made to validate these results with the fully relativistic code but failed.
These studies suggest that, for heavy elements, several factors come into play that the GRASP code may not have completely resolved.

Seniority or Quasi-Spin Formalism
The underlying theory [1] for GRASP assumes the expansion of the wave function of an atomic state in terms of orthonormal basis states, referred to as configuration state functions. The angular quantum numbers J, M, parity are not sufficient, and additional quantum numbers may be needed. Though the program for generating the basis states still refers to seniority, the latest GRASP angular library relies on the introduction of quasi-spin quantum numbers [14]. The need for quantum numbers beyond seniority arises in the f -shells and, by extension, for g-shells (see the book by Grant [15]). Though no errors have been found in GRASP's angular theory, this is something users need to be aware of. Moreover, extensions will be needed when dealing with the open g-shells.

Core-Core Correlation
Core correlation is a correction that represents the interaction between core-electrons and affects all levels of a spectrum. In fact, there are cases where only the total energy is affected and not the spectrum (determined by energy differences) but these are special cases.
In the CI+MBPT method described by Kozlov et al. [16], core-correlation is expressed as a correction to the potential of the orbitals that define the configuration interaction (CI) matrix for valence-valence and core-valence correlation. The modification is derived from perturbation theory. Thus, the CI wave function has captured core-correlation indirectly-there was no component of the wave function that represented the effect of core-correlation and thus the calculation of transition rates from CI wave functions may need further corrections as well.
The aim of a GRASP calculation is to determine wave functions that can be used to compute other atomic properties. Because correlation is a local effect, the optimum orbital basis for core-correlation is one that is not the same as valence correlation. This was demonstrated clearly in the PCFI method as investigated by Verdebout et al. [17] for Be 1s 2 2s 2 , where biorthogonal transformations were used to compute the CI matrix between non-orthogonal basis matrix elements. The latter is a transformation that, in general, would affect entire wave function expansion that might consist of the millions of CSFs. Further studies are needed. One possibility is to evaluate matrix elements between non-orthogonal basis states by orthogonalizing individual matrix elements, as needed, rather than the entire orbital sets and expansions.

Numerical Accuracy
Calculations for heavy and superheavy elements present numerical challenges. Clearly, the numerical accuracy of orbital properties of the 1s orbital with no nodes will be better than those for the 7s orbital with six nodes and a much larger range of significance. To our knowledge, there has never been a systematic study of accuracy for the procedures used in GRASP in determining orbital energies, for example.

Accuracy of Relativistic Calculations
Relativistic equations differ fundamentally from non-relativistic equations in that the speed of light (or the inverse fine-structure constant, α) is not known exactly. Thus, high-precision calculations will depend on the value that is used. According to the latest CODATA 2018 recommendations, the latter has a numerical value of c = 137.035999084 with a relative standard uncertainty of 1.5 × 10 −10 , or about 10 significant digits. Double precision calculations with 52 bits (about 15 significant digits) are used for representing the value of a variable can involve significant cancellation. In fact, for a 1s electron for hydrogen where n = 1, κ = −1, the formula for the orbital energy e 1s is In this case, = 0.99997 . . ., so that in floating point arithmetic, there are four digits of cancellation in the calculation of (1 − ) yielding an e 1s value with about eleven (11) digits of accuracy for a given c. Thus, numerical calculations should strive for at least 10 significant digits of accuracy for the main contributions to the energy, like the 1s orbitals .
In this paper, we present results of an analysis of numerical procedures for solving the multiconfiguration Dirac-Hartree-Fock equations for atoms [1], as implemented in the GRASP2018 [2] software package and earlier versions, all based on double precision arithmetic. We will refer to these generically as GRASP. In all tables the notation, e-nn, represents 10 −nn with nn often a rough indication of the number of digits of accuracy.
The present study is part of a larger project of rewriting GRASP in an object-oriented fashion adhering to more modern practices for software development. The new code should be easier to understand and modify. It is being developed with heavy and super-heavy elements in mind and will improve numerical accuracy, efficiency and facilitate parallel computing.

Underlying Theory
In the multiconfiguration Dirac-Hartree-Fock (MCDHF) method [1], as implemented in the GRASP2K package [2], the wave function Ψ(γPJ M J ) for a state labeled γPJ M J , where J and M J are the angular quantum numbers and P the parity, is expanded in antisymmetrized configuration state functions (CSFs) The labels {γ j } denote other appropriate information about the CSFs, such as orbital occupancy and coupling scheme. The CSFs are built from products of one-electron orbitals, having the general form where χ ±κ,m (θ, ϕ) are 2-component spin-orbit functions.
In spectrum calculations, where only energy differences relative to the ground state are important, wave functions for a number of targeted states are determined simultaneously in the extended optimal level (EOL) scheme. Given initial estimates of the radial functions, the energies E and expansion coefficients c = (c 1 , . . . , c M ) t for the targeted states are obtained as solutions to the configuration interaction (CI) problem, where H is the CI matrix of dimension M × M with elements Radial functions for orbital a, say, with quantum numbers n a κ a , are solutions of systems of differential equations that define a stationary state of an energy functional for one or more wavefunction expansions. It is possible to derive the MCDHF equations from the usual variational procedure by varying both the large and small component so that where w a is the occupation number of the orbital and V(a; r) = V nuc (r) + Y(a; r) +X(a; r) are potential consisting of nuclear, direct, and exchange contributions arising from both diagonal and off-diagonal Φ α |H DC |Φ β matrix elements [1]. In each κ-space, symmetric Lagrange-related energy parameters ab = n a n b are introduced to impose orthonormality constraints in the variational process. The above derivation, which leads to a system of coupled, non-linear differential equations, assumed the existence of large and small components (P a (r), Q a (r)) in an infinite Hilbert space and applied the variational principle for finding an optimized solution. Alternatively, if we assume the existence of a piecewise polynomial sub-space of order k, with continuous derivatives of order k − 1 over an interval, say (0, R) for which B-splines form a complete basis, the equation for the orbitals become matrix eigenvalue problems. This approach has been described in detail [18] and applied to the case of a single configuration [19,20]. An advantage of the B-spline method is that Lagrange multipliers can be eliminated by projection operators, but preliminary studies [21] indicated that these could lead to numerical instability. Consider a basis of (1s, 2s, 3s, 4s, 5s) orbitals for which there are ten (10) orthogonality constraints and hence four projections for each orbital. When applied to the study of correlation in He 1s 2 1 S using a natural orbital expansion, problems were encountered that have not been resolved.
Differential equation methods have been used successfully in GRASP2018 and, for this reason, are evaluated in this study. New methods will also be more efficient and are expected to be faster than B-spline methods for 10-12 digits of accuracy in the total energy.

Numerical Finite Difference Methods
In finite difference methods, a continuous function is represented by its values on a grid, preferably equally spaced. Formulas for differentiation and integration are derived by approximating the function locally by an interpolating polynomial and integrating (or differentiating) the local approximation.
To test our numerical procedures, we will use hydrogenic wave functions that are known exactly in the case of a point nucleus using the program written by Parpia [22]. In all our tests, we will assume that these functions are exact, that errors arise due to our procedures. Some tests will increase the speed of light so that results are those of a non-relativistic calculation.

The Numerical Grid
The range of the independent variable is (0 ≤ r ≤ ∞), with all orbitals decreasing exponentially at large r. In order to represent orbitals accurately, it is necessary to have a grid where points near the nucleus are closer together than when the points are far from the nucleus. Instead of introducing a variable grid, it is customary to transform the independent variable. In GRASP2018, the independent variable is where Z is the atomic number, c 0 = 2 × 10 −6 and t i the grid with h = 0.05. In this way, the one-electron hydrogenic case has the same t i -grid for all Zr. Then In particular, the array r i = (c 0 /Z)(e ih − 1) represents equally spaced values in the t variable. Thus, for an integral such as the normalization integral, we have A simple but surprisingly accurate integration formula is based on Simpson's rule that integrates over two adjacent intervals, namely When applied to many pairs of adjacent intervals, a simple formula can be derived in terms of sums over odd or even grid-points and a simple end-point correction, namely Here, we have assumed N is odd and that f N = 0.

A Slight Variation to the Grid
A small simplification to a completely exponential grid has significant implications for numerical algorithms, namely Zr = e t or t = log(Zr) so that with c 0 = 2 −19 ≈ 1.9 × 10 −6 and h = 2 −5 ≈ 0.03. Notice that this grid is a universal grid that can be applied to any element of the periodic table. In particular, the grid is the same for any hydrogenic equation. Figure 1 shows the different distributions near the origin for these two grids. This is the critical region. Plotted is every tenth point of the grid as a function of r 1/4 , so that r = 0 is on the graph and the region near the origin is expanded somewhat. The main difference appears to be that the fully exponential grid has smaller step-sizes over the entire region, although the earlier GRASP grid has r = 0 on the grid and the first non-zero value closer to the origin. Although the exponential transformation transforms the (0, ∞) range of r to the (−∞, ∞) range of t, we now have The advantages of this transformation will be reported, but it should also be noted here that by choosing h = 2 −5 to be an exact binary number with only one non-zero bit, the binary arithmetic involved in generating the numerical grid can be executed exactly in many cases. We will also see that arrays of r k will not be needed for computing Slater integrals.

The Fermi Nucleus
Though our numerical tests are for the point nucleus, in practice, most calculations are done for a finite nucleus such as the 'Fermi' nucleus. Figure 2 shows the grid in the region of the origin for a uranium atom Z = 92, A = 235 for both GRASP and present methods. The P(1s; r) orbital for a point nucleus in the same region is also shown. Note that electrons penetrate the nucleus.

Simple Quadrature
An example of straightforward integrations are the normalization or orthogonality integrals, such as (11) Because our transformation started with the small radius r 1 , the range of integration (0, r 1 ) cannot be dealt with in the t variable. Ideally, this contribution should be dealt with from a series approximation of the integrand. However, in many cases, the contribution is negligible, so our initial approach has been to ignore this contribution but monitor the error, which is easily deduced if the answer should be either unity or zero as in orthonormalization. Table 1 shows quadrature results for overlap integrals. Orthonormality results are excellent with about an extra digit of accuracy compared with GRASP2018. Table 1. Errors in present results for overlap integrals for hydrogenic orbitals from calculations with c = c au × 10 5 compared with GRASP and DBSR-HF. The notation, e-nn, represents ×10 −nn with nn often a rough indicator of the number of decimal digits of accuracy.

Integral
Exact DBSR-HF GRASP Present

The One-Electron Integral I (a, b)
A more challenging integral is the calculation of the orbital energy. By multiplying the DHF Equation (5) by the row vector (P a , Q a ) and integrating, it can be shown that the orbital energy is equal to the integral referred to as I(a, a) [15]. The general form consists of four different integrals where, for a point nucleus v nuc (r) = −Z/r. This formula requires the evaluation of the derivative of a function. For a numerical grid of equally spaced points, the derivative could be approximated using Stirling's formula for seven (7) adjacent grid points, three on each side, This symmetric formula is based on seven grid-points, whereas GRASP used a 13-point formula but with a different transformation. Table 2 shows some non-relativistic numerical results for hydrogenic orbitals in calculations where the speed of light was c = c au × 10 5 . In this table, present results are compared with those reported by GRASP and the B-spline Dirac-Hartree-Fock program, DBSR-HF [20]. Generally, the present procedures are more accurate than those of GRASP, but both are less accurate for the first few integrals. In using a larger value for c than the physical speed of light, the value of the function Q → 0. Any error in the calculation of Equation (13) will then be multiplied by the large value of c, resulting in larger than usual errors. With c = 10 10 c au , the relativistic integrals are so small that only the non-relativistic integral is of importance and the error is in Equation (15). Table 3 compares results for a few values from calculations, where c = c au . In this case, GRASP and present results have similar accuracy. In relativistic theory, the one-electron energies for 2s and 2p− are degenerate. Table 3 shows that results degenerate to 13 digits of accuracy. Table 2. The error in I (a, a) results for non-relativistic hydrogenic orbitals compared with GRASP and DBSR-HF. All calculations used c = c au × 10 5 as the value for the speed of light. The notation, e-nn, represents ×10 −nn with nn often a rough indicator of the number of decimal digits of accuracy.

Integral
Exact DBSR-HF GRASP Present

Potential Functions Y k (a, c; r)
In variational methods, the energy is expressed in terms of Slater integrals R k (ab, cd) which contribute to the potential functions in terms of functions Y k (a, c; r), defined as follows: where f (r) = P(a; r)P(c; r) + Q(a; r)Q(c; r). According to this definition, the calculation of Y k functions requires the integration of integrands with factors r k or r k+1 , where for f -electrons k = 0, 2, 4, 6.
For large k, the finite difference approximation may not be sufficiently accurate when approximated by polynomials of lower degree. Because Y k (a, c; r) functions occur frequently in atomic structure calculations, it is desirable to compute them as efficiently as possible. Hartree was developing equations with exchange when Fock published his paper and thus was pre-empted. However, Hartree was also trying to find efficient methods for computing these Y k functions using mechanical calculators before publishing his equations. These methods are computationally the most efficient and stable methods today.
Hartree showed that the Y k and Z k functions are solutions of a pair of first-order equations, namely The first equation can be integrated outward, the second inward. Transforming to the t = log(Zr) variable and multiplying by r, we obtain with boundary conditions Z k (a, c; ∞) = 0 and Y k (a, c; ∞) = Z k (a, c; ∞), respectively. It was shown by Worsley [23] that these equations have simple integrating factors, namely, e kt and e −(k+1)t , respectively. On integrating, we get the equations For our grid, since we have eh = (1 + h), parts of the calculations simplify and are essentially exact at the binary level. At the same time, if we also introduce the variable s such that t = t i + s, then these simplify the equations some more, namely Both require an integration over m intervals, where m is arbitrary. For the first non-zero interval (t 1 , t 2 ), a formula for m = 1 can be used since, in terms of the r variable, the grid-points are very close together and a low order approximation may be sufficient. Thereafter, the ATSP package [24] uses m = 2 along with the Simpson's rule with an error h 5 and can be improved (after the value at t 3 has been computed) with the consequence that the errors on the odd grid points are different from the error at the even grid point. This algorithm for solving for Y k (a, c; r) is stable for k > 0 in that, as the integration proceeds, the error is reduced by the integration factor. The most sensitive case is k = 0. For this value of k, the calculation at large r is an othonormality integral, so that Z k (a, c; ∞) = 1 for i = i and 0 otherwise. Thus, numerical errors are readily detected. The non-relativistic ATSP [24] package corrected the Z 0 (i, i; t) values, since in this case, the error for large values of r is easy to determine from the last few values. In relativistic calculations, with r 1 much closer to the nucleus than in ATSP, this correction was no longer needed.
For integrating these equations over two intervals, a five-point rule was used, namely, Notice that some points outside the intervals of integration improve the polynomial approximation over Simpson's rule.

Slater Integrals
The Slater integrals are defined as where f (r) = P(b; r)P(d; r) + Q(b; r)Q(d; r). Certain symmetries have special notations: R k (ab, ab) ≡ F k (a, b) and R k (ab, ba) ≡ G k (a, b). Symmetry rules apply: for example R k (cb, ad) = R k (ab, cd).
In many applications, such as configuration interaction, numerous Slater integrals are needed and often used many times, especially in large calculations. One effective strategy is to generate lists of all possible Slater integrals for a given orbital set in a canonical form that takes into account all possible symmetries, and evaluate all integrals in advance. This is particularly effective in parallel runs, since the task can easily be distributed among the processors.
In our present object-oriented approach, lists of integrals are ordered by k. Let orbitals be ordered and let i a be the position of orbital a in the list consisting of n orb orbitals. GRASP generates values of Slater integrals effectively "on demand". A more efficient algorithm would be:

For i_a =1, n_orb
For i_c = i_a, n_orb Compute Y^k(ac; r) For i_b = i_a, n_orb For i_d = I_c, n_orb Compute and store R^k(ab,cd) integrals The current practice in GRASP is to compute each R k integral using only integration which requires the generation of r k factors, if not already computed. With orbitals up to l = 6, k may be as large as 12. Such factors are not needed for differential equation methods, making the differential equation method more efficient computationally.

Accuracy of Slater Integrals
To our knowledge, there has been no effort made to evaluate the accuracy of Slater integrals for the relativistic hydrogenic, point nucleus case, where exact orbitals are available.
In non-relativistic theory, Slater integrals are rational numbers. By increasing the speed of light to c = c au × 10 10 , calculations were performed that are essentially non-relativistic. Table 4 shows the accuracy of some such integrals that we refer to as non-relativistic. The F k (a, a) integrals are all for k = 0, although the G k (a, b) integrals also include k = 1. The most accurate results are the B-spline results from DBSR-HF, by 2-3 digits over present finite difference methods, which in turn have about 2 more significant digits over GRASP. It has not been investigated whether the latter is the result of the binary arithmetic related to the grid and the calculations of r k , or the finite difference formulas used. In the present implementation, the tabulation of arrays of r k was not needed.
A piecewise polynomial approximation is based on a grid in which a function is approximated by a polynomial of degree k s in each interval between grid points. In addition, the appoximation and its first k s − 1 derivatives are continuous at interior grid points. Thus, the accuracy of spline approximations depend on two parameters, k s and the grid. In DBSR-HF, the transformation of its grid is not needed and Gaussian quadrature is used to integrate polynomials of degree less that k s exactly. Typically, equally spaced grid parameters are needed near the origin, followed by an exponential grid. By increasing the order of the piecewise polynomial approximation by three orders and halving the grid parameters, we obtained results that for the present purpose are "exact". Table 4. Present errors for non-relativistic Slater integral values for hydrogenic functions compared with GRASP and DBSR-HF. The notation, e-nn, represents 10 −nn with nn often a rough indicator of the number of decimal digits of accuracy.  Table 5 compares the numerical accuracy of present algorithms with those of GRASP and DBSR-HF using splines of order 9 (the default accuracy). All calculations are for the hydrogenic atom with Z = 10. The "exact" values are computed using splines of order 12 and a step-size that is half the default value. We assume that the latter is exact. It is noticed that present F 0 values are considerably more accurate than those of GRASP, such as F 0 (1s, 1s), and occasionally even more accurate than the default DBSR-HF values as for G 2 (1s, 3d−). It appears that the default DBSR-HF grid needs to be refined near the origin, particularly for s-orbitals.
However, Table 5 also draws attention to integrals with a higher value of k. For these cases, the integrand over an interval (see Equation (17)) has a factor r k , increasing the degree of the approximating factor for a given pair of orbitals, requiring a higher-order formula. With B-splines, it is easy to increase the order of the B-spline so that integrands of higher-order can be evaluated exactly. Figure 3 compares the accuracy of present, GRASP, and DBSR-HF values from the standard default order (k s = 9) with higher-order "exact" values. Clearly, the DBSR-HF results are the most accurate. For the point nucleus and the r-grid with default step-size (h s = 0.25), r = 0 is a point of singularity and results can be improved to double-precision accuracy of 15-16 decimal digits by decreasing the step-size near the origin. This is particularly true for 1s (l = 0) orbital and explains the increase in accuracy as l increases. Notice the excellent accuracy as l increases for all k. However, GRASP and present results have larger errors as k increases. GRASP expresses both Z k and Y k in terms of integrals and powers of r k or 1/r k+1 evaluated using five-point formulas and a final quadrature also based on a five-point formula. Present procedures integrate the differential equations over two intervals using a five-point symmetric formula and an exponential grid as explained earlier. The final integration of Equation (27) is computed using Simpson's rule, a three-point formula. Clearly, the error increases with l for both implementations, but the present method is more accurate by about two extra digits. These are deemed sufficiently accurate for current applications, but seven (7) point formulas for integrating the differential equations should reduce the increase with respect to l and k. Alternatively, formulas could be developed for integrating functions that have a weighting factor r k .

Timing
Accuracy often comes at the cost of more computational effort. As a final check, we computed the average time needed for calculating a Slater integral by running the code for numerous integrals. Table 6 shows our timing results on a standard, serial computer. It shows that the B-spline methods, though accurate to almost 16 digits, have times best measured in milliseconds, whereas both GRASP and present methods compute Slater integrals in microseconds (µs).

Solution of the DHF Equations
Since DHF equations are non-linear (the potential, for example, depends on the radial functions being computed), the equations are generally solved iteratively until results are converged or "self-consistent". Thus, generally, estimates of the solution are always available. Including corrections from a previous iteration will be referred to as a deferred correction.
Equation (5) expressed DHF equations as a non-linear integro-differential equation of eigenvalue type. Separating out theX(a; r) contribution to the potential and including this on the right-hand side simply as X(r), we obtain a nonhomogenous equation which, for r = e t , has the form, where, u(r) is a two component vector (P(r), Q(r)) t and In the r variable, S = I, the identity matrix, but in the transformed form the diagonal elements are S ii = r. The function −Z nuc (r) = −Z + rV nuc (r), where the latter term takes into account the effect of the finite nucleus, although it is omitted in the present study of numerical accuracy. In this form, it has been assumed that there are no off-diagonal Lagrange multipliers and the factor c for the derivative has been retained as in Ref [15]. Note also that both P(r) and Q(r) decay exponentially for large r.
In many numerical methods for differential equations, the derivative would be replaced by expressions in terms of finite differences on the grid and equations established that the solution needs to satisfy. In this work, the equation is integrated and integrals are approximated by expressions in terms of values of the unknown on a grid.
Consider the simple equation y = f (r) with the boundary condition y(0) = 0. The solution can be expressed as Needed is an approximation for the integral. In this work, we will use Simpson's rule with a correction, namely r i+1 In the DHF equations, the function f (r) is replaced by f (r)y(r), so that the unknown quantities also appear in the integral. This is accounted for in the terms that comprise Simpson's rule, but the correction becomes a deferred correction, computed from an estimate of the solution. The differential equation is then replaced by a set of linear equations with a three-term formula. Equation (30) applied at the midpoint of the interval, r i , i = 2, . . . , N, and boundary conditions at i = 1 and i = N + 1, leads to N − 1 equations in N − 1 unknowns.
In the solution proposed here, the equation at r i is obtained by integrating from r i−1 , r i+1 on an exponential grid, namely where the function rF(r) = rV(r) + S(r).
In the latter, the orbital energy is assumed to be positive. The matrices evaluated at a grid point i are denoted as V i , S i , and F i respectively and the 2-component vectors are u i . It is also useful to define the matrix where c is a constant representing the speed of light.

Outward Integration
GRASP integrates outward over a single interval (r i , r i+1 ) at a time, using the mid-point rule with higher-order corrections deferred (computed from previous estimates in the SCF process). In this work, integration is done over two intervals so that the highly accurate Simpson's rule can be used for computing the new values of u i+1 in terms of the previous values for (i − 1, i). A higher-order deferred correction can be included when necessary, as done in ATSP [24].
Expressing each of the integrals in terms of values at (i − 1, i, i + 1), then rearranging the equations so that the unknown values at i + 1 are defined in terms of the previous values (i − 1, i), we obtain the three-term relationship, which requires that u 1 and u N + 1 are either zero or known. In essence, the equations are systems of equations where all entries are 2 × 2 matrices and the solution is an array of 2-component vectors. For the DHF equations, the matrices arise from using Simpson's rule are The R-array includes three contributions.
• The integration of theXu(r) function. Since this contribution is computed entirely in terms of functions that are already known, a five-point integration procedure can be applied directly.

•
The deferred correction for the function V(r).

•
The deferred correction for the S(r) function. Note that, since r = e t , the integration of S(r) with respect to t is not exact.
When solutions are known at i − 1, and i, the solution at i + 1 is predicted as The process is stable in that F(r) < 0 near the origin. Outward integration is terminated when F 11 (r) → 0.

Near the Origin
The boundary conditions for the radial equations require that P(r) and Q(r) be zero at the origin and infinity. It is customary to analyze the behavior of the solutions near the origin in terms of power series expansions [15] in r. The results depend significantly on whether the nucleus is a point nucleus or finite nucleus and also on κ.
In the present work, the first grid point is not at the origin. Instead of applying a different analysis depending on whether or not there is a singularity at r = 0, we could make the assumption that the differential equation at r i , i = 1, 2 needs to be satisfied. Then, Since DHF equations are solved iteratively, estimates are always available. Often, these are from analytic expressions for hydrogenic orbital, using the program Dirac-Coulomb wave function program DCWF [22].

Inward Integration
Outward integration may become numerically unstable in the outer tail region. At the same time, for the homogeneous equation, solutions satisfying both boundary conditions only exist for discrete eigenvalues, so part of the numerical problem is to also find these eigenvalues.
Inward integration is similar to outward integration, the only difference being that unknown values are at the point r i−1 rather than r i+1 . However, the problem of deciding where inward integration should start can be avoided by treating the equations in the tail region as a system of linear equations. Then, the 2 × 2 matrix form of LU factorization can be used, that assures that the two regions join at a point, and asymptotic boundary conditions can be applied to determine whether solutions from L-factorization are satisfactorily small. When such a point has been reached, backward substitution can be applied.
The situation is somewhat different in the case of the three-term relation, in that the equation in the outward region at i = j predicts u j+1 , whereas the first equation in the tail region i = j + 1, uses the same value for u j , but also computes values u j+1 . Thus, there will be differences in both components at i = j + 1 that can be used to determine changes to the energy. Alternatively, given estimates of the radial functions, a new energy can be computed so that radial functions and orbital energies are determined iteratively.

DHF Equations as a System of Linear Equations
Let us introduce some 2 × 2 matrices, and define equation i as the equation for integrating over the interval (r i−1 , r i+1 ). Then, where and C i includes other corrections, including exchange and the deferred corrections. As r → ∞, these matrices become constant (see Equation (29)). At the join, u j is known. Values at j + i, i = 1, 2, 3, . . . M, where M is such with |P M | ≤ 10 −10 for j+ i > M need to be computed. Then, the first equation for the system of equations for the unknown variables is (since u 1 is known) and those for i > 2 Applying a matrix form of LU decomposition, we eliminate u i−1 from the next equation by computing A i−1 D −1 i−1 and subtracting a multiple of the equation for i − 1 from i'th equation to obtain Given the 2 × 2 matrix, the inverse can readily be computed. When where d = a 11 a 22 − a 12 a 21 . The LU decomposition transforms the tridiagonal system of equations to the upper diagonal form of To be determined is the last value of i. In the region for large r the matrices are approaching a constant value. The last equation for which the i − 1 variable has been omitted is of the form In the exponential region, both large and small components decay exponentially so that P(r i+1 = e − (r i+1 −r i ) P i , for example. Using this information, Equation i can be used to estimate the solution u i . If |u i | is sufficiently small, the equation defines M, the largest non-zero value and backward substitution determines the solution in the tail region. For a homogeneous equation, the join may not be perfect because was not the desired eigenvalue.

The Lagrange Multipliers
So far, we have not considered the energy parameters that enter into the DHF equations because of orthogonality constraints to other orbitals with the same value of κ. Let the orbital under consideration be u a and let u b another orbital of the same κ symmetry. Then, the DHF equations have the form Notice that, this form differs from Equation (5), in that the equations have been divided by the generalized occcupation number and include factors of r needed for the t variable. As a result of the divide, the matrix of energy parameters ab is no longer symmetric when w a = w b , whereas the original matrix of Lagrange multiplier is symmetric. Multiplying the equation by u b (t) and integrating, defines ab . Notice that the differentiation of orbital u a is required. In fact, the I(a, a) discussed earlier is such an integral for the one-electron orbital energy. The variational method requires that the total energy be stationary for all variations that satisfy the boundary conditions. The off-diagonal energy parameters are associated with orbital transformations that are orbital "rotations" [1]. When two orbital subshells are both filled, such rotations change neither the wave function nor the energy, but the "best" solution is one where ab = 0, in that each orbital has one exponential tail. There are others, like 1s2s 3 S, where the off-diagonal energy could be set to zero, but these are special cases that can only be detected by a rotational analysis. Another case is the case of a complete expansion where the wave function is invariant under rotation. The best solutions in those cases are in terms of natural orbitals, namely orbitals from CSF expansions, such as 1s 2 + 2s 2 + 3s 2 . . ., etc. In effect, the degrees of freedom for rotations are used to eliminate CSFs rather than setting Lagrange multipliers to zero.

Tests for Accuracy of the Hydrogenic Equation
The differential equation methods for the hydrogenic equations can readily be evaluated to the point nucleus. As for previous tests, accuracy is assessed by observing how well the method predicts solutions when provided with estimates of machine accuracy and comparing the computed results with exact results. Table 7 reports errors for some one-electron cases of Be 3+ (Z = 4) for several different values of κ. For a point nucleus, the behavior of the large radial component ( [15]) is The process is one where the estimated solution is essentially exact, and so deferred corrections can immediately be determined. The first step is to compute the orbital energy from the function which is similar to computing the I(a,a) integral that we already discussed in an earlier section. Using the computed energy and the two initial estimates, inward and outward integration is performed. The maximum difference in this solution from the initial exact function determines the error in the computed function. The deviation from unity of the normalization integral is another test. In an iterative process, the square root of this integral is used to normalize the radial functions, and the process is repeated. For the present method, only one iteration was applied. Table 7 shows that, for the s-orbitals, the 1s energy is computed the most accurately. There is a slight decrease in accuracy as the energy decreases and the range expands with more grid-points included in outward integration. As κ increases, there is a slight improvement in the energies for higher values, possibly because the omitted contribution from the (0, r 1 ) is less important as κ increases. The accuracy of other parameters can be expected to decrease as the number of points for outward integration increases.
Most important is how these parameters compare with the current GRASP code. Table 8 compares the present methods with those of GRASP2018, based on a mid-point rule for integration along with deferred corrections. Unlike the present method where the orbital energy was obtained in the same way as the off-diagonal energy parameters, the GRASP code uses energy adjustment procedures, outward-inward integration, variational equations that define how the radial functions change as the orbital energy changes, and series expansions at the origin. In Table 8, three quantities are compared-the orbital energy associated with the new radial functions, the difference between the new radial function for the large component, and the normalization integral of the new functions. For some reason, the least accurate orbital energy is that for 1s, which is somewhat surprising. The normalization integral for 1s, 2p, 3d orbitals is an order of magnitude too small resulting in the large error in d NORM for the first iteration. This suggests that there might be errors in obtaining the correct starting values but all calculations terminated after three (3) iterations with only small improvements. The errors in the large component after normalization have better accuracy but are not as accurate as the present values.
After three iterations, dP max (1s) < 1 × 10 −7 whereas dP max for 2p, 3d did not improve significantly with additional iterations. Table 7. Results showing the accuracy of the one-electron equation for Z = 4 from the results based on Simpson's rule. Shown are the absolute errors in the orbital energy computed from the "exact" wave functions (dE), the maximum change in the computed radial function (dP max ) before normalization, the error in the normalization integral (d NORM ) for the solution of the differential equations, the location of the join (n join ) of outward and inward integration, the range of the solution (n max ) and difference in the large component at the join (dP join ). The notation, e-nn, represents ×10 −nn with nn often a rough indicator of the number of decimal digits of accuracy.

Conclusions
New methods for performing relativistic calculations based on Dirac theory have been proposed. The three-point Simpson's rule is shown to provide excellent results when applied to equations based on the universal exponential grid Zr = e t or t = log(Zr), along with deferred corrections that are readily applied in an iterative process. For bound state solutions, the range of integration is partitioned into the outward integration range and a tail region in which a solution is obtained that maintains continuity at the join and determines the last point. The method could also be applied to obtaining continuum solutions that require only outward integration. The range of integration would now be partitioned into the exponential grid region joined with a uniform grid at large r. A solution in the latter would be obtained that maintains continuity and applies a boundary condition at the end of the uniform region. The procedure has not yet been implemented. Further studies are also needed to solve the MCDHF equations for multi-electron systems.