Numerical variational solution of hydrogen molecule and ions using one-dimensional hydrogen as basis functions

The ground state solution of hydrogen molecule and ions are numerically obtained as an application of our scheme to solve many-electron multi-center potential Schrödinger equation by using one-dimensional hydrogen wavefunctions as basis functions. The all-electron sparse Hamiltonian matrix for the given system is generated with the standard order finite-difference method, then the electronic trial wavefunction to describe the ground state is constructed based on the molecular orbital treatment, and finally an effective and accurate iteration process is implemented to systematically improve the result. Many problems associated with the evaluation of the matrix elements of the Hamiltonian in more general basis and potential are circumvented. Compared with the standard results, the variationally obtained energy of H2+ is within 0.1 mhartree accuracy, while that of H2 and H3+ include the electron correlation effect. The equilibrium bond length is highly consistent with the accurate results and the virial theorem is satisfied to an accuracy of −V/T = 2.0.


Introduction
One-dimensional hydrogen atom (1D H) [1][2][3][4] has found many applications in high magnetic fields [5], semiconductor quantum wires [6], carbon nanotubes [7], and polymers [8]. Nevertheless, the application of the 1D H function to molecular calculation has not yet much investigated. Unlike the 1D H functions, three-dimensional Gaussian type orbitals [9][10][11][12] and Slater type orbital (STO) [13,14] are more popularly used functions. The former uses the analytical and the recursive integration scheme of the two-electron integrals calculation, while the latter has an extensive analytical energy expressions table [15] as well as basis size efficient. In fact, the novelty of a computational procedure is seen from the use of 1D Cartesian functions. The 1D functions can be considered as more general basis functions because the 1D functions can generally reproduce the result of the 3D counterparts but not otherwise, as in the Gaussian orbitals case where e −r 2 is composed of three 1D Gaussians in x, y, z. If the 1D H function is to be explored for molecular calculations, the studies need to be proceeded within Rayleigh-Ritz variational framework. Moreover, the studies will be numerical in nature due to the inseparability of the Coulomb attractions and repulsions in this coordinate system which make analytical integrations infeasible. As the general solution, the trial wavefunction can be expressed as the linear sum of the 1D H functions. Then the normalizability of the three-dimensional trial wavefunction relies on the variational coefficients which need to keep the wavefunction small at large distances. The energy is allowed to be zero or positive like the boundary condition of the 1D free particle [16]. This particular step is different than the common procedure which usually constructs the wavefunction as a product of wavefunction in each dimensions like the particle in a box problem. The available additional variational space in the use of the linear combination of the 1D H functions is capable of capturing more information about the behavior of the electrons in each dimension.
Subsequently, the optimization process is implemented to provide reasonable approximate ground state solutions of the many body Schrödinger equation.
The linear combination of atomic orbital molecular orbital [17] has been successfully applied for many-body Hamiltonian molecular calculations. In the early H 2 calculation, the molecular orbitals [18][19][20][21] treatment performed by Coulson [17] approximated the wavefunction of H 2 molecule as the sum of 1s STO electron atomic orbitals. The treatment results in the analytical expressions of most energy components of both electrons [15]. Yet the development is rather slow due primarily to the multicenter integral problem [22]. The great progress is made through the one-electron effective Hamiltonian in which many numerical methods have been developed [23][24][25]. For example, in order to overcome the multicenter integral computational difficulty, Becke [26] introduced the partitioning scheme to divide an integral over the entire molecular volume into a sum of atomic orbitals where the weights are defined in terms of atomic cell function. Nonetheless, it is desirable to have a molecular calculation technique which solves the original many-electron molecular Schrödinger equation with more general basis functions in the form of the 1D H functions and which does not lead to the integrals of great difficulty with unnecessary partition scheme required.
In section 2, we present a numerical variational and an optimization method for the molecular calculations. The implemented multidimensional numerical integration technique makes it possible to use more general basis function and to avoid the integration difficulties inherent in ordinary molecular calculations where nonspherical potentials are present. As an example, the method is applied to the energy calculation of hydrogen molecule and ions by using the 1D H functions. In section 3 the detailed calculations of the energy, the equilibrium bond length, the virial ratio, and the density of hydrogen molecule and ions demonstrate the feasibility and the flexibility of our method which is considerably easy and effective without the need of the quadrature formula for the integration of the spherical harmonics and the need of the standard quadrature for the radial integrations. The method [27,28] is free from the decomposition of molecular functions into single-center components and is of potential to be developed to a certain accuracy and to larger size molecules.

