Highly accurate numerical solution of Hartree–Fock equation with pseudospectral method for closed-shell atoms

The Hartree–Fock (HF) equation for atoms with closed (sub)shells is transformed with the pseudospectral (PS) method into a discrete eigenvalue equation for scaled orbitals on a finite radial grid. The Fock exchange operator and the Hartree potential are obtained from the respective Poisson equations also discretized using the PS representation. The numerical solution of the discrete HF equation for closed-(sub)shell atoms from He to No is robust, fast and gives extremely accurate results, with the accuracy superior to that of the previous HF calculations. A very moderate number of 33 to 71 radial grid points is sufficient to obtain total energies with 14 significant digits and occupied orbital energies with 12 to 14 digits in numerical calculations using the double precision (64-bit) of the floating-point format.The electron density at the nucleus is then determined with 13 significant digits and the Kato condition for the density and s orbitals is satisfied with the accuracy of 11 to 13 digits. The node structure of the exact HF orbitals is obtained and their asymptotic dependence, including the common exponential decay, is reproduced very accurately. The accuracy of the investigated quantities is further improved by performing the PS calculations in the quadruple precision (128-bit) floating-point arithmetic which provides the total energies with 25 significant digits while using only 80 to 130 grid points.


Introduction
The Hartree-Fock (HF) approximation [1][2][3][4][5] provides a simple description of a many-electron system in terms of one-electron spin-orbitals. It explicitly accounts for the fermionic character of electrons by using the determinantal form of the manyelectron wavefunction, though it neglects the correlations due to their mutual Coulomb interactions. The Slater determinants (or their combinations, for open-shell systems [6]) obtained for ground and excited states in the HF approximation are usually a good starting point in more sophisticated quantum-mechanical calculations which include electron correlations, like the Møller-Plesset perturbation theory [7], the coupled-cluster approach [8] or the quantum Monte-Carlo (QMC) method [9].
The HF one-electron spin-orbitals and their energies satisfy the HF equation including the Hartree (direct Coulomb) potential and the non-local (integral) exchange operator which both depend on the orbitals themselves. The equation takes the simplest form in the restricted HF method for the closed-shell systems where the orbitals are doubly occupied and the valence shell is fully filled. The accurate numerical solution of this eigenvalue equation for atoms with closed (sub)shells is the main focus of the present paper. The extension of the restricted HF method to open-shell systems [6,10,11] results in the Hamiltonians (Fock operators) which are different for each atomic subshell. As a result, the respective HF equations include off-diagonal Lagrange multipliers so they are not eigenvalue equations. Alternatively, these equations can be transformed into an eigenvalue problem by introducing a unified coupling operator, common for all the atomic orbitals with the same orbital number. However, this operator is not unique and its dependence on orbitals is more complex than for the Fock operator.
The HF integro-differential equation is usually solved by representing molecular orbitals in an atomic orbital basis (of Gaussian or Slater type) and subsequent diagonalization the resultant Roothan algebraic eigenvalue equations [6]. However, in the case of atoms and small molecules, alternative numerical methods based on finite-difference [12][13][14][15][16][17][18], spline [19][20][21][22][23][24], finite-element [25] or wavelet [26] representations are also feasible. The so obtained orbitals and their energies can be highly accurate and can serve for comparison with less accurate results obtained in calculations with orbital bases. The numerically obtained orbitals have also been used when a high accuracy is required, like in studies of the node structure of the HF orbitals [27]. This structure has an important role in the QMC calculations due to the fermion sign problem [28] since it is often assumed that the nodal surfaces of the correlated wavefunction are fixed at those of the trial HF wavefunction [9].
A very accurate solution is also possible in the case of the many-electron Schrödinger equation, with a numerical effort much smaller than in the full configuration-interaction variational approach. The iterative complement interaction (ICI) method proposed by Nakatsuji [29][30][31] takes advantage of the analytical form of the many-electron Hamiltonian and finds an analytical expansion of the many-electron wavefunction, with potentially arbitrarily high accuracy, by an iteration using of a suitably scaled Schrödinger equation. The total energies with 7 to 16 significant digits have thus been obtained for small atoms and molecules Journal of Mathematical Chemistry (2020) 58:1571-1600 [31][32][33] and the remarkable 43-digit accuracy has been reached for the total energy of the He atom [34]. However, the ICI method is not suitable for determining the HF many-electron wavefunction (though this function can be used as a initial guess [34]) since the applied iteration would not preserve its determinantal form required in the HF approximation. Also, an application of this method to solve the herein investigated HF equation for one-electron orbitals is hindered by the presence of the non-analytical (Hartree and exchange) terms in the Fock operator.
Numerical solution of differential equations can also be obtained with the pseudospectral (PS) method which is based on a combination of the spectral and real-space representations. It has been found to be very accurate in solving onedimensional Schrödinger equation [35][36][37][38][39][40], in particular the Kohn-Sham (KS) equation for atoms [41,42] described with the density functional theory (DFT). The KS equation includes only the local (multiplicative) effective potential which makes the application of the PS method straightforward once suitable scaling is applied. The presence of the non-local orbital-dependent exchange operator in the HF equation makes the application of the PS method more complex. The approach proposed by Friesner [43,44] over three decades ago combines expressing molecular orbitals in an atomic basis with the use of the PS representation for calculating the Coulomb integrals with a real-space quadrature. Although this approach was successfully applied to polyatomic molecules [45,46], it does not provide the PS solution of the HF equation using the real-space representation of the orbitals and its accuracy depends on the applied basis. The purely PS representation of the HF equation has previously been formulated for atoms in intense magnetic fields using two-dimensional real-space grid in cylindrical coordinates [47] but its application was limited to helium and lithium atoms [48].
The present paper derives the one-dimensional (radial) PS representation of the HF equation for atoms at zero magnetic field and applies it to closed-(sub) shell atoms with the atomic number up to Z = 102 . The developed method is different from Friesner's approach since the HF orbitals are now represented directly by their values on a spatial grid which are found by solving the discretized (algebraic) HF equation. The use of the PS representation provides very high accuracy of the resulting orbitals for a moderate number of grid points. This accuracy is further enhanced when the precision of the floating-point format used in the calculations is extended from double to quadruple. The obtained results (orbitals and their energies, total energy and its components, the virial ratio, and the quotient defining the Kato cusp condition at the nucleus) are compared with the results of previous numerical HF calculations, in particular by Bunge et al. [49] and Tatewaki et al. [50] who have obtained the total energies of atoms with ten-digit accuracy using the HF code by Froese Fischer [13,14]. The comparison is also done with the results of the Roothan-Hartree-Fock (RHF) calculations based on the expansions of the orbitals in the Slater-type basis [51,52] and in the B-spline basis [21]. The latter basis is formed of specific splines which are piecewise functions composed of polynomials of a small order M and each of these basis functions vanishes outside a radial interval including exactly M neighbouring knots while different functions correspond to different intervals.

