A Program for Single-center Expansion in Laguerre-type Orbitals for the Hydrogen Molecular Ion

Abstract

A Fortran program is given for calculating wave functions of the molecular hydrogen ion, expanded in terms of single-center Laguerre-type orbitals. Using this program, the radial quantum number has been extended to 203, and accuracy has been attained in energy of the order of 10⁻⁶ a.u. The electron-nucleus Coulomb integrals are evaluated numerically by applying Gaussian quadratures.

1 INTRODUCTION

In modern quantum chemical calculations, Gaussian-type orbitals (GTOs) [1] are most commonly used as basis functions. To perform accurate calculations for atoms or small molecules, Slater-type orbitals (STOs) [1] are also utilized. GTOs provide an advantage in efficient multi-center 2-electron integral evaluation.

As a result of overlap integrals among these basis functions, however, problems associated with linear dependence [2,3] occasionally prevent accurate computation. This may be a severe impediment when extending the basis toward completeness. In this paper we apply Laguerre-type orbitals (LTOs) [1] which do not suffer from linear dependence issues.

In Hartree-Fock calculations with basis-set expansions, the wave functions obtained often have artificial nodes [4,5,6], which should not appear in the exact wave function. Artificial nodes have been found for GTOs [5,6] and STOs [4,5]. It is of interest whether artificial nodes can occur even for orthogonal basis sets such as LTOs having no overlap integrals. This question motivated us to write the present program involving LTOs.

The Laguerre-type orbitals give rise to a complete orthonormal basis (CONS) [1,7,8]. It is well known that LTOs are complete sets for bound states of the hydrogen atom. By adding increasing numbers of LTOs, wave functions derived using this method approach the exact solution ever more closely.

LTOs nevertheless have shortcomings. The biggest problem, common to all exponential-type orbitals (ETOs), is in evaluating multi-center Coulomb integrals. For this reason, LTOs are applied only to single-center expansion, and actual single-center expansion calculations have been limited to small systems such as the He atom [2,3], H₂⁺ [9,10], or H₂ [11].

The hydrogen molecular ion (H₂⁺) is the simplest molecule and has been studied from the early days of quantum chemistry, beginning with the pioneering works of Burrau [12] and Hylleraas [13]. The exact solution was obtained. Peek’s paper [14] presently provides the most precise value of the exact electronic energy, −1.1026342144949 a.u., for the internuclear distance R_bond = 2.0 a.u. using nonrelativistic theory and the Born-Oppenheimer approximation. This internuclear distance, which is quite close to the equilibrium distance of (R_e = 1.9971933199 a.u.) [15], is used throughout this paper. Bates et al. [16] gave details of the exact wave function including contour maps and the peak value (0.458) at the nucleus. Even today H₂⁺ is an important molecular system used as a testing ground for novel numerical methods and approaches.

Single-center expansion is usually based on the molecular midpoint, but its convergence is slow. Because no basis functions are placed on the nuclear positions, it will be difficult to give satisfactory results for the nuclear cusp condition [17] if a finite number of basis functions is used. This explains the slow convergence.

Single-center expansions with LTOs have been made before for H₂⁺; the first was by Howell and Shull [9], and the second by Kranz and Steinborn [10]. In the latter, the maximum number (Nmax) of the radial quantum number (n) in the basis is 25 and the resulting electronic energy (−1.1022 a.u.) is accurate to the order of 10⁻⁴ a.u. This contrasts with the rapid convergence of the multi-center expansion method (linear combination of atomic orbitals, LCAO) using GTOs. Wells and Wilson [18] quote an energy accuracy of 2 × 10⁻⁷ a.u. with a reasonably compact basis set of [23s 9p 4d 4f 2g] (232 primitive GTOs in total).

In spite of the above shortcomings of LTOs, they are still valuable for their capability to construct a complete orthogonal basis. In testing newly developed methods on the hydrogen molecular ion, it is desirable to lift the limit for the radial quantum number (n) and attain accuracy in energy of the order of 10⁻⁶ a.u.

In calculations with LTOs, the most time-consuming step is the evaluation of electron-nucleus Coulomb integrals. Although analytical integration is possible, there is serious loss of accuracy, as will be described in section 3.1. To overcome this problem a multiple-precision arithmetic package such as MPfun90 [19] or FM [20] is required.