Theory
We consider the variational solution [16] of the molecular Schrödinger equation with an all-electron potential so that, in atomic units and the Born-Oppenheimer approximation and non-relativistic limit, the Hamiltonian expression for nuclei a, b and electrons i, j iŝ although it is important to note here that our method depends in no way on the specific form of the potential.
In the linear variational framework, the approximate solution of equation (1) for one-electron systems is sought by expanding the eigenfunction as ψ = j c j χ j .
where the coefficients c j are parameters to be determined variationally and the basis functions χ j are the 1D H functions centered on the corresponding nuclei, for example the 1D H [29] function of electron 1 centered on x-coordinate of the nucleus a is where the principal quantum number n, the generalized Laguerre polynomial L 1 n−1 , and the wavenumber k related to the corresponding 1D H energy E x as k = √ 2mE x / . To describe the ground state of, for example H 2 + , a trial variation function which uses as a basis three 1D H functions of the n = 1 states centered at each nucleus and follows the MO treatment can be constructed. The trial wavefuction has six variational coefficients to be determined We note that each of the basis functions in equation (4) and their linear combination are generally not normalizable in 3D space since they do not vanish at infinity particularly in the other two dimensions which the 1D H function does not depend explicitly on. Therefore, to make the probability amplitude finite the variational coefficients need to keep the wavefunction small at large distances.