3
The present method for the PS solution of the HF-like equations was already used for investigating the DFT exchange potential in the previous works by the present author [53,54]. In particular, the very high accuracy of this method was crucial for determining the correct dependence of the orbital-specific HF exchange potentials at large distances r from the nucleus in atoms. Surprisingly, this dependence is given by c 1 + c 2 ∕r with constants c 1 < 0 and c 2 > 0 for most orbitals, with the sole exceptions of the highest occupied orbital (HOMO) and atoms with s occupied orbitals only [53,55]. However, other results for the HF atomic ground states and the actual PS formalism for the HF equation were not reported therein.

Hartree-Fock equation
In the unrestricted HF theory the many-electron wave function describing a system of N el interacting electrons moving in an external potential v ext ( ) is approximated by a single Slater determinant built of one-electron spin-orbitals a ( ) ( a = 1, … , N el , N ↑ el + N ↓ el = N el ) dependent on the electron spin ( ↑ or ↓ ) and the electron position . In the ground state, these occupied (occ) orbitals a ( ) and their energies a satisfy the HF equation obtained by minimizing the mean value ⟨ �Ĥ� ⟩ of the system Hamiltonian Ĥ over the subspace of normalized determinants . This equation includes the Hartree potential v H ( ) and the external potential v ext ( ) which are due to the electrostatic interaction between the electrons and their interaction with the external field of the atomic nuclei, respectively. The non-local Fock exchange operator v F x ( ) is built of the occupied HF orbitals a ( a = 1, … , N el ) with spin and its action on an oneelectron wave function ( ) yields [5] The Hartree potential v H ( ) = ∫ d � n tot ( � )∕| � − r| is found with the total electron density n tot ( ) = n ↑ ( ) + n ↓ ( ) , given by the sum of the spin-projected densities n ( ) = ∑N el a=1 � a ( )� 2 . The HF equation needs to be solved selfconsistently. For sake of simplicity, the present paper considers closed-(sub)shell systems for which the orbitals a ( ) = a ( ) are the same for both spin directions and they are doubly occupied: N ↑ el = N ↓ el = N el ∕2 . This corresponds to the restricted HF method where the HF equation is common for the two spins so that the spin label can be skipped hereafter.
For a closed-(sub)shell atom in its ground state, the Fock exchange operator, acting on an atomic orbital a ( ) = r −1 nl (r)Y lm ( , ) ( a ≡ nlm ), yields

3
Journal of Mathematical Chemistry (2020) 58:1571-1600 where Y lm ( , ) is the spherical harmonic and n, l, m are the principal, orbital, and magnetic quantum numbers, respectively. The factor (where the step in the summation over l ′′ is 2 which is marked with the prime next to the sum symbol) is defined [5] with the functions where r < = min(r, r � ) , r > = max(r, r � ) , and the Gaunt coefficients which are expressed in terms of the 3j Wigner symbols [5,56]. The Hartree potential depends on the total radial density The resulting HF equation for the radial orbitals nl (r) reads Note that both the angular and radial parts of the atomic orbitals a ( ) are normalized to unity so that we have ∫ ∞ 0 dr 2 nl (r) = 1 and, accordingly, the radial density (8) present in Eq. (7) integrates to the total number of electrons, In numerical calculations on a finite grid the functions v l �� (n � l � , nl;r) defining the Fock exchange operator and the Hartree potential v H (r) can be conveniently found by solving the corresponding Poisson equations These equations, which originally included the radial part r −1 [d 2 ∕dr 2 ]r of the Laplacian, have been multiplied by r on both sides, and, simultaneously, the modified quantities and ṽ H (r) = rv H (r) have been introduced. In a similar way, the Eq. (9) for the radial orbitals is obtained from the HF Eq. (1) after multiplying its both sides by r.
In the present paper, the Poisson and HF equations are solved numerically with the PS method in a similar way as it is done for the KS equation in the previous works [35,41,42].

Pseudospectral method-general considerations
The PS method [57,58] is based on the (approximate) representation of a sought quantity f(x), depending on a variable x, as a linear combination of cardinal functions g j (x) where the expansion coefficients are values of f(x) at the grid points x j ( j = 0, … , N ). For boundary value problems the Lobatto-type grid is usually adopted so that the points x 1 , … , x N−1 are defined as the zeros (nodes) of the derivative P � N (x) of a polynomial of a high order N. In present work, the Legendre polynomials P N (x) with N up to 200 are used and the zeros of P � N (x) are found numerically with the subroutine ZELEGL by D. Funaro [59]; the two end (boundary) grid points used in Eq. (13) The cardinal functions are the polynomials of order N that satisfy the relations (i, j = 0, … , N ). This makes the representation (13) accurate at all the grid points x i while it is approximate for a general point x ∈ [−1, 1] . The cardinal functions are in fact Lagrange basis polynomials but they can be simply expressed in terms of the polynomial P N (x) as [35,41,57] The PS solution of the HF problem also requires the values of the second derivative of the cardinal functions at the grid points [35,60], r .
(12) v l �� (r) =ṽ l �� (n � l � , nl;r) = rv l �� (n � l � , nl;r) . Their first derivative is needed in the calculation of the kinetic energy (see Sect. 2.5) and it is expressed as follows [57] If the value of calculated function f(x) needs to be determined also at x points between the grid points, the original definition of g j (x) , Eq. (15), can be replaced with an equivalent expression which is the expansion of the cardinal function in terms of the Legendre polynomials and does not include the denominator x − x j , troublesome at x close to x j . This expression is obtained from the summation formula [61] (valid for any x and a, in particular for a = x j ) combined with the identity [57] applied for any x, and, separately, at the grid points x = x j where the function (1 − x 2 )P � N (x) vanishes for all j = 0, … , N . To the author's best knowledge, the formula (20) has not been reported before. Its computational cost is similar to that of the original expression (15) for g j (x) since the recurrence relations that are used to find P � N (x) simultaneously yield P k (x) , k = 1, … , N. The grid points x j can also be used as the nodes of the highly accurate Gauss-Legendre-Lobatto (GLL) quadrature where the weights are defined as [57] The GLL quadrature is used for numerical integration in this work. For this purpose the integral of a function q(r) over r is first transformed to ∫ 1 −1 f (x)dx , where f (x) =̇r(x)q(r(x)) and ̇r = dr∕dx , by applying the following mapping where L is a constant. Accordingly, the grid points x j are mapped to the points r j = r(x j ) ( j = 0, 1, … , N ) of the radial grid with r 0 = 0 and r N = ∞ . There are formally N + 1 points r j , but, since the outermost point r N is at the infinity, the actual number of the radial grid points is N and it is so denoted in this paper.

