Solution of two-electron Schrödinger equations using a residual minimization method and one-dimensional basis functions

In conventional self-consistent field (SCF) calculations such as density-functional theory (DFT)^1,2 and Hartree–Fock (HF),^3,4 the two-electron old problem of the helium atom and its isoelectronic ions is proceeded with separating the two-electron problem into a set of accessible single-electron self-consistent eigenvalue equations. Consequently, the many-body effects of electron correlation are excluded in HF. To describe such effects, DFT necessitates exchange–correlation functional.^5–7 These approaches with separation of particles are different from the non-SCF accurate variational calculations of energy levels of the helium and its isoelectronic ions, which have been effectively performed with either Hylleraas choice of basis set or hyperspherical coordinate method.^8–14 The corresponding accurate variational calculations have tackled the Schrödinger equations of the helium and its isoelectronic ions in their full dimensionality and as full two-electron problems.

In this work, to solve the time-independent Schrödinger equation of the helium atom and its isoelectronic ions, we follow neither the former nor the latter approaches. We employ the one-dimensional basis functions (1D functions), so the separability ansatz by the Cartesian components is assumed rather than by the particle, as stated in recent papers.^15–18 The primary significance of the 1D function approach over any such particle-separability approaches lies in the electron correlation, which can be exactly treated. The approach can include many-body effects directly; hence, more accurate electron correlation, as we shall demonstrate when the repulsion energy comparison is performed. Another key benefit is that the basis that is expanded in all Hilbert space directions at once. In the present work, we use the real-space grid representation of the Hamiltonian and the wavefunction. The former is represented as sparse matrices and the latter as vectors. This brings to the next benefit that allows well-known iterative methods, such as the residual minimization method in the direct inversion in the iterative subspace (RMM-DIIS),^19–22 to be implemented. Iterations by the repeated applications of the Hamiltonian operator $Ĥ$ on a trial energy wavefunction $ψ$ are conveniently carried out since the action of an operator on a wavefunction is then implemented numerically as a straightforward matrix–vector product. Unlike the plane wave basis¹⁷ integral evaluations that require another type of approach, the integrals are easily performed in the real-space basis. The needs of integral evaluation for matrix elements are met with the presence of a convenient numerical integration scheme as a row-vector matrix column-vector product times the power of 6 of the grid-spacing size. The residual norm is progressively made small until the convergence is reached. Hence, our approach is similar to the previous accurate variational calculations for helium and its isoelectronic ions where SCF is not needed. Furthermore, the behavior of numerical convergence is known.

In the present work, we calculated the ground state of the helium atom and its isoelectronic ions but mostly concern with the novel algorithmic development necessary to implement the pair-wise electron approach. Accordingly, this straightforward calculation is considered using primarily Cartesian real-space grids; hence, neither the coordinate system nor the basis set is customized. The chosen crude basis has the wavevector itself as 6D and the corresponding full basis representation Hamiltonian matrix as effectively 12D. If the number of basis functions per dimension N is 30, the calculations will require storage of vectors with N⁶ = 30⁶ $≈$ 10⁸ and matrices with N¹² = 30¹² $≈$ 10¹⁷ elements for their explicit representations. Despite the simple choice of basis, we are able to obtain the two-electron Schrödinger solution exactly using around 10¹¹ (100 GB) of RAM. This suggests that the explicit storage of such vectors and matrices is generally not necessary.

Our two-electron approach brings significance for many-electron problem. First, all the relevant quantities of the exact two-electron electronic structure are available within our approach including kinetic energy, ionic energy, and electron–electron repulsion energy. Second, the obtained wavefunction stores information to plot the radial correlation function and common 1D single-electron functions. Third, if the present approach can reach to four electron calculations, a large electronic structure problem may be considered as solved to reasonable accuracy via a corresponding contracted Schrödinger equation.²³ Specifically, the reduced density matrix N-representability conditions can be formulated directly for four electrons.^24,25 Even if this is not going to happen, the explicit two-electron Schrödinger solutions still serve as a useful and accurate approximation for the many-electron problem through the second-order density matrix.²⁶ 

In atomic units, the nonrelativistic Schrödinger wave equation of the helium atom and its isoelectronic ions is given by²⁷ 

where r_i, i = 1, 2 is the radial distance of the ith electron from the nucleus, E is the energy of the system, Z is the nucleus atomic number, and $ψr1,r2$ represents the two-electron wavefunction.