Constructing the matrix H ij and S ij
and solving the familiar generalized eigenvalue equation for the minimum eigenvalue E (0) and the corresponding eigenvector ψ (0) , we obtain the approximate ground state energy and the variational coefficients, respectively. In order to minimize further the error in the energy, the well-known iterative process of the residual minimization method-direct inversion in the iterative subspace (RMM-DIIS) [30][31][32] is used, started with the evaluation of the residual vector δ (0) = Ĥ − E (0) ψ (0) . Then a trial step ψ (1) = ψ (0) + C δ (0) along this direction is performed, and the new residual vector δ (1) is calculated. A linear combination of the initial ψ (0) and the new trial wavefunction ψ (1) is sought to minimize the norm of the linearly assumed residual vector which is equivalent to finding the lowest eigenvector or eigenvalue from the secular equation Overall in the m + 1 trial step, ψ (m+1) which equals to ψ (m) along the direction δ (m) and a new residual vector δ (m+1) are added to the iterative subspace. If the minimization goes correctly, the trial step C becomes considerably small and the process of iteration is terminated once the difference between the successive eigenvalue E in equation (7) is smaller than 10 −5 a.u. Depending on the number of electrons n, we need to generate uniform orthogonal 3n-dimensional grids of equal spacing h. In the case of H 2 , for example, the points are defined in a finite domain by x 1i , y 1j , z 1k , x 2l , y 2m , z 2n . The finite-difference method is implemented to create approximately the H 2 kinetic energy operator matrix representation where − 1 2 ∇ 2 is expanded with the standard order (c ±1 = 1, c 0 = −2) second derivative finite-difference expressions although higher order may be used [33]. So, we approximate ∂ 2 ψ/∂x 2 1 at x 1i , y 1j , z 1k , x 2l , y 2m , z 2n by assuming the wavefunction ψ can be approximated accurately by the Taylor power series expansion in h and by using 'three-points' expressions. The corresponding potential energy operators are sampled on the grid points After all quantities are put in the grid representations, the multidimensional numerical integration over the Cartesian grids is performed to obtain the matrix elements. For example, integrals over six-dimensional space to obtain the H 2 Hamiltonian matrix elements in the 1D H functions basis become sum over grid points m, n of grid spacing h in the position of grid point r m , r n . The computation looks like a simple row vector-matrix-column vector product times the volume per grid points: The integration scheme equation (11) offers an alternative from the well-known analytical Gaussian product theorem or STO tables. Unlike the plane wave basis [30,[34][35][36][37] integral evaluations which require reciprocal space, the integrals with the 1D H function are easily performed in real space. The accuracy of the integration scheme above largely depends on the number of grid points and the grid spacing size.
The computational scaling in each iteration is the same as that of the RMM-DIIS [30]. The evaluation of the residual vector which is a matrix-vector multiplication of the Hamiltonian on a trial function scales as O N 2 for N system size. Then O N 3 operations such as the exact diagonalization and orthogonalization are negligible for small matrix. Hence, the total scaling depends on the number of iterations required to converge times the scaling of each iteration.

Results and discussion
Our goal is to apply the scheme in the calculation of the ground state energy, the bond length, the virial ratio, the wavefunction, and the density of hydrogen molecule and ions by using the 1D H functions as the basis. The energy of H 2 + is calculated with the grid spacing size h = 0.03 and the number of the grid points N g = 400 3 , while the energy of H 2 and H 3 + are evaluated with h = 0.43 and 0.42 respectively for the N g = 26 6 . The compromise between the computational convenience, the energy accuracy, and the virial ratio determine the choice of the number of the grid points and the grid spacing size. In the present setup, we use considerably more number of grid points and finer grid spacing size than those of the common density functional theory (DFT) codes [38,39]. Furthermore, the 3n-dimensional box size ensures that the studied system is fully confined. The eigenvalue E (0) obtained from equation (6) is generally inferior to that of other popular basis functions, which is understandable since the multicenter nonspherical potential is inseparable in Cartesian coordinates. The variational coefficients of the H 2 + trial wavefunction in equation (4) are given in table 1 where the coefficients in the symmetry axis c 3 and c 6 are more dominant than the other axis coefficients and tend to cancel each other with their opposite sign. This is reasonable because for large values of x and y but with z near the two nuclei, the χ a (z 1 ) and χ b (z 1 ) contributions to the trial wavefunction is still large. In order to keep the wavefunction small in these regions, the two terms have to cancel each other out hence the opposing coefficients. The other coefficients behave as expected. The coefficients c 1 and c 4 are considerably small although they do not cancel each other, while c 2 and c 5 are small and opposing each others. Before this stage, the major difficulty in applying the Rayleigh-Ritz method is the evaluation of the Hamiltonian matrix elements, particularly with the multicenter potential and general basis function, but in our scheme as the kinetic operator construction is performed in the Cartesian system and the potential is put on the grids, convenient multidimensional integration like matrix-vector multiplication can be easily carried out, resulting in quantities to be used in the next steps. Additionally, our numerical method has another advantage for inexpensive preliminary studies with a small number of sampling points in order to reveal the gross features of the energy.
To do the correction, the variational coefficient C in equation (8) is relatively large at several first iterations. Subsequently the value of C becomes small and displays the 'see-saw' pattern, a typical minimum discovery process such as in the steepest descent or in the conjugate gradient method until convergence is reached. Its behavior at each iteration is shown in figure 1.
After correction, the typical convergence energy of the hydrogen molecule and ions obtained by using the 1D H functions is shown in table 2. The energy of H 2 + is already very good, to within 0.0001 a.u. of the exact energy -0.6026 a.u., and more improved than table 1 in reference [23]. The result is explainable if we see table 3 where the comparison between our scheme and the Hartree-Fock method performed with large basis set cc-pVQZ for the evaluation of the energy components integrals are presented. The use of the fine grid spacing size in H 2 + accurately matches the integrals evaluation of HF/cc-pVQZ to two decimal places. The kinetic and the potential energy satisfy the molecular electronic virial ratio −V/T [40] to an accuracy of 2.0. Moreover our method approaches the true energy in variational behavior from rigorous upper bound to the converged energy, unlike the 'anti-variational' behavior shown by the previous method [23], as seen in figure 2 where the iterative energy E of equation (7) and the virial ratio are plotted at each iteration. As expected, both undergo correction from the early iteration until the convergence is achieved.
Similar to H 2 + , the behavior of the coefficients of the H 2 trial wavefunction keeps the wavefunction small at large distance either by having almost zero values or same yet opposite signs values. The correction coefficient C of H 2 is now two orders of magnitude larger than the previous one although they still share   [42] where the orbital exponent was optimized to 1.32. b Exact solution [43][44][45]. c Variational method [46,47]. d Configuration interaction expansion method [48]. the same behavior. This is also the reason the time taken for H 2 + to reach convergence is longer than that of H 2 as shown in table 2 where the total computational time for the given systems to converge is presented. Smaller correction coefficients in H 2 + cause more loops to converge hence longer computational time. As seen in table 4, the choice of the grid specification for H 2 at the present level development yields the integrals evaluation of the kinetic and attraction energy which match those of the HF/cc-pVQZ to one decimal place although we need to treat operator H with Coulomb singularities at each nucleus with the uniform point density. The results of H 2 energy components also imply the electron correlation effect is included although no explicit r 12 dependence shown in the trial electronic wavefunction equation (2) which has now thirty six variational coefficients to be determined as the product of the electron 1 and 2 H 2 + trial wavefunction based on the MO treatment subject to symmetry and spin constraints. Avoiding the term r 12 comes with the cost of more set of generalized Laguerre polynomials to describe the electronic wavefunction, in parallel with the conclusion obtained earlier [49]. As commonly known, the HF underestimates the electron correlation effect, shown with their electron-electron repulsion component in   [50] where the former is seen from the use of the atomic and ionic parts of the MO constructed wavefunction whereas the latter is accounted by the explicit z dependence of the wavefunction. The kinetic and potential energy component keep on satisfying the molecular electronic virial theorem to an accuracy of 2.00 while that of HF and DFT have slight deviation as seen in table 4. If more terms or some modifications are not present to approach the correct asymptotic form of H 2 , we anticipate our result fails at dissociation limit since our wavefunction contains polynomial ∼xe −kx . As shown in table 5, the equilibrium bond length calculated with the present scheme is in well agreement with the results of the standard methods and the accurate one.
Another major advantages of our method is the ease with which wavefunctions and charge densities are obtained. If the Slater exchange in DFT framework is used, the obtained wavefunction thus the charge densities can be employed to generate self-consistent local and nonlocal potential in an iterative calculation [25]. Figure 3 shows the electron density distribution obtained from the square of the wavefunction along the molecular axis of hydrogen molecule H 2 developed at several choices of iteration. In the early iterations  the density is asymmetric with the values of the cusps minor compared with those of higher iterations. The DFT based density is also presented to show the difference between the computed density and the standard method density.
The residual plot which shows quantitatively the absolute difference between the computed density and the DFT based density is given in figure 4. As seen, the computed density is still less symmetric compared with that of the DFT one. The density in the negative z-axis is less developed than that in the positive, observed from the difference in region A, B, and C. The correlation energy contributes to the difference in region A, while the cusp at the nuclei position causes the difference in region B. The vanishing density away from the nuclei contributes to the residual in region C. Figure 5 shows the surface plot of the electron density of hydrogen molecular ion H 2 + at the equilibrium bond length. The density can reproduce the familiar singularity or the cusp at the nucleus position without any difficulties to recognize the nearly vanishing one far from the nuclei.

Conclusions
The combination of the numerical Rayleigh-Ritz variational technique and the RMM-DIIS optimization method is presented to solve the molecular Schrödinger equation. The scheme is capable of applying one-dimensional Cartesian function as the basis, and is specialized for obtaining the ground state of many-electron hydrogen molecule and ions. Detailed results of the energy, the equilibrium bond length, the virial ratio, and the density are given. At the present level of development, the energy of H 2 + is within 0.0001 a.u. accuracy from the exact solution, while the electron correlation effect is included in the case of H 2 and H 3 + . The virial ratio is satisfied to an accuracy of at least 2.0. The results of the calculation indicate that the scheme is flexible and easily implemented without any partitioning of molecular systems into single-center terms and without any Fourier transform required. Furthermore, the results demonstrate potential development of the scheme to calculate larger size molecules with improved accuracy and with the electron correlation effect included. This can be proceeded with the use of simply more powerful computing resources or the technical improvements in the underlying theories such as HF or DFT framework. Another development includes the application of our method to do the excited states calculations. Both the linear variational principle and the RMM-DIIS optimization method are capable of performing such calculations, but with some modifications in the latter, for example by using larger iterative subspace.