Fock exchange operator and the Hartree potential in the pseudospectral method
The present paper determines the PS representation of the Fock exchange operator and Hartree potential. The mapping (25), previously used in the PS solution of the KS equation, [35,41,42,62] is now also applied in solving the HF and Poisson Eqs. (9), (10) and (11) with the PS method. It is accompanied by suitable scaling of the sought functions (multiplying them by (̇r) −1∕2 ) to avoid the occurrence of their first derivative over x in the scaled HF and Poisson equations which would lead to non-symmetric matrices in the final discretized forms of these equations. Subsequently, the PS representation of the scaled functions is used in the transformed equations. In particular, for the functions ṽ l �� (n � l � , nl;r) that define the exchange operator through Eq. (4) the scaling reads and the scaled functions f (x) = p l �� (n � l � , nl;x) are further represented as in Eq. (13). As a result, the Poisson Eq. (10) is transformed into the algebraic equation The matrix A = (A ij ) is defined as and it is symmetric according to Eq. (16). Equation (27) is solved for = (y 1 , … , y N−1 ) and its solution = A −1 is found for each of the several right-hand sides = (b 1 , … , b N−1 ) corresponding to different quantum numbers {nl, n ′ l ′ }. This multiple operation does not increase significantly the computational cost of finding F x;nl (r) at the points r = r i ( i = 1, … , N − 1 ) since the matrix A is independent of nl and n ′ l ′ . However, the solutions still need to be obtained separately for different values of In practice, prior to actually solving Eq. (27) for , the matrix A can be inverted for l �� = 0, … , 2l max where l max is the highest orbital number l of the occupied HF orbitals �nl⟩.
Equation (27) accounts for the boundary conditions at the end points x 0 = −1 and x N = 1 . As a result, the j = 0 and j = N terms are removed from the sum of A ij y j for 1 ≤ i ≤ N − 1 since A i0 y 0 = 0 and A iN y N = 0 . The latter relations are obtained by noting that v l �� (n � l � , nl;r) ∼ O(r l+l � +2 ) + O(r l �� ) for r → 0 and v l �� (n � l � , nl;r) ∼ 1∕r l �� +1 for r → ∞ and accounting for the variation of ̇r(x) near the two end points: ̇r = L∕2 at x = −1 (i.e., at r = 0 ) and ̇r = (2L) −1 r 2 for x → 1 (i.e., for r → ∞ ). These boundary conditions also lead to vanishing of ∑ N j=0 A ij y j and b i for i = 0 and i = N so that the equations for the two values of i are not included in Eq. (27). Thus, the (N − 1) × (N − 1) linear problem, given by Eq. (27), finally emerges.
The PS representation of the Poisson Eq. (11) for the Hartree potential v H (r) also leads to Eq. (27), with The value of the Hartree potential at the nucleus ( r = 0 ) can be found directly from its definition using the GLL quadrature: where the i = 0 and i = N terms are excluded since the integrand vanishes at r = 0 and r = ∞.

Pseudospectral representation of the Hartree-Fock equation
The radial atomic orbitals nl (r) are scaled in a similar manner [35,41,42,62] (28) as in Eq. (26) and the function f (x) = nl (x) is represented for each (nl) as in Eq. (13). This PS representation is used in Eq. (1) and combined with the PS solutions for v H (r) and F x;nl (r) obtained from Eq. (27) at the grid points r = r i . It results in the following algebraic equation which is the sought PS representation of the restricted HF Eq. (9) for closed-shell atoms. The dimensionless quantity (i = 1, … , N − 1 ) represents the atomic orbitals at the grid points. The exchange matrix and the Hartree potential present in Eq. (33) are expressed in terms of u (n � l � ) i for the occupied orbitals and the matrices A(l �� ) as where and Note that both nl (r i ) and u (nl) i vanish at the end points r 0 = 0 ( x = −1 ) and r N = ∞ ( x = 1 ) so that the terms with j = 0 and j = N are absent in the PS representation of the HF Eq. (33).
Since the matrices (g �� ij ) and A(l �� ) are symmetric the matrix expressing the HF hamiltonian (the total Fock operator) in the PS representation (33) is also symmetric. This allows for its diagonalization with standard numerical techniques, which gives the energies nl and the radial wavefunctions nl (r i ) = P N (x i )(̇r i ) −1∕2 u (nl) i of the atomic HF orbitals. The exchange matrix F (l) x;ij and the Hartree potential v H (r i ) enter- ing Eq. (33) depend on the orbitals, thus the iteration of the HF solution is needed to achieve its convergence.