To solve a two-electron problem in Eq. (1), separation by particle, by space-vector component, or by both can be considered. Nonetheless, our strategy is to use the space-vector components for the sake of simplicity. Consequently, we need a set of optimal bases. For the 6D problem, a wavevector that is separable by the component but not necessarily by the particle has the general form of

We assume that $ϕx1,x2$ are eigenstates of some 4D single-component two-electron effective Hamiltonian $ĥdx1′,x2′;x1,x2$¹⁷ The efficiency of this basis can be seen from the formation of the strongly separable approximate Hamiltonian $Ĥ0$⁠,

where $Î$ is the 2D single-component two-electron identity operator and $ĥx′;x$ is the 2D single-component two-electron effective Hamiltonian. Then, $ψr1,r2$ in Eq. (2) yields the lowest eigenfunction of a separable approximate Hamiltonian where the separability is component-wise. In the present work, we choose a 2D helium atom and its isoelectronic ions’ wavefunction as the spatial components $ϕx1,x2$⁠. A 2D single-component two-electron of helium and its isoelectronic ions’ wavefunction $ϕx1,x2$ are obtained by solving their corresponding 2D wave equation given by

where β = 0.000 01 has been added to avoid the infinite solution at x_i = 0, i = 1, 2.²⁸ 

To represent the Hamiltonian in Eq. (1), we generate uniform orthogonal 6D grids of equal spacing h. Each side of the uniform 6D hypercubic grids is of length 2r₀. In principle, there are two convergence parameters: (a) the maximum range of each position coordinate r₀ and (b) the number of grid points per Cartesian component N. For choosing r₀, we follow the work of Hawk and Hardcastle²⁹ who stated that the best value for r₀ is inversely proportional to the atomic number Z. Additionally, we choose r₀ such that the energy is lowest because the optimal choice of r₀ corresponds to the numerical integration scheme. r₀ is fixed depending on the given systems throughout so that N is the only convergence parameter. The choice of r₀ for the helium atom and its isoelectronic ions is presented in Table I. As N increases, the grid-spacing size h, which is defined as $h=2r0N−1$⁠, decreases proportionately.

The great advantage of this approach is that all the potential energy contributions to the Hamiltonian are represented as diagonal matrices as a consequence of a certain quadrature approximation.³⁰ At $x1i,y1j,z1k,x2l,y2m,z2n$ grid position, the potential energy operators are sampled as

To avoid effectively Coulombic singularity, even number of grid points N is chosen so a grid point at the origin can be circumvented. The electron repulsion singularity is handled by restricting the two electrons not to occupy the same space. To evaluate ionic energy integrals more accurately without requiring very fine grid-spacing size, the centers of the two electrons are slightly displaced in the x direction by the grid-spacing size distance. At the continuous limit or h → 0, the two centers will coalesce.

The finite-difference method is implemented to represent the kinetic energy operator matrix representation. The Laplacian is expanded with the standard order $c±1=1,c0=−2$ second derivative finite-difference.³¹ Hence, we approximate $∂2ψ/∂x12$ at $x1i,y1j,z1k,x2l,y2m,z2n$ by

assuming that the wavefunction ψ can be approximated accurately by the Taylor power series expansion in h and by using “three-point” expressions.

In the grid representation, the 6D two-electron representation of a trial wavefunction takes the form

where N is the number of real-space grid basis functions for each component of each particle and coefficients c_ijklmn. This results in 6D hypercubic grids, so the required storage for the two-electron wavefunction scales as N⁶.

In the grid representation, the operators of all quantities are represented as matrices and the wavefunctions as vectors. Then, the convenient multidimensional numerical integration over the Cartesian grids can be performed to obtain the matrix elements. The action of an operator on the wavefunction is simply a straightforward matrix–vector product. Hence, the accuracy of the implemented integration scheme largely depends on the number of grid points and the grid-spacing size.

After obtaining the electronic structure Hamiltonian, we can now move to calculate the lowest energy level and eigenfunction. In this situation, iterative Krylov subspace methods such as residual minimization method-direct inversion in the iterative subspace (RMM-DIIS)^19–22,32 are implemented since they are faster than traditional methods, particularly if the computational cost of a matrix–vector product is considerably low, as is the case here. First, the residual function $δ0r1,r2$ is obtained from the proposed wavefunction $ψ0r1,r2$ and its energy E₀ by using the relation

An improved trial wavefunction $ψ1r1,r2$ can be obtained from the linear combination of the residual function $δ0r1,r2$ and the initial result $ψ0r1,r2$⁠,

The trial steps C₀ and C₁ minimize the norm of the residual vector. They can be determined from another generalized eigenvalue problem in the m step,