To evaluate integrals efficiently and precisely we wrote a program that implements numerical integration by Gaussian quadratures [21,22,23]. We examined the accuracy of the integrals and found that, even for n = 203, at least 20 decimal digits are maintained as significant figures. Using this program we have successfully extended n up to 203 and obtained energy accuracy of the order of 10⁻⁶ a.u. The dimension of the Hamiltonian matrix is 10404. We provide this program with the present paper.

The program is written in Fortran 2003. It is self-contained, so that no external library programs are needed. The computation is performed with quadruple precision arithmetic based on the IEEE Standard for Floating-Point Arithmetic [24], which is supported by most modern Fortran compilers. As a result, the program has good portability.

Studies of H₂⁺ continue and recent theoretical work can be found in the literature [25,26,27].

2 FORMULATION

The wave function for H₂⁺ is expanded in terms of basis functions φ_nlm in polar coordinates, as   

$\Psi (r,\theta ,\varphi )=\sum _{n=1}^{\infty }\sum _{l=0}^{n-1}\sum _{m=-l}^l{C}_{nlm}{\phi }_{nlm}(r,\theta ,\varphi ),$(1)

  

${\phi }_{nlm}(r,\theta ,\varphi )=R_{nl}\left(r\right)Y_{lm}(\theta ,\varphi )$(2)

Here, Y_lm is a spherical harmonic function. The expansion coefficients C_nlm are obtained by the variational principle,   

$\delta \int {\Psi }^{*}H\Psi dv=0,$(3)

  

$H=-\frac{1}{2}{\nabla }^2-\frac{Z}{\left|r-R_A\right|}-\frac{Z}{\left|r-R_B\right|}+\frac{Z^2}{\left|R_A-R_B\right|}$(4)

Here H is the Hamiltonian, and R_A and R_B are the position vectors of nuclei A and B. Z is the nuclear charge. For H₂⁺, Z is 1. The fourth term, the nuclear repulsion energy, is constant in the Born-Oppenheimer approximation and will be omitted henceforth. As a result we consider the electronic energy (EE) rather than the total energy (TE). Atomic units are used throughout.

We adopt a single-center expansion based on the molecular midpoint for the wave function. Laguerre-type orbitals (LTOs) [1] are used for the radial functions,   

$R_{nl}\left(r\right)={\left(2z_l\right)}^{3/2}\frac{\sqrt{(n-l-1)!}}{\sqrt{(n+l+1)!}}{\left(2z_{l}r\right)}^l{L}_{n-l-1}^{2l+2}\left(2z_{l}r\right)e^{-z_{l}r}$(5)

Here $L_{n-l-1}^{2l+2}\left(2z_{l}r\right)$ is an associated Laguerre polynomial; z_l is called the exponent. The polynomial is written as   

$L_n^m\left(x\right)=\sum _{k=0}^n\frac{(n+m)!{(-x)}^k}{(n-k)!(m+k)!k!}.$(6)

These polynomials satisfy the following recurrence relation:   

$\begin{array}{c}n{L}_n^m\left(x\right)=(2n+m-1-x)L_{n-1}^m\left(x\right)\\ -(n+m-1)L_{n-2}^m\left(x\right),\end{array}$(7)

  

$L_0^m\left(x\right)=1,L_1^m\left(x\right)=m+1-x.$(8)

The Hamiltonian element consists of three terms,   

$\begin{array}{c}\int {\phi }^{*}{}_{nlm}H{\phi }_{n^{\prime }{l}^{\prime }{m}^{\prime }}dv=T(nlm,n^{\prime }{l}^{\prime }{m}^{\prime })\\ +V_A(nlm,n^{\prime }{l}^{\prime }{m}^{\prime })+V_B(nlm,n^{\prime }{l}^{\prime }{m}^{\prime })\end{array}$(9)

The first term, representing the kinetic energy, is evaluated according to the following formula, given in an Appendix (Eq.7) of the paper by Hagstrom and Shull [11],   

$\begin{array}{c}T(nlm,n^{\prime }{l}^{\prime }{m}^{\prime })=\left(nlm\right|-\frac{1}{2}{\nabla }^2\left|n^{\prime }{l}^{\prime }{m}^{\prime }\right)\\ =z_l^2{\delta }_{l{l}^{\prime }}{\delta }_{m{m}^{\prime }}n\times \left\{\frac{2n+1}{2l+3}-\frac{{\delta }_{n{n}^{\prime }}}{2}\right\}{\left\{\frac{(n+l+1)!(n^{\prime }-l+1)!}{(n-l-1)!(n^{\prime }+l+1)!}\right\}}^{\frac{1}{2}}\\ (n\le n^{\prime })\end{array}$(10)