Hartree-Fock energy and its components
The energy of the external (electron-nucleus) interaction v ext (r) = −Z∕r , the electrostatic (direct Coulomb) energy and the exchange energy are found from their integral definitions by the application of the GLL quadrature: where the i = 0 and i = N terms are excluded since the respective integrands vanish at r = 0 and r = ∞.
The kinetic energy, defined as the sum of the integrals of −(1∕2) * a ( )∇ 2 a ( ) for the occupied orbitals a = nlm of both spins, can also be found with the GLL quadrature where i 0 = 0 for l = 0 and i 0 = 1 for l > 0 . The derivatives � nl = d � nl ∕dr are found by expressing them as � nl = (̇r) −1∕2 � nl + (2L) −1∕2 nl with Eq. (32) and applying the PS representation (13) (18). Alternatively, once the energy components E ext , E es and E x are found, the kinetic energy E kin can be immediately calculated as This known formula is obtained by multiplying the HF Eq. (1) by * nlm ( ) , followed by the integration of the resultant equation over and the summation over all occupied orbitals of both spins.
The total energy is finally found as the sum of all its components The accuracy of E tot can be estimated by checking how well the virial theorem E pot = −2E kin , or equivalently E tot = −E kin , is satisfied. This theorem involves the total potential energy E pot = E ext + E es + E x .

Kato cusp condition
For atoms the Kato cusp condition is defined for s orbitals by ratio of the derivative of the wavefunction n00 (r) = (4 ) −1∕2 r −1 n0 (r) to itself at the position of the nucleus: This quotient should be equal to −Z according to the Kato theorem This theorem is also valid for the ratio q tot = n � tot (0)∕n tot (0) which, again, should be equal to −Z. The values of the total electron density n tot (r) and its derivative at the nucleus ( r = 0 ) are calculated as follows where � n0 (0) (here replacing n00 (0) ) is found with Eq. (43) for l = 0 and i = 0 . The second derivative of n0 (r) at r = 0 is calculated with the following formula at i = 0 and l = 0 , with g ′′ ij defined in Eq. (16). This expression is obtained from the scaling relation �� nl (r) = (̇r) −3∕2 �� nl (x) and the PS representation (13) of f (x) = nl (x) . Let us note that the formula (49) and a similar expression for d 2ṽ l �� (n � l � , nl;r)∕dr 2 at r = r i in terms of g ′′ ij and y j ( j = 1, … , N − 1 ) stand behind the transformation of the Laplacian terms present in the HF and Poisson Eqs. (9), (11), (10) to the PS representations of these terms in Eqs. (30) and (33).
The PS solution of the HF equation is subject to the Dirichlet boundary conditions: nl (0) = 0 and nl (∞) = 0 which are taken into account in deriving the j discretized representation (33) of this equation. Thus, the Kato cusp condition, satisfied by the exact HF orbitals with the s symmetry at r = 0 due to the −1∕r singularity of the Coulomb potential, is not directly imposed on the approximate orbitals in the PS method, unlike in the finite-difference methods that involve outward integration of the HF equation starting at a finite radius very close to r = 0 with the initial condition set using the cusp condition [5,13]. Therefore, the ratios q n0 and q tot obtained with the PS method are not exactly equal to −Z though their deviations from this value are expected to rapidly decrease as the PS solution becomes more accurate with increasing the number of grid points N.

PS solution for closed-(sub)shell atoms and the iteration method
The HF equation is solved with the aid of the derived PS formalism for the atoms with closed shells or subshells, from He to No. The reported results include the orbitals and their energies, the total and exchange energies, as well as the quantities that are the subject of the virial theorem and the Kato cusp condition.
The radial wavefunctions nl obtained as the PS solution of the HF equation are further modified at each iteration step k by the linear mixing of input and output values of the HF orbitals at all grid points r = r i , with a fixed mixing parameter . Here the output (k)out nl is the solution of the PS HF Eq. (33) including the Fock matrix calculated with the input (k)in nl obtained in the previous iteration step. The speed at which the iterated HF orbitals converge strongly depends on which is tested by choosing = 0.1, 0.2, … , 2.0 and calculating |q vir − 2| in consecutive iteration steps (see Fig. 1). The optimal value of and the corresponding minimal number of iteration steps k min after which the solution approaches its convergence limit are shown in Table 1. Other iteration schemes, like the Pulay method [63], could be also be applied in order to speed up the convergence. However, this is actually not needed since the presented PS method of numerically solving the HF equation is very undemanding computationally. Indeed, the calculations with N = 100 grid points and 200 iterations, run as a Fortran programme on a lowend processor (Intel Core i3-2310M, 2.10 MHz, installed in a notebook), require the CPU time of no more than 3 s even for the heaviest atoms.
This estimate of the CPU time is valid for the standard version of the developed programme which assumes the double precision (DP), i.e., the 64-bit floating-point format and is used to obtain most of the reported results. The code makes use of the Intel Math Kernel library (MKL) [64] which contains linear-algebra routines from the LAPACK and BLAS libraries [65] and is optimized for the applied Fortran compiler and hardware (processor). In particular, the solution of the discretized Poisson equation (27) is done with the DPOTRF and DPOTRI routines (and other routines Fig. 1 Absolute difference |q vir − 2| versus iteration step k for the Ar and Yb atoms, obtained using the double-precision (DP) and quadruple-precision (QP) floating-point formats. The virial ratio q vir = −E pot ∕E kin , with E kin found from Eq. (42), is obtained with the N = N min grid points and various values of the mixing parameter (two left panels), and with the value of = 0.9 optimal for Ar and various N ≤ N min (right panel). If |q vir − 2| obtained with the DP becomes less than 10 −13 after k min = k min (N, ) ( Table 1) iteration steps the parameter = 0.1 is used in steps k ≥ k min + 10 to reduce oscillations of q vir . In the QP calculations, = 0.1 is used after 45 steps at N = 90 for the Ar atom Table 1 The number of grid points N = N min for which the converged value of |q vir − 2| becomes less than 5 × 10 −14 , the optimal value of the mixing parameter that gives the fastest convergence in the linear mixing scheme and the respective minimal number of iteration steps k min after which |q vir − 2| becomes less than 10 −13 The last column shows N Kato min which is the minimal N at which the Kato index k n0 for the outmost occupied s orbital becomes less than 10 −10 for He and Be atoms and 10 −11 for other atoms. The results are obtained using the DP floating-point format Journal of Mathematical Chemistry (2020) 58:1571-1600 used by them) and the routine DSYEVR is used to find the numerical solution of the HF eigenvalue problem with the symmetric matrix defined in Eq. (33). Since the numerical precision of the DP floating-point calculations is limited to 15-16 digits, this also limits a maximum accuracy of the obtained results, with an expected loss of one or a few significant digits due to roundoff errors.
To increase the accuracy beyond this limit, the programme has been converted (with suitable compiler options) to its alternative version which uses the quadruple precision (QP), i.e. the 128-bit representation of the floating-point numbers.
Since the MKL library package bundled with the compiler does not provide such an extended precision, the original LAPACK and BLAS libraries [65] have been used instead after compiling them with options automatically providing the QP for all floating-point variables and constants. These modifications resulted in a significant increase of computation time, by around two orders of magnitude, due to the extended precision and the use of unoptimized library routines. However, the CPU time required, using the above specified computer, for the QP calculation with 130 grid points and 200 iterations does not exceed 320 s for heaviest atoms, like No, while the accuracy of the results is improved by many orders of magnitude, as it is presented below. Thus, the LAPACK and BLAS libaries compiled in the QP have been successfully applied in the present calculations although these two libraries are officially dedicated only for the use in the single (32-bit) and double precisions.
The constant L used in the mapping (25) is found to give highly accurate results for a wide range of its values, from 0.3 a.u. to 2.0 a.u.. The applied value L = 1.5 a.u (1.2 a.u for Ra, Rn, No) is optimal for the accuracy of the total energy and the satisfaction of the virial and Kato-cusp theorems.