where the residual vector matrix $δij=δiδj$⁠, the overlap matrix in the iterative subspace $Sij′=ψiSψj$⁠, the optimized energy E, and the trial steps C_j, j = 0, m. For the m step, the wavefunction $ψm−1r1,r2$ and the residual vector $δm−1r1,r2$ form the basis set in the iterative subspace. They are combined to give the new wavefunction,

The coefficient C₁ becomes smaller with more steps as commonly found in the optimization process. The iteration is terminated once the difference of the consecutive energies E in Eq. (10) achieves 10⁻⁵ a.u.

The main goal of the calculations here is to compute the repulsion energy of the helium atom and its isoelectronic ions from their corresponding 1D functions. Subsequently, the energy components, the virial ratio, the radial distribution function, and the radial correlation function of the helium atom and its isoelectronic ions are presented. The numerical convergence for the ground state energy of the helium atom as a function of the single-component basis size is shown with its error relative to the accurate results plotted against the grid-spacing size.

Since the two-electron wavefunction [Eq. (7)] is the eigenstate of the strongly separable approximate Hamiltonian $Ĥ0$ of Eq. (3), the eigenvalue E₀ is less superior than that of the eigenvalue of the helium atom and its isoelectronic ions. The optimization process in Eq. (11) is then implemented to obtain more accurate energy. The plots of C₀ and C₁ are shown in Fig. 1. The C₀ is negative, which means that the correction takes place in the steepest descent direction rather than in the steepest ascent. Its value, which is close to 1, ensures that the overlap matrix in Eq. (10) is equal to unity; hence, the correction occurs along the residual vector. Additionally, the orthonormality with the previous wavefunction is maintained. The coefficient C₁ varies largely in the several early iterations, which corresponds to a large gradient vector trial step to find the minimum. The correction leads to the correct ground state and prevents from going to the wrong state. After close to the minimum point, the value of the coefficient C₁ is relatively small, which shows the stability of the algorithm. The zig–zag nature like in the common optimization process is also evident as finer and finer correction occurs until the convergence is achieved.

The energy and the virial ratio as a function of the number of iterations are shown in Fig. 2. The energy obtained using our method approaches the true energy variationally,³³ which means that the corrected energy starts from the upper bound to the converged energy. As refinement of energy occurs, the virial ratio also shows improvement.

After optimization, the typical repulsion energy convergence of the helium atom and its isoelectronic ions for the fixed number of grid points N, which is at 30 per one dimension, is shown in Table I where the grid-spacing size and the maximum distance of each dimension are given. For comparison purpose, the repulsion energy of the helium atom and its isoelectronic ions obtained from Hartree–Fock calculation with very large basis sets HF/cc-pVQZ³⁴ and from the accurate calculation with 27-term Hylleraas wavefunction^8,9 is presented.

As seen in Table I, for all the considered systems, the calculated repulsion energies with the 1D function approach are more accurate than the HF results. Furthermore, the error of the HF for H⁻ can reach 0.1 hartree, while the error of the 1D function approach remains within 0.05 hartree or better. The situation is explainable from the electron–electron interaction that is exactly treated with the use of 1D function approach, which separates the dimensions yet contains electron-pairwise nature. Once considerably finer grid-spacing size and more number of grid points are implemented, the numerical integration scheme yields an accurate result for $1r1−r2$⁠. The fact that the repulsion energy is more improved than the HF calculation, which excludes electron–electron correlation, brings broader significance of the one-dimensional basis function approach for not facing difficulties to account for electron correlation energy such as the separation by particle method encounters.

The electron–electron repulsion energy, the kinetic energy, the ionic energy, the total energy, and the virial ratio of the helium atom and its isoelectronic ions are given in Table II. The increasing error of the kinetic energy, the Coulomb ionic energy, and the total energy as Z increases is understandable due to the heavier systems computed with the same number of grid points. On the other hand, the heavier nucleus requires finer grid-spacing size; hence, a larger number of grid points are needed to achieve the same accuracy. The scaling coordinate of r′ = r/Z where Z is the nuclear charge may be implemented to treat the increasing error of the heavier species, stemming from the same reason above. Nonetheless, the virial ratio is satisfied in our calculation to an appreciable accuracy of 1.99 or better.

If we plot the 1D function error relative to the accurate results vs the grid-spacing size as shown in Fig. 3, a simple power relation is found, seen from the high value of R². This is expected given the employed finite-difference method, the choice of one-dimensional basis functions, and the Coulombic singularity.