The second term in the Hamiltonian,   

$V_A(nlm,n^{\prime }{l}^{\prime }{m}^{\prime })=Z\int {\phi }^{*}{}_{nlm}\left(-\frac{1}{\left|r-R_A\right|}\right){\phi }_{n^{\prime }{l}^{\prime }{m}^{\prime }}dv,$(11)

represents the electron-nucleus Coulomb energy. The Coulomb term is expanded in terms of Legendre polynomials as   

$\frac{1}{\left|r-R_A\right|}=\sum _{s=0}^{\infty }\frac{r_{<}^s}{r_{>}^{s+1}}{P}_s\left(\cos {\theta }_A\right),$(12)

where r_< = min (r, R_A), r_< = max (r, R_A). The second and third terms in the Hamiltonian are the most difficult to calculate.

The Hamiltonian is nonzero only for m = 0. For H₂⁺ (1s σ_g), both l and l’ are even. By substituting R and V for R_A and V_A, V can be written as follows:   

$\begin{array}{c}V(nl,n^{\prime }{l}^{\prime };R)=Z\sum _{s=0}^{\min (l,l^{\prime })}\frac{\sqrt{(2l+1)(2l^{\prime }+1)}}{2l+2l^{\prime }-2s+1}\\ \times \frac{A_{l-s}{A}_s{A}_{l^{\prime }-s}}{A_{l+l^{\prime }-s}}\times K(n,l,n^{\prime },l^{\prime };s,R)\end{array}$(13)

  

$A_m=\frac{(2m-1)!!}{m!}$(14)

  

$K(n,l,n\text{'},l\text{'};s,R)\text{ =}\underset{0}{\overset{\infty }{\int }}{R}_{nl}\left(r\right)R_{n^{\prime }{l}^{\prime }}\left(r\right)\frac{r_{<}{}^{l+l^{\prime }-2s}}{r_{>}{}^{l+l^{\prime }-2s+1}}{r}^{2}dr$(15)

The integral K consists of two parts: F of finite interval [0, R], and G of infinite interval [R, ∞), as follows:   

$\begin{array}{c}F(n,l,n\text{'},l\text{'};s,R)=\\ \underset{0}{\overset{R}{\int }}{r}^{2(l+l^{\prime }-s+1)}{L}_{n-l-1}^{2l+2}\left(2z_{l}r\right)L_{n^{\prime }-l^{\prime }-1}^{2l^{\prime }+2}\left(2z_{l^{\prime }}r\right)e^{-(z_l+z_{l^{\prime }})r}dr\end{array}$(16)

  

$\begin{array}{c}G(n,l,n\text{'},l\text{'};s,R)=\\ \underset{R}{\overset{\infty }{\int }}{r}^{2s+1}{L}_{n-l-1}^{2l+2}\left(2z_{l}r\right)L_{n^{\prime }-l^{\prime }-1}^{2l^{\prime }+2}\left(2z_{l^{\prime }}r\right)e^{-(z_l+z_{l^{\prime }})r}dr\end{array}$(17)

Then K becomes   