Total energies and the virial theorem
The overall accuracy of the HF solution can be conveniently examined by calculating the deviation of the virial ratio q vir = −E pot ∕E kin from its theoretical value 2. The absolute difference |q vir − 2| decreases exponentially in consecutive iteration steps k and converges to values less than 10 −13 in the calculations with the DP, provided that the number of grid points N is sufficiently large, as seen in Fig. 1. This limit is reached after a moderate number of iterations k min : 14 for He, 23 for Ca, 33 for Hg, 42 for Yb, and between 19 and 40 for other atoms (Table 1), using the linear mixing method. A very similar number of iterations is needed for the convergence of the total energy: its value converges to 13 digits after k min iteration steps. The energy values with 14 significant digits are achieved by lowering the mixing parameter to = 0.1 for k > k min + 10 which reduces very small oscillations of E tot during iteration. Averaging of E tot over many iteration steps ( k min + 20 ≤ k ≤ k min + 120) finally improves its convergence to 15 significant digits, the ultimate numerical accuracy in computations using the double-precision floating-point format. Thus maximally converged values of E tot = E tot (N) are obtained for each number of the grid points N.
However, the HF solutions and the converged values of E tot (N) can have much lower accuracy if the number of the grid points N is too small. This is illustrated in Fig. 1 (bottom panel)) where the absolute difference |q vir − 2| converges to values a few orders larger than 10 −13 if N is not large enough. The numerical errors in the present calculations are the combined effect of the approximate PS representation which becomes more accurate for larger N and the round-off errors growing with increasing N which defines the size of the matrices involved in the solution of the HF and Poisson equations. These counteracting trends are clearly seen in the dependence of |q vir − 2| on N in Fig. 2 and result in a crossover at some N = N min so that for larger N the rapidly decreasing error of the applied PS method is dominated by the unavoidable numerical (round-off) errors. The optimal number N min of grid points depends on the type of atom and varies from around 35 to around 70, being smaller for lighter atoms though this is not a strict trend ( Table 1). The total energy E tot (N) converges after N reaches N min : for N > N min the energy E tot (N) oscillates as a function of N around an average value E av tot and the amplitude of these oscillations (affecting only the 14-th and higher digits of the energy value) slowly grows with N due to increasing round-off errors; see Fig. 3. The energy E av tot is well approximated by the mean value of E tot (N) in the interval N min + 5 ≤ N ≤ N min + 5 + N of the width The optimal fixed values of the mixing parameter, = 0.9 for the Ar atom and = 0.6 for the Yb atom, are taken in iteration steps k < k min + 10 in the DP calculations and k less than 45 and 105, respectively, in the QP calculations, while = 0.1 is used for larger k. The DP and QP results coincide for N < N min (Table 1) Table 2 versus the number of grid points N for the Ar and Yb atoms. The chosen vertical scale directly represents changes of the 14-th digit in the value of E tot . The results are obtained using the DP floating-point format 1 3 Journal of Mathematical Chemistry (2020) 58:1571-1600 N = 30 and it is taken as the final value of the total energy E tot given in Table 2. By examining the convergence of E av tot with increasing N , the accuracy of the so obtained E tot = E av tot ( N = 30) is found to be 14 significant digits for all investigated atoms. This conclusion is supported by the fact that the values of |q vir − 2| are less than 10 −14 if the virial ratio q vir = −E pot ∕E kin is calculated with E pot and E kin obtained using the similar averaging scheme as for E tot = E av tot . The same values of the total energy with 14 significant digits and the virial ratio differing from 2 by less than 10 −14 are obtained if the kinetic energy is found from the alternative formula (44). If the described averaging over N is not applied the energy E tot = E tot (N) shows minute oscillations for N ≥ N min and thus its accuracy in the calculations with the DP is reduced to 13 significant digits (Fig. 3).
The roundoff errors in the calculated energies due to finite numerical precision also manifest themselves as the rapid tiny oscillations of |q vir − 2| (with 10 −15 to 10 −13 amplitude for the DP results) versus both the number of iterations k and the number of grid points N after this deviation reaches the range of its lowest possible values for k > k min and N > N min ; see the flat portions of the plots in Figs. 1 and 2. In addition, the absolute deviation |q vir − 2| of the virial ratio shows another type of non-smooth dependence before it reaches the convergence, i.e., for k < k min and N < N min , with significant lowering of |q vir − 2| in some iteration steps k and for some specific N. In this regime, the dependence of q vir − 2 on k is smooth in fact but this deviation oscillates around 0 during iteration and thus changes its sign at some steps k. Consequently, its absolute value has a local minimum in each of such iteration steps which leads to the obtained irregularities in the otherwise linear decay of ln |q vir − 2| with increasing k as seen in Fig. 1. A similar argument explains the quasi-periodic irregularities in the roughly linear dependence of ln |q vir − 2| on N (Fig. 2) which arise because the deviation q vir − 2 also oscillates around 0 and changes its sign for some N when the number of grid points is gradually increased.
The total energy values, obtained with the DP version of the code, have a very high accuracy close to the precision of the DP floating-point representation (15-16 digits) used in the calculations. The 14-digit values of E tot given in Table 2 are fully reproduced (with the sole exception of the last digit for the Sr atom) in the even more accurate numerical calculations done with the QP version of the code which provides further 11 significant digits of E tot . Indeed, as shown in Table 3 and Figs. 1 and 2 , the absolute deviation |q vir − 2| of the virial ratio found in the QP calculations reaches the limit of 10 −25 to 10 −27 with increasing N. Simultaneously, the values of E tot (N) converge up to 25 or more digits as illustrated for the atom Yb in the lower part of Table 3. Thus, it can be claimed that the total atomic energies calculated with the QP have at least 25 significant digits. This remarkable accuracy is achieved with a moderate number of grid points, from N = 80 for the Ne atom to N = 130 for the heaviest atoms like Ra and No. The number of grid points needed for convergence of E tot in the QP calculations is roughly doubled compared to the respective number (denoted by N min in Table 1) required to obtain the values of E tot with 14 significant digits in the calculations using the DP. The number of iterations needed to reach the convergence of the total energy is also doubled compared to the respective iteration number k min (Table 1) required in the DP calculations. Such doubling results from the fact that the deviation |q vir − 2| decays roughly exponentially with both the number of iterations k and the number of grid points N so that the logarithm log 10 |q vir − 2| decreases roughly linearly with k and N (Figs. 1 and 2) while its minimum value drops almost twice, from around −14 to around −26 , upon extending the floating-point precision from DP to QP.
The obtained values of E tot fully agree with the 9-and 10-digit values of the total HF atomic energies obtained, respectively, by Tatewaki et al. [50] in the numerical HF calculations with a finite-difference method and by Koga [52] in the RHF calculations with a modified Slater orbital basis. The present results also show that the last figure in the 10-digit values of the numerical HF total energies reported by Bunge et al. [49] is slightly incorrect for most of the closed-(sub)shell atoms (usually differs by 1 from the correct figure) while the accuracy of E tot found by Bunge et al. in the RHF calculations with a Slater basis [51] is limited to 8 digits for most of these atoms. The present 14-digit and 25-digit values of the total energy fully agree with its 12-digit values obtained by Saito [21] in the RHF calculations based Table 3 Total energy E tot with the 25-digit accuracy, the virial ratio deviation q vir − 2 , and the HOMO energy HOMO with the 22-to 25-digit accuracy for the closed-(sub)shells atoms from He to No The results are obtained with the specified number N of grid points, using the QP floating-point format, and the virial ratio q vir = −E pot ∕E kin is calculated with E kin found from Eq. (42). The lower part of the on the expansion of the orbitals in the B-spline basis. Lastly, the present value of E tot for the He atom agrees in all 14 digits with the result obtained by Ozaki and Toyoda [25] with a finite-element technique. Although those authors do not give the values of the total energy for other atoms they report that their HF calculations required up to a few thousand of grid points to satisfy the virial theorem with the 13-to 14-digit accuracy for the noble gas atoms. The PS method proves clearly superior in this respect as only 33 to 71 grid points (Table 1) are needed to achieve the same accuracy level using the standard DP format of the floating-point numbers in numerical calculations. Furthermore, it is sufficient to double the number of grid points to obtain the energy values with 25 significant digits though it is then necessary to use the QP floating-point representation to implement the calculations at such an extreme accuracy level.