The wavefunction $ψmr1,r2$ of the helium atom and its isoelectronic ions contains such vast and elaborate information that the current available programs cannot carry out its entire visualization. For clarity, $ψmr1,r2$ is the two-electron wavefunction in its full 6D from which all single-electron information such as one-electron wavefunction, densities, and radial distribution function can be easily obtained, but not the other way round. Furthermore, the information retrieved from this exact 6D two-electron wavefunction is able to build the radial correlation function which one-electron densities cannot. In the following, several 1D one-electron functions are presented from the 6D two-electron wavefunction to give useful understanding. The first is the 1s1s wavefunction vs radial distance plot and its surface plot characterized by the presence of a cusp at the nucleus position. The second is the radial distribution function that is the probability for finding an electron at a distance r from the nucleus, irrespective of the other electron. The third is the radial correlation function, which shows the probability that the two electrons may be found at a distance r₁₂ from each other, irrespective of the nucleus. Before detailed discussion of the obtained wavefunction begins, the probability density of the He atom extracted from the 6D wavefunction at the last iterations m is shown in Fig. 4 where the DFT-based density is also given. Compared to the DFT-based density, the peak of the optimized wavefunction is relatively lower with the location being slightly more to the right. There is another negligible difference at the position where the density vanishes. The situation is understandable due to the electron correlation effect in our calculated wavefunction. In the one-electron distribution, the electron correlation increases the probability of finding electron within 0.4 bohr and larger than 1.5 bohrs but decreases at the intermediate distances.³⁵ 

The scaled 1s1s wavefunction is plotted as a function of the radial distance in Fig. 5. The cusp at the nucleus position at the origin is reproduced using the spline function. Relative to other lower Z cases, the cusp of Be²⁺ is the highest resulting from its steepest slope although the maximum radius of Be²⁺ for the wavefunction to vanish is the shortest. This is explainable from the potential energy and the kinetic energy values in Table II. Higher potential energy causes electrons to be more localized near higher Z nucleus. To conserve the Heisenberg uncertainty principle, the momentum of higher Z nucleus needs to be more distributed. Hence, the kinetic energy that is proportional to the momentum square is larger for higher Z nucleus. The inset of Fig. 5 shows the surface plot of the 1s1s wavefunction with respect to coordinates x₁ and x₂. The resultant 2D single-component distribution shows the well-known 1s1s wavefunction, which is spherically symmetric. The cusp at the nucleus position is well reproduced with the spline function. The vanishing wavefunction at large distance is seen.

The scaled radial distribution function for the electronic ground states of the helium atom and its isoelectronic ions is presented in Fig. 6. The ground state radial distribution function has a single peak in accordance with the 1s1s character of the ground state. Consistent with the previous results, the peaks are closer to the nucleus for heavier atoms.

In Fig. 7, we present a scaled radial correlation function. The ground state distribution is again singly peaked. From Fig. 7, we can find that for higher Z nucleus, the two electrons persist the tendency to have higher kinetic energy and to move in more localized region. Comparison with the HF result,^36,37 the effect of the electron correlation in the present work is seen in the slight shift to the right of the peak at r₁₂ = 1.03 bohrs and r₁₂ = 0.64 bohr for He and Li⁺ atoms, respectively.

Unlike the standard electronic structure methods of particle-separability-ansatz, we have solved the Schrödinger wave equations of the helium atom and its isoelectronic ions by employing one-dimensional functions to separate components. The single-component two-electron wavefunctions of the two-electron systems are chosen to be the optimal basis for a strongly separable Hamiltonian. The two-electron 6D Schrödinger wavefunction is obtained after further refinements to the true value using a Krylov based iterative method until a stable energy value up to five decimal places is achieved. The approach proves to be simple and efficient, which is applicable to neutral and charged atomic systems. As expected, the many-body effects in the form of electron correlation are more accurately described with the component separation or one-dimensional basis function approach, rather than the single-electron Hartree–Fock or density-functional method. For larger systems that realistically involve hundreds or thousands of electrons and nuclei, the present approach is not a likely candidate to solve their large N-electron Schrödinger equation. This distinct characteristic of our approach, which allows for numerical convergence, is completely unavailable in all standard electronic structure methods.

F.U.R. and Y.P.S. contributed equally to this work.

This work was supported by the NSAF (Grant No. U1930402). We also acknowledge the Beijing Computational Science Research Center, Beijing, China, for allowing us to use its Tianhe2-JK computing cluster.

The data that support the findings of this study are available from the corresponding author upon reasonable request.