The logarithmic grid

Numerical methods are greatly simplified when the grid points defining the discretized problem are equally spaced with respect to the independent variable. At the same time, in order to achieve a desired accuracy, the points must be sufficiently close together that a polynomial passing through adjacent points (r_j, P_ja) is a reasonable approximation to the function P(a;r). In regions where the function is changing rapidly, or where some higher derivative is large, the points must be close together.

The natural choice for the independent variable is the radius r. But in terms of this variable the radial functions change fairly rapidly near the origin, which is also a singular point. As a result, it is necessary to have points close together near the origin and further apart for larger r. Early methods would use equal step sizes hbetween points with an occasional doubling of step size. The process of doubling h is not difficult but it does increase the complexity of a numerical algorithm. A similar result may be achieved analytically by transforming to a new variable.

A transformation that has been used successfully in the present atomic structure program is

$$\rho = \log (Zr) \Rightarrow r = \frac{e^{\rho}}{Z}$$ (45)

along with the transformation of the dependent variable

$$\overline{P}(\rho) = P(r)/\sqrt{r}.$$ (46)

With this transformation the Hartree-Fock equations become

$$\left[ \frac{d^2}{d\rho^2} + 2Zr - (l_a+\frac{1}{2})^2-r^2\varepsilon_{aa}\right]\overline{P}(a;\rho) 2r\left[ Y(a;\rho)\overline{P}(a;\rho)+ \overline{X}(a;\rho)\right]$$ =

$$+ \sum_{b\ne a} \delta_{l_al_b} r^2\varepsilon_{ab}\overline{P}(a;\rho)$$ (47)

with the boundary conditions $\overline{P}(a;\rho) \rightarrow 0$as $\rho \rightarrow \pm \infty$. Here $\overline{X}(a;\rho)) =X(a;r)/\sqrt{r}$.

This transformation has two advantages. As will be shown in the next subsection, the equations defining Y^k(ab;r) now have a simplified form and can be computed more rapidly. For large atoms this is extremely important since the number of exchange terms increases rapidly with the number of functions. Another advantage arises from the fact that to zero order within perturbation theory, the radial functions are hydrogenic functions with Zras independent variable. The transformation makes it possible to use the same set of $\rho$ values for all atoms in the periodic table. Only in exceptional circumstances is the grid changed.

A disadvantage of the transformation is the fact that now the range $(0,\infty)$ is transformed to $(-\infty, \infty)$. In order to avoid the point at $-\infty$ it is customary to divide the range of r into two region, (0,r₁) and $(r_1,\infty)$and transform only the latter to the $\rho$ variable. The value of r₁ is chosen so that $\rho_1 = \log(Zr_1)$ is a predetermined number, sufficiently small so that a series expansion of three terms may be used near the origin. Numerous tests on the hydrogen equation suggest that appropriate values are $\rho_1 = -4$ and h=1/16 for an accuracy of 9 digits in a hydrogenic 1s orbital energy. Thus the grid points become

$$\rho_j = -4 +(j-1)/16 \mbox{\quad or \quad } r_j = e^{\rho_j}/Z, j =1,2,\ldots, M.$$ (48)

In a many-electron system, since the energies of other electrons are smaller, not only because of the principal quantum number but also because of the smaller $Z_{\mbox{eff}}$, the total energy is expected to be limited by the accuracy of the 1s orbital.

In the interval (0,r₁) all functions are approximated by a series expansion about the origin. Some important series are

$$\mbox{$P(nl;r)$ } a_0(nl) r^{l + 1} \left\{ 1 -\frac{Zr}{l+1} +\alpha r^2 + {\cal O}(r^3)\right\}$$ = (49)

$$vr + {\cal O}(r^3)$$ Y(nl;r) = (50)

$$a_0(nl)wr^{l+2} + {\cal O}(r^{l+3})$$ X(nl;r) = (51)

as $r \rightarrow0$ where

$$\alpha = \frac{Z^2 + (l+1)( \varepsilon/2 + v + w)}{(2l+3)(l+1)}$$ (52)

$$\varepsilon= \varepsilon_{nl,nl} + \sum_{n^\prime l \ne nl} \varepsilon_{nl,n^\prime l}\frac {a_0(n^\prime l)}{a_0(nl)}.$$ (53)

The constants v and w depend on the configuration and a₀(nl) is the AZ(nl) parameter included in the program. It is clear from these expansions that the least accurate case always occurs for l=0, where the series converge least rapidly.

2001-01-09