Orbital energies
The energies nl of the occupied HF orbitals obtained in the DP calculations are also found to converge at N = N min given in Table 1. They show very small oscillations for larger N so that the accuracy of nl can be improved by one order of magnitude by the described averaging over N. The values of HOMO are given for all closed-(sub) shell atoms in Table 2, while the energies of all occupied HF orbitals for four chosen atoms (Be, Ar, Zn and Yb) are given in Table 4. The accuracy of the DP values of the orbital energies reaches 14 digits for innermost shells, and is lowest, with 12 significant digits, for the HOMO energy HOMO . Thus, the obtained orbital energies are much more accurate, up to 7 orders of magnitude, than in the previously reported Table 4 Energies of the occupied HF orbitals in the Be, Ar, Zn and Yb atoms The shown values are found as the averages of nl over the interval of N min + 5 ≤ N ≤ N min + 35 and are obtained using the DP floating-point format calculations. Their values fully agree with the nl values of 10 −6 a.u. accuracy (corresponding to 6-11 significant digits) given by Saito. [21] The occupied orbital energies nl obtained by Bunge et al. in the RHF calculations are slightly higher than their present more accurate values, by 10 −6 a.u. for Be and Ar, and 10 −5 or 2 × 10 −5 a.u. for Zn. The present values of HOMO also fully agree with the 6-digit values of this energy reported for He, Be, Ne, Mg and Ar by Tatewaki et al. [50] Performing the calculations in the QP with the doubled number of grid points N ≈ 2N min provides not only the 25-digit values of E tot but also the orbital energies with the greatly enhanced accuracy -their values converge to 22 to 25 significant digits, as for the HOMO energies given Table 3. The QP values of HOMO fully confirm the 12-digit values of this energy obtained for closed-(sub)shell atoms in the DP calculations (Table 2), apart from the Yb, Hg and No atoms for which the 12-th digit of the respective DP values differs by -1 from the correct digit obtained in the QP calculations.