$\begin{array}{c}K(n,l,n\text{'},l\text{'};s,R)=\\ \left(\frac{F(n,l,n^{\prime },l^{\prime };s,R)}{R^{l+l^{\prime }-2s+1}}+R^{l+l^{\prime }-2s}G(n,l,n^{\prime },l^{\prime };s,R)\right)\\ \times {\left(}^2\sqrt{\frac{(n-l-1)!}{(n+l+1)!}}{\left(}^2\sqrt{\frac{(n^{\prime }-l^{\prime }-1)!}{(n^{\prime }+l^{\prime }+1)!}}\end{array}$(18)

3 NUMERICAL INTEGRATION

3.1 Term-by-term analytical integration

By inserting Eq.6 into Eq.16 and Eq.17, the integrals F and G can be written as the sum of the individual terms:   

$\begin{array}{c}F(n,l,n\text{'},l\text{'};s,R)\text{= }\\ (n+l+1)!(n\text{'}+l\text{'}+1)!\\ \times \sum _{k=0}^{n-l-1}\frac{{(-2z_l)}^k}{(n-l-1-k)!(2l+2+k)!k!}\\ \times \sum _{k^{\prime }=0}^{n^{\prime }-l^{\prime }-1}\frac{{(-2z_{l^{\prime }})}^{k^{\prime }}}{(n^{\prime }-l^{\prime }-1-k^{\prime })!(2l^{\prime }+2+k^{\prime })!k^{\prime }!}\\ \times {\int }_0^R{r}^{2(l+l^{\prime }-s)+k+k^{\prime }+2}{e}^{-(z_l+z_{l^{\prime }})r}dr\end{array}$(19)

  

$\begin{array}{c}G(n,l,n\text{'},l\text{'};s,R)\text{= }\\ (n+l+1)!(n\text{'}+l\text{'}+1)!\\ \times \sum _{k=0}^{n-l-1}\frac{{(-2z_l)}^k}{(n-l-1-k)!(2l+2+k)!k!}\\ \times \sum _{k^{\prime }=0}^{n^{\prime }-l^{\prime }-1}\frac{{(-2z_{l^{\prime }})}^{k^{\prime }}}{(n^{\prime }-l^{\prime }-1-k^{\prime })!(2l^{\prime }+2+k^{\prime })!k^{\prime }!}\\ \times {\int }_R^{\infty }{r}^{k+k^{\prime }+2s+1}{e}^{-(z_l+z_{l^{\prime }})r}dr\end{array}$(20)

Analytical integration is possible for the individual terms [28] in this expression. However, F and G both contain a series of terms which alternate in sign as k(or k’) is incremented. This causes severe cancellations of significant digits.

We examined the cancellations for the case of n = n’ = 151, l = l’ = 0, s = 0, z_l = z_l’ = 37.65. The partial sum of F grew to 1.2357857664 × 10¹¹⁶ at the maximum, and F(151,0,151,0; 0,1) finally converged to 9.3719206312 × 10⁻². It is immediately apparent that cancellation of 118 decimal digits was involved. Furthermore, the partial sum of G grew to 3.2394517384 × 10¹³⁸ although G(151,0,151,0; 0,1) finished up as 1.1259694138 × 10⁻². In this case the cancellation amounts to 140 digits.

These facts make clear that the quadruple precision arithmetic supported by current Fortran systems, ensuring 33 significant figures, is not sufficient. Accordingly a multiple-precision arithmetic package (such as MPfun90 [19], FM [20]) is required. In fact, 166 decimal digits were needed for the examination above. Such multiple-precision arithmetic requires much CPU time, however. Moreover the computational cost of evaluating one integral of V is of the order of n³. For these reasons we abandoned term-by-term analytical integration for large-scale calculations, using it only for testing purposes. Instead, we used selected numerical integration methods. Details are set out in the next subsection.

3.2 Gaussian quadratures

For numerical integration we adopt the Gaussian quadrature formula [21,22,23].   

${\int }_a^{b}h\left(x\right)dx={\int }_a^{b}w\left(x\right)f\left(x\right)dx=\sum _{k=1}^M{W}_{k}f\left(x_k\right)+\text{error}$(21)

Here w(x) is a weighting function, x_k are abscissas, W_k are weights, and M is the number of abscissas.

The integrands of Eq.16 and Eq.17 are plotted in Figure 1 and Figure 2, respectively. These figures reveal that the integrands oscillate frequently for large n or n’, and that as many abscissas as the frequency should be taken in order to provide sufficient accuracy for evaluation of the integrals.

Figure 1.

 The integrand of F(151,0,151,0; s = 0, R = 1) over [0,1].

Figure 2.

 The integrand of G(151,0,151,0; s = 0, R = 1) over [1,10].

Since the error in Gaussian quadrature can be made zero for G, but not for F, we describe the numerical integration of G(Eq.17) first, to simplify the discussion. Later we shall consider F(Eq.16).

Equation 17 can be changed to the following form by transformation of variable:   

${\int }_0^{\infty }g\left(x\right)e^{-x}dx,$(22)

to which we apply the Gauss-Laguerre quadrature rule [21] on (a,b) = (0,∞) with the weighting function w(x) = e^−x.   

$G(n,l,n^{\prime },l^{\prime };s,R)=\sum _{k=1}^M{W}_{k}g\left(x_k\right)+\text{error}$(23)

  

$\begin{array}{c}g\left(x\right)=\frac{1}{z_l+z_{l^{\prime }}}{e}^{-(z_l+z_{l^{\prime }})R}{\left(\frac{x}{z_l+z_{l^{\prime }}}+R\right)}^{2s+1}\\ \times L_{n-l-1}^{2l+2}\left(2z_l\left(\frac{x}{z_l+z_{l^{\prime }}}+R\right)\right)\\ \times L_{n^{\prime }-l^{\prime }-1}^{2l^{\prime }+2}\left(2z_{l^{\prime }}\left(\frac{x}{z_l+z_{l^{\prime }}}+R\right)\right)\end{array}$(24)

  

$J=2s+n-l+n\text{'}-l\text{'}-1.$(25)

Here J is the order of the polynomial. To make the error zero, J should be ≤ 2M − 1. Then M takes the form,   

$M=\lceil (2s+n-l+n^{\prime }-l^{\prime })/2n\rceil$(26)

We next checked the accuracy of the integrals. The results of the test are listed below to 33 decimal places, taking the case of n = n’ = 151, l = l’ = 0, s = 0, z_l = z_l’ = 37.65 as an example. For this check we wrote a test program by modifying subroutines in the mathematical library NUMPAC [29,30] (Y.H. is one of its authors). The numerical integration was performed with quadruple precision arithmetic. Results from the Gaussian quadrature are compared with those obtained by analytical integration using the multiple-precision arithmetic package MPfun90 [19]. The result from the analytical integration with 158 significant digits specified for MPfun90 is quoted below as scheme (A). In the scheme (B)/(C), 166/173 digits are specified.    G for n = n’ = 151, l = l’ = 0, s = 0, z_l = z_l’ = 37.65.   

(a) Gauss 151 points

1.125969413770606624522211291992090E-02

(A) MPfun90 158 digits

1.125969413770606624522211268404800E-02

(B) MPfun90 166 digits

1.125969413770606624522211291991500E-02

(C) MPfun90 173 digits

1.125969413770606624522211291991500E-02

This test shows that M defined by Eq.26 works well. Since the order of the polynomial (J) reaches its maximum when l = l’ = s, the maximum M is given by Nmax. Thus, for Nmax = 203, M is 203.

Using the Golub and Welsch algorithm [22], the abscissas (zeros of the Laguerre polynomials L_M(x)) are calculated efficiently from the eigenvalues of a symmetric tridiagonal matrix that is constructed from the three-term recurrence relation (Eq.7). The weights [21] are:   

$W_k=\frac{{(M!)}^2{x}_k}{{\left[L_{M+1}\left(x_k\right)\right]}^2}.$(27)

We turn now to the evaluation of F(Eq.16). The Gauss-Legendre quadrature rule [23] is used, with the weighting function w(x) = 1.   

$F(n,l,n^{\prime },l^{\prime };s,R)=\sum _{k=1}^M{W}_{k}f\left(x_k\right)+\text{error}$(28)

  

$f\left(x\right)=x^{2(l+l^{\prime }-s+1)}{L}_{n-l-1}^{2l+2}\left(2z_{l}x\right)L_{n^{\prime }-l^{\prime }-1}^{2l^{\prime }+2}\left(2z_{l^{\prime }}x\right)e^{-(z_l+z_{l^{\prime }})x}$(29)

  

$J=l\text{ + }l\text{'}-2s+n+n\text{'}$(30)

Since the integrand involves an exponential function, the error cannot be zero in this case. It is not easy to estimate the number of abscissas necessary by theory alone. We therefore performed trial calculations to find the number of abscissas that best balances accuracy and efficiency. We first tested a value of M obtained from J ≤ 2M − 1, namely   

$M=\lceil (l+l^{\prime }-2s+n+n^{\prime }+1)/2\rceil$(31)

We tested also the double of this value,   

$M=l+l\text{'}-2s+n+n\text{'}+1$(32)

We list below one of the results of the trial calculations.    F for n = n’ = 125, l = l’ = 20, s = 5, z_l = z_l’ = 114.98528.   

(a’) Gauss 141 points

2.878459649515099324510714516949794E-16

(b’) Gauss 281 points

2.878459649515099317492932366032297E-16

(A’) MPfun90 158 digits

2.878459649515099317492932366230900E-16

(B’) MPfun90 166 digits

2.878459649515099317492932366230900E-16

(C’) MPfun90 173 digits

2.878459649515099317492932366230900E-16

The scheme (a’) gives only 17 significant digits, whereas the scheme (b’) gives 28 digits. To obtain the desired accuracy, scheme (b’) is needed. We therefore adopt scheme (b’), i.e., we set M to J+1. Since J reaches its maximum when l = l’ = Lmax and s = 0, the maximum M is given by 2 × Nmax + 2 × Lmax +1. For Nmax = 203 and Lmax = 202, M takes the value 811. Here Lmax is the maximum of l or l’ in the basis.

The abscissas (zeros of the Legendre polynomials P_m(x)) are calculated using Newton’s iterative method, starting from approximate roots given by   

$x_k=\sin \left(\frac{M+1-2k}{2M+1}\right)\pi (k=1,2,\dots ,M).$(33)

The weights [23], W_k, are calculated from the following formula:   

$W_k=\frac{2(1-x_k{}^2)}{{\left[M{P}_{M-1}\left(x_k\right)\right]}^2}$(34)

To conclude this subsection, we compare the Coulomb integral (V) obtained from the Gaussian quadrature scheme (a and b’) with that evaluated from Eq.13 using MPfun90 [19].

V for n = n’ = 151, l = l’ = 0, z_l = z_l’ = 37.65.

Gaussian quadratures

0.439338507119362656787621563083154

MPfun90 158 digits

0.439338507119362656787621559152698

MPfun90 166 digits

0.439338507119362656787621563540349

MPfun90 173 digits

0.439338507119362656787621563540349

V for n = n’ = 151, l = l’ = 150, z_l = z_l’ = 37.65.

Gaussian quadratures

0.253398692737192741982923026971287

MPfun90 158 digits

0.253398692737192741982923026971291

MPfun90 166 digits

0.253398692737192741982923026971291

MPfun90 173 digits

0.253398692737192741982923026971291

Further tests confirm at least 20-figure accuracy for Nmax ≤ 203. The cost of the numerical integration is O(Nmax⁶).

4 PROGRAM STRUCTURE

4.1 Calling sequence

The program system provided is self-contained, so that other library routines are not necessary. The calling sequence of the major routines is shown in Figure 3. For Gauss-Legendre quadrature, subroutine TGLEGQ of NUMPAC library [29,30] is incorporated. For diagonalizing the Hamiltonian matrix to obtain EE, subroutine HOBSVQ of the NUMPAC library [29,30] is used.

Figure 3.

 Calling sequence.

4.2 Array size

The memory size required to run the program is determined by one large array named HMATRX, where the Hamiltonian elements are kept. The size of HMATRX is NNQH × NNQH, while NNQH = Nmax × (Nmax+1)/4 + Nmax. With Nmax = 203, the size of HMATRX exceeds 2 GB assuming quadruple precision of 16 B.

5 PROGRAM MANUAL

5.1 Installation

First, expand the downloaded file. Here we take its file name as scelto_gauss.tar.gz [31].

tar -zxvf scelto_gauss.tar.gz

Compile the source programs in the directory src using the make command.

cd src

make

cd ..

An executable program named sscelto_gauss.exe will be generated. Move it into the directory bin.

mv src/scelto_gauss.exe bin/

Run the program.

cd wshop.N19

../bin/scelto_gauss.exe

The Intel Fortran compiler is the default. For the Fujitsu compiler, Makefile.Fujitsu is available. To use other compilers, the command name and compile options must be changed. In particular, some compile options may have to be changed; the large array HMATRX) should therefore be placed in a memory region beyond 2 GB.

The value of the parameter constant NQMAX, which defines the upper limit for Nmax, is set to 211. By compiling the source programs with smaller NQMAX, the program size can be reduced.

Test input data c5nqg.inp is stored in the directory wshop.N19. This data is for the case Nmax = 19. It gives EE of− 1.10227888655218872575 a.u.

5.2 Input data

We describe the input data parameters, taking the input file in the directory wshop.N19 as an example.

Line: Variable or array / Meaning

1: comment

2: Nmax / maximum of n.

3: comment

4: Lmax / maximum of l.

5: comment

6: R / nuclear position (a.u.) from the midpoint.

7: comment

8: ZNUC / nuclear charge. (1.0 for H₂⁺)

9: comment

10–19: ZL / exponents for even l numbers.

5.3 Output files

The following files will be generated after execution of theprogram.

c5nqg.outlist

c5nqg.log

c5nqg.out9.txt

c5nqg.out.txt

For each value of l, the Hamiltonian matrix is diagonalized and the exponents, eigenvalues, and eigenvectors are written to the file c5nqg.out.txt. Immediately after the diagonalization for a particular l is finished, the number l is output to the file c5nqg.log, so that users can observe the progress. To the file c5nqg.out9.txt, the eigenvalues alone are output.

6 PERFORMANCE

Table 1 shows the computational results. The fourth column, named Ndim, gives the dimensions of the Hamiltonian matrices. For the method titled ‘Scaling’, identical exponents were used for all the l numbers. For Nmax = 203, the scaling method yields the energy as −1.10263287 a.u, which differs from the exact energy by 1.34 × 10⁻⁶ a.u. The exponents were determined by extrapolation from the preceding runs of n = 1–70.    z_l = 1.6193375636179528 + 0.2782894519200024n – 0.0005804140783906168n² + 0.000002103469132713343n³    (35)   

Table 1.  Electronic energy of hydrogen molecular ion at R_bond=2.0 a.u.

Method	Nmax	Lmax	Ndim	Exponent(zl)	EE(a.u.)	CPU-time(sec)*
Scaling	51	50	676	14.59721	−1.10258505	271
Scaling	71	70	1296	19.26503	−1.10261435	2235
Scaling	101	100	2601	25.97300	−1.10262669	22493
Scaling	131	130	4356	32.84360	−1.10263056	126906
Scaling	203	202	10404	51.18760	−1.10263287	2521290
Full opt.†	131	130	4356		−1.10263323
Emp. opt.‡	203	202	10404		−1.10263395
Exact (by Peek [14])					−1.10263421

* CPU-time on an Intel Core i7-3970X (3.5 GHz) in serial mode.

^† The exponents are fully optimized.

^‡ The exponents are optimized empirically. See text.

By optimizing the exponents for l individually, EE has been decreased further to −1.10263395 a.u. The difference between the final energy for Nmax = 203 and the exact value is less than 2.6 × 10⁻⁷ a.u.

For Nmax = 203, the exponents were optimized in an empirical manner. First, the exponent for l = 0 (denoted as z₀) is optimized, with Lmax set to 0. Next, the exponent for l = 2 (z₂) is optimized from an initial value given by multiplying z₀ by 1.1, with Lmax set to 2. By repeating this procedure the exponents are determined, until l = 50. The exponent for l = 202 (z₂₀₂) is fixed to 203.5 by n+0.5. The exponents for l = 52–200 are calculated from a line connecting z₅₀ and z₂₀₂ (The optimization routines are not provided.)

The electronic energies are plotted versus Nmax in Figure 4. The blue curve marked with squares shows the energy obtained by the scaling method. The red curve marked with circles shows the energy obtained with optimized exponents. Solid circles indicate that the exponents are fully optimized. The open circles for Nmax = 161–203 signify that the exponents are optimized empirically. The calculated energies approach the exact value.

Figure 4.

 Electronic energy vs Nmax.

NOTE ADDED IN PROOF

To observe the l-dependency of the eigenvalue, the Hamiltonian was diagonalized for each l-block (l = 0,1,...,Lmax). The CPU-time in Table 1 includes all the diagonalization steps for l. If the Hamiltonian is diagonalized only once for the maximum l(Lmax), the CPU time is drastically reduced. Moreover, we have found that the Hamiltonian can be diagonalized in the double-precision without loss of accuracy. By considering these 2 factors, we have reduced the CPU-time for Nmax = 203 to 681190 sec (about 1/4 of the original). Finally we succeeded to extend n to 241 and obtained EE = −1.10263406 a.u.

For the case of Nmax = 39 and Lmax = 38, EE was calculated to be −1.10253135 a.u. We have performed additional calculations by locating the expansion center at one nucleus and have obtained EE = −1.10244106 a.u. We call this computational scheme one-side expansion against midpoint expansion described in text. The midpoint expansion gives slightly better EE. In the one-side expansion, we have also calculated the nuclear cusp condition to be −2.34. Since in the midpoint expansion, the peak position differs from the nuclear position, the calculated cusp condition loses its exact physical meaning. However, we estimated its value to be −0.33 by numerical differentiation along the Z axis. Along the y axis (perpendicular to the molecular axis), it was − 1.0 × 10⁻⁵. In all these cases, the calculated values deviate considerably from the true value (−1). Although the LTOs themselves are considered to be superior to the GTOs for computing the cusp conditions, their single-center expansion gives unsatisfactory results.

Acknowledgements

We are grateful to Professor Emeritus Hiroshi Tatewaki of Nagoya City University for valuable suggestions. Y.H. thanks Professor Hiroshi Sugiura at Nanzan University and Professor Emeritus Takemitsu Hasegawa at Fukui University for their helpful suggestions regarding Gaussian quadratures. The main computation was performed using a cluster of Intel Xeon-based computers located at the Institute for Advanced Studies in Artificial Intelligence (IASAI) of Chukyo University. The Fujitsu PRIMEHPC FX100 supercomputer at the Information Technology Center of Nagoya University was used for portability testing.

REFERENCES

• [1]   T. Helgaker, P. Jørgensen, J. Olsen, Molecular Electronic-Structure Theory, John Wiley & Sons, Chichester, 2000.
• [2]   M. W. J. Bromley, J. Mitroy, Int. J. Quantum Chem., 107, 1150 (2007).
• [3]   L. G. Jiao, Y. K. Ho, Int. J. Quantum Chem., 115, 434 (2015).
• [4]   M. T. Marron, N. C. Handy, R. G. Parr, H. J. Silverstone, Int. J. Quantum Chem., 4, 245 (1970).
• [5]   J. Hachmann, P. T. A. Galek, Y. Yanai, G. K.-L. Chan, N. C. Handy, Chem. Phys. Lett., 392, 55 (2004).
• [6]   N. Takeda, Y. Hatano, S. Yamamoto, H. Tatewaki, Comput. Theor. Chem., 1045, 99 (2014).
• [7]   R. Courant, D. Hilbert, Methods of Mathematical Physics, Vol. I, Interscience, New York, 1953.
• [8]   G. Szegö, Orthogonal Polynomials, American Mathematical Society, Providence, Rhode Island, 1939.
• [9]   K. M. Howell, H. Shull, J. Chem. Phys., 30, 627 (1959).
• [10]   H. H. Kranz, E. O. Steinborn, Phys. Rev. A, 25, 66 (1982).
• [11]   S. Hagstrom, H. Shull, J. Chem. Phys., 30, 1314 (1959).
• [12]   Ø. Burrau, Det. Kgl. Danske Vid. Selskab., 7, 1 (1927).
• [13]   E. A. Hylleraas, Z. Phys., 71, 739 (1931).
• [14]   J. M. Peek, J. Chem. Phys., 43, 3004 (1965).
• [15]   R. M. Fernández, J. Chem. Phys., 103, 6581 (1995).
• [16]   D. R. Bates, K. Ledsham, A. L. Stewart, Phil. Trans. Roy. Soc. London A, 246, 215 (1953).
• [17]   T. Kato, Commun. Pure Appl. Math., 10, 151 (1957).
• [18]   B. H. Wells, S. Wilson, J. Phys. At. Mol. Opt. Phys., 22, 1285 (1989).
• [19]   D. H. Bailey, ACM Trans. Math. Softw., 21, 379 (1995).
• [20]   D. M. Smith, J. Comput. Sci. Eng., 5, 88 (2003).
• [21]   F. B. Hildebrand, Introduction to Numerical Analysis, 2nd ed., Dover, New York, 1987.
• [22]   G. H. Golub, J. H. Welsch, Math. Comput., 23, 221 (1969).
• [23]   P. Davis, P. Rabinowitz, J. Res. Natl. Bur. Stand., 56, 35 (1956).
• [24]   IEEE Computer Society, (August 29, 2008), IEEE Standard for Floating-Point Arithmetic, IEEE,doi: .
• [25]   J. A. Kempe, S. P. Goldman, J. Chem. Phys., 108, 7679 (1998).
• [26]   A. Ishikawa, H. Nakashima, H. Nakatsuji, J. Chem. Phys., 128, 124103 (2008).[CrossRef]
• [27]   Y. Ning, Z.-C. Yan, Phys. Rev. A, 90, 032516 (2014).[CrossRef]
• [28]   K. Matsubara, Y. Horino, T. Murai, Bull. Nara Univ. Educ., 17, 13 (1969).
• [29]   I. Ninomiya, Y. Hatano, “Performance of the Vector Version of NUMPAC”, in: High Performance Computing: Research and Practice in Japan, R. Mendez (Ed.), John Wiley & Sons, Chichester, 1992. The source programs and the manual can be downloaded from http://www2.itc.nagoya-u.ac.jp/center/ja/numpac/ index.html.
• [30]   T. Hasegawa, T. Sandoh, Y. Sato, Y. Hatano, I. Ninomiya, Memoirs of the Faculty of Engineering, Fukui Univ., 48, 193 (2000).
• [31]   scelto_gauss.tar.gz