BornOppenheimer potential for H
Abstract
The BornOppenheimer potential for the state of H is obtained in the range of 0.1 – 20 au, using analytic formulas and recursion relations for twocenter twoelectron integrals with exponential functions. For small distances JamesCoolidge basis is used, while for large distances the HeitlerLondon functions with arbitrary polynomial in electron variables. In the whole range of internuclear distance about precision is achieved, as an example at the equilibrium distance au the BornOppenheimer potential amounts to . Results for the exchange energy verify the formula of Herring and Flicker [Phys. Rev. 134, A362 (1964)] for the large internuclear distance asymptotics. The presented analytic approach to Slater integrals opens a window for the high precision calculations in an arbitrary diatomic molecule.
pacs:
31.15.ac, 31.50.BcI Introduction
In order to accurately obtain vibrational and rotational spectra of molecules not only the nonrelativistic BornOppenheimer (BO) potential has to be calculated, but also the adiabatic, nonadiabatic, relativistic and quantum electrodynamics effects. An excellent agreement of theoretical results h2 () with experimental values achieved for dissociation energies of the hydrogen molecule H h2_exp () and D d2_exp () indicates good understanding of all physical effects up to precision level. Such a precision has not yet been achieved for molecular systems with more than two electrons. It is because numerical calculations with explicitly correlated functions are very demanding. The commonly used explicitly correlated Gaussian functions require global minimization of thousands of nonlinear parameters, and the accuracy achieved for relativistic effects is difficult to estimate, due to the improper analytic properties of these functions. The largest system considered so far is the helium dimer, where about precision is achieved for the interaction energy he2 ().
Herewith, we apply explicitly correlated exponential functions (for a recent review see Ref. correlated ()), in order to increase numerical accuracy obtained so far for two electrons diatomic systems. The use of exponential functions may allow for a significant improvement not only in rotational and vibrational energies, but also in other important properties like magnetic shielding, spinrotational or spinspin coupling constants. It is because evaluation of relativistic effects is very sensitive to correct analytic properties of the basis functions, a good example being the relativistic correction to the magnetic shielding rudzins ().
When precision is achieved for vibrational transitions, the electronproton mass ratio can be determined with better accuracy, than it is known presently nist (). Moreover, any nonelectromagnetic longrange interactions between nuclei, which cannot in principle be excluded, may be visible not only in vibrational spectra, but also in the spin interaction between nuclei spindep (), for example in HD molecule. The electromagnetic spin interaction is extremely small, of order , about 43 Hz for the ground molecular state chemrev (), and any deviation between experimental values and theoretical predictions will signal existence of a nonelectromagnetic longrange interactions between nuclei.
In this work we present an effective computational method with the use of exponential functions and numerical results for the BornOppenheimer energy of the H molecule. We demonstrate that numerical precision can be achieved for BO energy with about 22 000 basis functions using the analytic formulae and recursion relations for Slater integrals, which have been obtained in our previous work in rec_h2 (). The obtained BO potential is at least two orders of magnitude more accurate, than anything published so far, and at the equilibrium distance our calculations confirm the most accurate so far result of Cencek and Szalewicz in Ref. cencek () well within their uncertainties. Our results for large internuclear distances verify the asymptotic formula of Herring and Flicker herring (), although we observe significant contributions from the off leading terms. Finally, we analyze possibility of the extension of this method to more than twoelectron systems.
Ii Short range of internuclear distances
The Schrödinger equation for H molecule in the BornOppenheimer approximation is
(1) 
where
(2) 
The most efficient basis set to solve this equation is the one introduced by Kołos and Wolniewicz in kw_basis (). However at small nuclear distances , we can use much simpler JamesCoolidge basis james_cool (), namely the functions of the form
(3)  
such that
(4) 
with . This basis have recently been used by Sims and Hagstrom SH06 () for precision calculation of BO potential in the range au, by a traditional way of numerical evaluation of the corresponding twocenter integrals. Here, using Ref. rec_h2 () we derive analytic expression for all these integrals up to . They have very simple form and involve only exponential integral (Ei), exponential and logarithmic functions. In this derivation we make use of a differential equation, which is satisfied by the most general twocenter twoelectron Slater integral, see Ref. rec_h2 ().
For the presentation of analytic formulas we consider an integral
(5)  
at the internuclear distance , since arbitrary distances can be obtained from the relation
(6) 
Positive powers of and can be obtained by differentiation with respect to and correspondingly, so without loss of generality one can assume . Let us introduce auxiliary functions defined by
(7)  
(8)  
(9)  
(10)  
(11) 
Then, using the equation (69) from rec_h2 (), all nonvanishing integrals with are the following:
(12)  
(13)  
(14)  
(15)  
(16)  
(17)  
(18) 
(19)  
(20)  
(22)  
(23)  
(24)  
Integrals with higher powers of are similar, with simple polynomial form in and . We have tabulated all integrals, such that what is sufficient for matrix elements involving exponential functions with maximum value of .
Iii Long range of internuclear distances
For the large internuclear distance the JamesCoolidge basis is not appropriate as it does not include the HeitlerLondon function. Instead, we employ generalized HeitlerLondon functions, which are the product of a Heitler London function with an arbitrary polynomial in all electron distances,
(25) 
such that
(26) 
and call them the explicitly correlated asymptotic (ECA) basis. Matrix elements of the nonrelativistic Hamiltonian can be expressed in terms of direct integrals of the form
(27) 
with nonnegative integers , and exchange integrals which coincide with that from the previous section. The direct integrals are calculated as follows. When all the so called master integral is given by rec_h2 ()
and its first derivative with respect to
(29) 
where
(30) 
This one dimensional integration can easily be obtained with the sufficient accuracy and for example at the value of is
(31) 
All integrals with higher powers of electron distances can be obtained from recursion relations, which were derived in Ref. rec_h2 (). Since
(32) 
we can present formulae at . Let us introduce auxiliary functions
(33)  
(34)  
(35)  
(36)  
(37) 
then all integrals with are the following
(38)  
(39)  
(40)  
(41)  
(42)  
(43)  
(44)  
(45)  
(46)  
(47)  
(48)  
(49)  
(50)  
(51) 
Integrals with higher powers are of analogous form, they are all linear combinations of , , and with coefficients being polynomials of . We have generated a table of integrals with , which corresponds to a maximum value of . It is less than in the case of JamesCoolidge functions, but anyway, requires the use of octuple precision arithmetics due to near linear dependence of these basis functions.
Iv Numerical results
The matrix elements of the nonrelativistic Hamiltonian between exponential functions are obtained as described by Kołos and Roothan in kolroth (). The resulting expression is a linear combination of various Slater integrals which are calculated using analytic formulas as presented in the previous sections. This evaluation is fast and accurate, thus allowing the use of a large number of basis functions.
Eigenvalues of the Hamiltonian matrix are obtained by inverse iteration method for various length of the basis set. Following Sims and Hagstrom SH06 () we use double basis set with two different nonlinear parameters, which were obtained by minimization of energy at . For the calculations with basis functions up to , the second nonlinear parameter is additionally multiplied by 1.5 in order to improve the numerical stability for the large values of . Numerical calculations with JamesCoolidge functions are performed in general using quadruple precision arithmetics, and for checking the numerical accuracy, the one point at au is calculated using the octuple precision. We observe a significant loss of digits for large values of . The quadruple precision arithmetics was not always sufficient for the whole range of internuclear distances. In many cases values with , and sometimes even had to be disregarded due to numerical instabilities in the inverse iteration procedure. Numerical results at and au for various sizes of basis length are presented in Table IV.
N  

