Born-Oppenheimer potential for H{}_{2}

Born-Oppenheimer potential for H

Krzysztof Pachucki Institute of Theoretical Physics, University of Warsaw, Hoża 69, 00-681 Warsaw, Poland

The Born-Oppenheimer potential for the state of H is obtained in the range of 0.1 – 20 au, using analytic formulas and recursion relations for two-center two-electron integrals with exponential functions. For small distances James-Coolidge basis is used, while for large distances the Heitler-London 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 Born-Oppenheimer 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.50.Bc
preprint: Version 3.0

I Introduction

In order to accurately obtain vibrational and rotational spectra of molecules not only the nonrelativistic Born-Oppenheimer (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, spin-rotational or spin-spin 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 electron-proton mass ratio can be determined with better accuracy, than it is known presently nist (). Moreover, any nonelectromagnetic long-range 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 long-range interactions between nuclei.

In this work we present an effective computational method with the use of exponential functions and numerical results for the Born-Oppenheimer 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 two-electron systems.

Ii Short range of internuclear distances

The Schrödinger equation for H molecule in the Born-Oppenheimer approximation is




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 James-Coolidge basis james_cool (), namely the functions of the form


such that


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 two-center 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 two-center two-electron Slater integral, see Ref. rec_h2 ().

For the presentation of analytic formulas we consider an integral


at the internuclear distance , since arbitrary distances can be obtained from the relation


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


Then, using the equation (69) from rec_h2 (), all nonvanishing integrals with are the following:


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 James-Coolidge basis is not appropriate as it does not include the Heitler-London function. Instead, we employ generalized Heitler-London functions, which are the product of a Heitler London function with an arbitrary polynomial in all electron distances,


such that


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


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




This one dimensional integration can easily be obtained with the sufficient accuracy and for example at the value of is


All integrals with higher powers of electron distances can be obtained from recursion relations, which were derived in Ref. rec_h2 (). Since


we can present formulae at . Let us introduce auxiliary functions


then all integrals with are the following


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 James-Coolidge 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 James-Coolidge 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.

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 James-Coolidge 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.

Table 4: Numerical values for BO potential at different internuclear distance . Results are obtained by extrapolation to complete set of basis functions

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

Figure 1: Rescaled exchange energy fitted to numerical points in Table IV in the range au, assuming Herring and Flicker herring () asymptotics, straight line.

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 Summary

We have demonstrated the applicability of analytic formulas for two-center two-electron integrals in high precision calculations of Born-Oppenheimer potential for the ground electronic state of the hydrogen molecule. The symmetric James-Coolidge basis with as much as 22 000 functions provided energies with relative precision of about for internuclear distances up to 10 au. The extended Heitler-London 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łos-Wolniewicz basis kw_basis ()).