Dirac-Fock atomic electronic structure calculations using different nuclear charge distributions

L. Visscher and K.G. Dyall, Dirac-Fock atomic electronic structure calculations using different nuclear charge distributions, Atom. Data Nucl. Data Tabl., 67, (1997), 207.

Numerical Hartree-Fock calculations based on the Dirac-Coulomb Hamiltonian for the first 109 elements of the periodic table are presented. The results give the total electronic energy as a function of the nuclear model that is used. The main purpose of the calculations is to provide a set of benchmark results for basis set calculations that use a finite size nuclear model.

INTRODUCTION

The use of finite size nuclear models in molecular relativistic electronic structure theory is becoming popular. Apart from presenting a more physical model for the nuclear charge distribution than the usual point charge model, the use of a finite size model has also important advantages in calculations that solve the Hartree-Fock or Dirac-Fock equations in the algebraic approximation.

An expansion in a finite basis is preferably done using Gaussian-type functions, since multi-center integrals over Gaussian-type functions are more readily evaluated than integrals over Slater functions. This computational advantage outweighs the slower convergence of the atomic energy with the number of basis functions. It remains, however, essential for computational feasibility to keep the number of basis functions employed as small as possible. Particularly in relativistic calculations a large number of basis functions is required to accurately represent the electron spinors in the core region. This slow convergence is due to the weak singularity at the nucleus for the 1s and 2p spinors. It takes many Gaussian-type functions to model this behaviour and get an energy close to the Dirac-Fock limit. As this weak singularity is a consequence of the point charge nuclear model, one observes that the convergence is better¹ when employing a finite size nuclear model.

The quality of an expansion set is usually judged from the errors made in atomic calculations. These errors can be determined by comparing the results obtained in the algebraic approximation with accurate finite difference results. For the methods that employ a finite nuclear model such a set of reference figures has not yet been tabulated in a consistent way. We here present a series of calculations that can serve this purpose and lead to standarization in the parameters used by different groups in the field.

We have done four sets of calculations using both the point nucleus approximation and three popular finite size nuclear models as they are implemented in the freely distributed atomic structure code GRASP². Below we summarize the important features of these models and define the formulae used to calculate the necessary parameters. In these formulae Z is the atomic number (number of protons) of the nucleus, A is the mass number and R and r represent the nuclear and electronic radial coordinate, respectively. The root-mean-square radius (RMS) of the nucleus that is used in these models can be approximately related to the isotope number via the empirical formula³

Note that the isotope or atomic mass number rather than the actual mass is used in this formula. The accuracy³ of this fit function relative to the measured RMS values is better than 0.05 fm. We chose to use this fit function rather than the actual measured RMS radii to keep the number of parameters in the models minimal and to facilitate computational implementation.

Model 0: Point Charge.

For completeness we start by giving results for the point charge approximation. The potential is in this case solely determined by the nuclear charge Z.

Model 1 :The homogeneously charged sphere.

Here we assume that the charge of the nucleus is uniformly distributed in a sphere of radius R₀. This simple model gives a harmonic potential within the sphere, which means that the exact solutions near the origin will be Gaussian-like. This is an obvious advantage when using a Gaussian expansion for the spinors. Integrals over such functions are relatively easy to calculate and the model is therefore also used in molecular electronic structure calculations^4-7. A compilation of atomic calculations using this type of model has been given earlier by Desclaux⁸ employing a different formula for R₀. R₀ is easily related to the RMS in formula (1) via

The nuclear charge density is given by

and the potential can be written in analytical form as

It is clear that the potential from a uniform nucleus has a discontinuity in its second derivative at the nuclear boundary, and as a consequence care must be taken in any numerical integration scheme. The present version of GRASP, used in this work, ensures that the nuclear boundary is a grid point, thereby eliminating numerical instabilities. The convergence with the numerical grid spacing remains slower, however, than with the other models.

Model 2 : The Fermi 2-parameter charge distribution.