Electron density at nucleus and Kato's theorem
To verify how accurately Kato's condition is satisfied by the presently obtained HF orbitals with s symmetry the quantities k n0 = | − q n0 ∕Z − 1| are calculated. The magnitude of the index k n0 obtained with the code using the DP becomes smallest, of the order of 10 −10 or 10 −11 , for N significantly larger (even more than twofold) than the grid number N min optimal for the total energy. The number of grid points leading to lowest values of k n0 (and thus the best accuracy of n00 (r = 0) ) grows with increasing n, reaching N ≥ 100 for the outmost occupied s orbital, since more expanded orbitals require more grid points for their accurate representation in the real space with the PS method. The same order of magnitude is found for the index k tot = | − q tot ∕Z − 1| referring to the total electron density n tot (0) at the nucleus. It should be noted that the deviation of −q tot ∕Z from 1 comes mainly from the inaccuracy of the derivative n � tot (0) which converges up to 10 or 11 significant digits. However, a further increase of this accuracy and the corresponding decrease of k tot and k n0 are limited by the finite precision of the applied DP floating-point format, similarly as for the total and orbital energies.
The value of the density n tot (0) , after averaging over several tens of iteration steps for each N, converges up to 12 digits in the DP calculations. The accuracy of the s orbitals and density at the nucleus can be improved to 13 significant digits (except for Yb) by further averaging of n0 (r) over N within the interval of N Kato min + 5 ≤ N ≤ N Kato min + 35 where N Kato min corresponds to the lowest N for which k n0 becomes less than 10 −11 (or 10 −10 for the He and Be atoms) for the outmost occupied s orbital; see Table 1. The so obtained values of n tot (0) and the corresponding Katocusp indices k n0 and k tot are included in Table 5. These indices are of the order of 10 −13 ( 10 −14 for Ne, Ar, Kr and Pd) for most atoms, apart for some atoms (He, Be, Mg, Ca and Yb) with HOMO of the s symmetry, for which k tot is one or two orders larger.