3  23  
4  51  
5  98  
6  180  
7  306  
8  501  
9  781  
10  1182  
11  1729  
12  2471  
13  3444  
14  4712  
15  6324  
16  8361  
17  10887  
18  14002  
19  17787  
20  22363  
Cencek cencek () 
The most accurate result obtained at the equilibrium distance au is compared in Table IV to all the previous results obtained so far in the literature.
Authors  type of w.f.  N  

1933 James and Coolidge james_cool ()  JC  5  
1960 Kołos and Roothan kolroth ()  JC  50  
1968 Kołos and Wolniewicz kolwol ()  KW  100  
1994 Rychlewski, Cencek, Komasa rych ()  ECG  700  
1995 Wolniewicz wol ()  KW  883  
2006 Sims and Hagstrom SH06 ()  JC  7034  
2007 Nakatsuji et al.nakatsuji ()  ICI  6776  
2008 Cencek, Szalewicz cencek ()  ECG  4800  
2010 this work  JC  22363 
In performing extrapolation, we observe the exponential convergence. In other words the of differences in energies for subsequent values of fits well to the linear function. This allows for a simple and reliable extrapolation to infinity. For example at the equilibrium distance the parameter is about , and convergence in the JamesCoolidge basis preserve its exponential behavior for all the distances.
The calculations for the internuclear distances au. are performed in ECA basis set using octuple precision arithmetics up to . This basis is probably the most effective one at large distances. We also observe exponential convergence for BO and even exchange energies, what makes quite simple the extrapolation to infinity, see Table IV with detailed results for various length of the basis set.

It was found by Herring and Flicker in herring () that the asymptotic limit is with . The fit to our numerical data, see Fig. 1, gives what can be regarded as a first numerical confirmation of the Herring and Flicker result. When assuming their constant, the fit of first three terms in expansion gives: , what indicates very slow convergence of this expansion at typical interatomic distances . V SummaryWe have demonstrated the applicability of analytic formulas for twocenter twoelectron integrals in high precision calculations of BornOppenheimer potential for the ground electronic state of the hydrogen molecule. The symmetric JamesCoolidge basis with as much as 22 000 functions provided energies with relative precision of about for internuclear distances up to 10 au. The extended HeitlerLondon basis with about 10 000 functions give similarly accurate description at internuclear distances au and greater. The extension of the analytic approach to excited states of H and other diatomic molecules, such as HeH requires the evaluation of Slater integrals with arbitrary nonlinear parameters (KołosWolniewicz basis kw_basis ()). 