A more physical description of the nuclear structure is given by a Fermi-type distribution³. Here the potential is defined by the "half charge radius" C and a skin thickness parameter T. C is defined as the value of R for which rho(R) = 0.5 rho₀, while T gives the interval in which rho(R) falls from 0.9 rho₀ to 0.1 rho₀. T is chosen to be 2.30 fm for all nuclei. C can be related to R₀ via the formula

For atomic masses < 5 mu (atomic mass units), this formula can not be applied. We then use the formula

The nuclear charge density for the Fermi distribution is

where is determined from the total charge. The potential cannot be expressed in simple analytic form.

Model 3 : The Gaussian charge distribution.

This model is popular due to the easy computational implementation. Since the electronic wave function is usually also expanded in Gaussian-type functions one may evaluate both nuclear-electron attraction and electron-electron repulsion integrals by the same efficient primitive integral routines. Different choices for the exponential parameter are in use¹. In order to come to a standarization we propose to use the model given below. In this parametrization the exponential parameter, xi, is obtained by equating the RMS value of the Gaussian radial distribution to the empirical formula in Eq. (1).

(11)

(12)

(13)

The potential V is related to the error function via⁹

(14)

We have done a series of reference Dirac-Fock calculations based on the four nuclear charge density models described above. The Hamiltonian that we used is the Dirac-Coulomb Hamiltonian. The speed of light, in Hartree atomic units, is taken equal to 137.0359895, the conversion factor a₀ to fm used is 52917.7249. As the value for the mass number, A, we chose the most abundant or most stable isotope of the element. In table I the resulting values of , C and xi are given in atomic units.

Radial electronic wave functions were obtained using the Average Level (AL) scheme in GRASP. Details about this scheme can be found in the description of the GRASP program². For most elements the electronic energies calculated are the weighted average energies of the ground state configurations. For the transition metals and some of the lanthanides and actinides, we also calculated the weighted average energies of other low-lying configurations. The logarithmic grid spacing for the solution of the radial equations was chosen to as 0.01 ln(a₀). For a nuclear charge Z=10 the total electronic energy as a function of this computational parameter is converged to 10^-9 E_h (where E_h is the atomic unit Hartree) for all the models. For Z=100 a precision of 10^-7 E_h in the total electronic energy is obtained for all models, except the uniform charge distribution where the energy is accurate to 10^-5. The energies are given in table II, displayed to the accuracy reached in the uniform charge distribution model.

Spinor energies and radial expectation values of the spinors that were calculated using the different models for the nuclear charge distribution are available as supplementary data. Since they comprise a large quantity of data, they are presented in electronical html-format via the reference "http://dirac.chem.sdu.dk/doc/FiniteNuclei/FiniteNuclei.shtml". Normalized differences between the models as a function of the nuclear charge are shown in figs. 1 and 2 for the total electronic energy and the binding energy of the 1s spinor, respectively. The analogous plots for the other spinors would look very similar as well.

Figure 1. Deviations in the finite nuclear size correction to the total electronic energy. Plotted is the difference [E(Model) - E(Fermi model)] / [E(Fermi Model) - E(Point Nucleus model)].

Figure 2. Deviations in the finite nuclear size correction to the 1s binding energy. Plotted is the difference [E(Model) - E(Fermi model)] / [E(Fermi Model) - E(Point Nucleus model)].

From the Tables it is clear that all finite-size nucleus models give a higher total energy than the point-charge model. The differences between the different finite size models are small, for the lighter elements the Uniform and the Gaussian charge distributions give similar energies; for the heavier elements the influence of the nuclear skin becomes smaller and the Uniform and Fermi energies are very similar. Most of the energy difference with the point-charge model arises from the change in energy of the innermost s-spinors. The influence of the chosen model on the spin-orbit splitting of the levels is negligible. For Z=100 the largest difference between the point charge and distributed charge models is 0.2 % in the spin-orbit splitting of the 2p shell.