3
The accuracy of the electron density and the s orbitals at r = 0 can be further enhanced by using a larger number of grid points which however this requires performing calculations in the QP to reduce roundoff errors (since these errors grow with increasing N at a fixed floating-point precision). In particular, it is found for the Yb atom that the cusp index k tot decreases down to 10 −22 for N = 200 and the corresponding indices k n0 are of the order of 10 −26 for the 1s, 2s, 3s orbitals, 10 −25 for the 4s, 5s orbitals and 10 −18 for the 6s orbital.
The obtained DP values of the electron density at the nucleus, calculated up to 13 significant digits ( Table 5), show that the five-digit values of n tot (0) found by Koga et al. [66] with the RHF method, are accurate apart from their last digit which differs by 1 from the present values of n tot (0) for some atoms. However, this slight inaccuracy of those previous results is consistent with the magnitude of k tot = | − q tot ∕Z − 1| of the order of 10 −5 given by Koga et al. [66]. It should also be noted that the RHF values of n tot (0) reported by Bunge et al. [51] are less accurate, limited to four significant digits, due to an error in those calculations which was identified Koga et al. [66]. Lastly, the present PS solution determined with the DP provides the electron density at r = 0 with the accuracy better by 5 orders than the 8-digit values of n tot (0) obtained by Saito with the B-spline expansion [21]. Table 5 HF electron density n tot (0) at the nucleus and the quantities k tot = | − q tot ∕Z − 1| and k n0 = | − q n0 ∕Z − 1| which determine the deviation from the Kato cusp condition in closed-(sub)shells atoms from He to No The shown values are found with n00 (0) averaged over the interval of N Kato min + 5 ≤ N ≤ N Kato min + 35 (cf.

Node structure of the HF atomic orbitals
The applied PS method provides very accurate representation of the orbitals in real space throughout the whole atomic volume. In particular, it correctly reproduces the node structure of the exact HF orbitals and avoids artificial nodes arising in the RHF calculations with finite atomic bases. [27] However, there are additional nodes, which are usually considered spurious though they are due to the oscillations in the tails of the exact HF radial wavefunctions [16,67]. Such additional nodes are found with the PS method in all closed-(sub)shell atoms heavier than Mg. This is illustrated in Fig. 4 which shows the occupied orbitals in the Ar atom and their radial nodes as the negative peaks in the plots of log 10 nl . The 2s, 3s, 2p and 3p orbitals have the expected number of n − l nodes, if we disregard distant atomic regions ( r > 20 a.u.) where the values of nl (r) drop below 10 −14 and become completely inaccurate due to the limited precision (here the DP) of the numerical calculations. An additional node in the Ar atom is found for the 1s orbital at r = 1.0932657894 a.u. in its tail region where 100 (r) is very low (around seven orders lower than at r = 0 ). Similar additional nodes are found for the Kr atom, at r = 0.4483419357 a.u. and r = 1.255565714 a.u. for the 1s orbital and at r = 1.506285094 a.u. for the 2s orbital. The latter node, at the same position ( r = 1.51 a.u.), has also been found in the recent numerical HF calculations, with two different finite-difference methods, by Mendez et al. [67] (the 1s orbital in Kr was not discussed by them). The present solution also predicts a node for the 1s orbital in the Zn atom, at the position r = 0.5588720448 a.u. which nearly coincides with the node at r = 0.549308 a.u. found for this orbital by Takeda et al. [68] in the RHF calculations. Hence, this particular node is a feature The values of nl (r) outside the grid points r i are obtained with the PS representation (13) for the function f (x) = nl (x) defined in Eq. (32). The plotted orbitals are obtained using the DP floating-point format of the exact HF 1s orbital and is not the result of its approximate representation in a finite atomic orbital basis as it was assumed by those authors.

Asymptotic variation of the HF orbitals and electron density at large r
The radial wavefunctions nl (r) of the occupied HF atomic orbitals have the common exponential decay with the rate H = √ −2 HOMO governed by the HOMO energy and the leading term of their asymptotic expansion [69][70][71] also includes an orbital-dependent power prefactor r nl and the constant coefficient b nl whose approximate value can be found by fitting asymp nl (r) to the numerical nl (r) at some large r. The exponent nl takes the highest value H = 1∕ H for the HOMO and is lower than H for other occupied orbitals. It is equal to nl = H − 3 for atomic orbitals with the same orbital number l = l H ≠ 0 as the HOMO, while for orbitals with l ≠ l H this exponent is nl = H − |l − l H | − 1 . In the specific case of l = l H = 0 , the value of nl = H − 2(l min + 1) is determined by the lowest nonzero orbital number l min among the occupied orbitals, unless for atoms with s orbitals only ( l min = 0 ), for which each n0 (r) decays with its own rate n0 = √ −2 n0 and n0 = 1∕ n0 . Thus. apart from the last case, the HF orbitals do not follow the orbital-specific exponential decay exp(− nl r) (where nl = √ −2 nl ) characteristic for the solutions of the Schrödinger equation with a local multiplicative spherical potential vanishing at r → ∞ . In particular, such a decay is found for atomic orbitals satisfying the Kohn-Sham equation within the density-functional theory [72], but it is not valid for the solution of the HF equation which includes the nonlocal (integral) exchange operator defined by Eqs. (3)-(5).
The very high accuracy of the PS solution for the HF equation leads to excellent reproduction of the described asymptotic behaviour of the HF atomic orbitals nl (r) . In particular, the common exponential decay is confirmed numerically for the Ar atom in Fig. 5 where the calculated ln | nl (r)| shows an almost linear dependence on r, with the slope similar for all HF orbitals and close to the theoretical value − H = − √ −2 3p ≈ −1.0872 a.u. In fact, all the numerically obtained orbitals nl (r) nearly coincide with their asymptotic form (51) within in a wide interval of r, up to a large radial distance at which the orbitals reach very low values and stop to rapidly decrease with increasing r. This boundary radius, possibly orbital-dependent, determines the spatial region where the PS solution is valid, and the radius increases with increasing the number of grid points N while the corresponding values of | nl (r)| decrease. However, the possibility of such improvement by increasing N is limited by the finite precision of the float-point format used in the calculations.
In particular, the asymptotic dependence of the HF orbitals is very well reproduced for each orbital with a moderate number of grid points ( N = 120 ) in the DP calculations, as shown for the Ar atom in Fig. 5a. In this case, the PS solution starts to break down at a radial distance between 20 a.u. (1s orbital) and 32 a.u. (3p orbital) where the numerical wavefunctions stop to decay exponentially while reaching values of the order of 10 −15 to 10 −13 a.u. (corresponding to ln | nl | from -35 to -30). (51) asymp nl (r) = b nl r nl exp(− H r) These limit values cannot be substantially lowered and, consequently, the validity of the PS solution cannot be extended to larger r by merely increasing the number of grid points beyond N = 120 in the calculations using the DP. This is so because the minimal values of the computed orbitals cannot be smaller than inevitable numerical inaccuracies (roundoff errors) due to the finite precision of the applied floating-point representation (15 to 16 significant digits for the DP).
Thus, to increase the accuracy of the calculated orbitals and by that means to extend the region where the PS solution is valid, a higher precision of the floatingpoint arithmetic is required alongside a larger number of grid points. Indeed, the computations done for the Ar atom with the QP and N = 150 lead to an even more accurate numerical solution which follows the correct asymptotic dependence of the HF orbitals up to r = 50 a.u. In this case, the boundary radius is similar for different orbitals while their respective values range from 10 −29 a.u. (1s orbital) to 10 −22 a.u. (3p orbital), which corresponds to the values of ln | nl | from -67 to -51 as seen in Fig. 5b. If the number of grid points is increased to N = 200 , the boundary of the region where the numerical solution is valid is extended to r = 60 a.u at which radius all the Ar occupied orbitals become less than 10 −27 a.u.
The QP solution facilitates an even closer look at the large-r variation of the HF orbitals determined with the PS method. In particular, it is revealed for the Ar atom that these numerically obtained wavefunctions reproduce very accurately not only the theoretically predicted exponential decay of the HF orbitals but also the power prefactors r nl present in their asymptotic expansion asymp nl (r) . Indeed, the quantities S nl (r) = ln | nl (r)∕ asymp 3p (r)| plotted in Fig. 6 almost exactly follow linear dependencies on ln(r) for 20 a.u. ≤ r ≤ 45 a.u. , which have the form of −2 ln r + c nl for the s orbitals and −3 ln r + c 2p for the 2p orbital, while the relation S nl = 0 is well satisfied for the 3p orbital (HOMO), all this in excellent agreement with the asymptotic formula (51) if we take c nl = ln(b nl ∕b 3p ).
The asymptotics of the HOMO wavefunction also determines the asymptotic variation of the total radial density, which is thus given by where l H is the orbital number of the HOMO and b H is the coefficient b nl in its asymtotic expansion (51). It is found that the theoretical relation (52) is very well reproduced by tot (r) obtained in the PS calculations with the DP in a wide interval of r outside the occupied shells. In particular, this is true for radial distances up to 32, 42 and 58 a.u. from the nucleus for the Ar, Zn and Yb atoms, respectively; see

Conclusions
The HF equation for closed-(sub)shell atoms takes a simple form and can be easily solved numerically when one-electron wavefunctions are discretized on a radial grid with the PS method. Its solution involves standard problems in linear algebra like matrix inversion and diagonalization of a symmetric matrix. The boundary conditions are accounted for in the exact way in the derived formalism and no special treatment is needed to determine the numerical solution at points close to the nucleus and very distant from it. The PS solution of the HF equation allows for immediate calculation of various radial integrals (like those defining the total energy and its constituents) by using the highly accurate GLL quadrature based on the PS grid points. The results obtained in calculations using the DP floating-point format are of very high numerical accuracy which is achieved for grids including just a few tens of radial points In particular, the values of the total energy and the electron density at the nucleus are determined with 14 significant digits while the HOMO energies are found with 12-digit accuracy. The PS calculations also yield extremely accurate radial wavefunctions which are valid up to distances of few tens of atomic units from the nucleus and remarkably well reproduce the specific asymptotic variation of the HF orbitals, including their common exponential decay. The excellent accuracy of the obtained results is enhanced even further, to 25 significant digits for the total energies and 22 to 25 significant digits for the orbital energies, by extending the precision of the floating-point format to the QP while still using a moderate number of about 100 grid points.
Summarizing, the simple formalism and its straightforward numerical implementation combined with the unsurpassed accuracy of the calculated quantities make the PS method superior to other approaches to the numerical solution of the HF equation for closed-shell atoms. Further investigations are needed to extend this approach to open-shell atoms, described either with orbital-specific Fock operators or unified coupling operator in the restricted HF method.