We hope that these values may serve as a reference for future (basis set) calculations based on the Dirac-Coulomb Hamiltonian.

Acknowledgements

The authors like to thank consulting editor Dr. W. Scholz for constructive remarks.

LV was supported by a National Research Council postdoctoral fellowship and a Cray Research Grant. KGD was supported by a prime contract NAS2-14031 from NASA to Eloret and by contract BPNL 291140-A-A3 from Pacific Northwest Laboratory to Eloret. Pacific Northwest Laboratory is a multi-program laboratory operated by the Batelle Memorial Institute for the U. S. Department of Energy under contract DE-AC06-76RLO-1830; the work was performed under the auspices of the Office of Scientific Computing and under the Office of Health and Environmental Research, which funds the Environmental and Molecular Sciences Laboratory Project, D-384.

REFERENCES

¹O. Visser, P. J. C. Aerts, D. Hegarty, and W. C. Nieuwpoort, Chem. Phys. Lett. 134, 34 (1987).

²K. G. Dyall, I. P. Grant, C. T. Johnson, E. P. Plummer, and F. Parpia, Comp. Phys. Comm. 55, 425 (1989).

³W. R. Johnson and G. Soff, Atomic Data Nuclear Data Tables 33, 405 (1985).

⁴Y. Ishikawa, R. Baretty, and R. C. Binning Jr., Chemical Physics Letters 121, 130 (1985).

⁵A. K. Mohanty and E. Clementi, Int. J. Quantum Chem. 39, 487 (1991).

⁶A. K. Mohanty and E. Clementi, Int. J. Quantum Chem. 40, 429 (1991).

⁷O. Matsuoka, Chemical Physics Letters 140, 362 (1987).

⁸J. P. Desclaux, Atomic Data Nuclear Data Tables 12, 311 (1973).

⁹K. G. Dyall and K. Faegri Jr., Chemical Physics Letters 201, 27 (1993).

EXPLANATION OF TABLES

Table I. Isotope numbers, Nuclear Radii and Nuclear Model Exponents

This table presents the approximate nuclear radii used to define the Uniform Sphere, Fermi and Gaussian nuclear charge distribution models. For the latter model the exponential parameter xi is also given.

Z Atomic number

A Isotope number

RMS Radius of the nucleus for the finite nucleus models (in atomic units).

C Half charge radius in the Fermi Model (in atomic units).

xi Exponential parameter in the Gaussian Model (in atomic units).

Table II. Electronic configuration and total electronic energies for the different Nuclear Models

This table gives the total electronic energies for the atoms using four different nuclear models.

Z Atomic number

Electronic Configuration Electronic configuration of the atom given in abbreviated form. Only the groundstate or lowlying excited configurations are calculated.

Electronic energy Electronic energy (in Hartrees) calculated using the Nuclear Model indicated on the next line. The spinors are optimized by minimizing the average energy of all states arising from the configuration specified.

Table III. Orbital energies for the different Nuclear Models

This table gives the orbital binding energies for the atoms using four different nuclear models.

Z Atomic number

Electronic Configuration Electronic configuration of the atom given in abbreviated form.

Orbital energy Orbital energy (in Hartrees) calculated using the Nuclear Model indicated on the next line.

Table IV. Radial expectation values <r> for the different Nuclear Models

This table gives the radial expectation values of r of the orbitals calculated using four different nuclear models.

Z Atomic number

Electronic Configuration Electronic configuration of the atom given in abbreviated form.

< r > Expectation value of r (Bohr) for the different orbitals calculated using the Nuclear Model indicated on the next line.

Table V. Radial expectation values <r^-1> for the different Nuclear Models

This table gives the radial expectation values of r^-1 of the orbitals calculated using four different nuclear models.

Z Atomic number

Electronic Configuration Electronic configuration of the atom given in abbreviated form.

<r^-1> Expectation value of r^-1 (Bohr) for the different orbitals calculated using the Nuclear Model indicated on the next line.