Solution of the radial equation
The radial equation for a single electron outside a core is a linear integro-differential equation of the form

$$y^{\prime\prime} + \left[f(r) -\varepsilon\right]y(r) = \int_0^\inftyK(r,s)\,y(s)\,ds$$ (77)

where y(0)=0 and $y(r)\rightarrow 0$ a $r \rightarrow\infty$. The integral equation part includes the exchange effect as well as off-diagonal energy parameters, if any. A direct discretization would lead to an eigenvalue problem

$${\bss (T -\varepsilon I -K)y} = 0$$ (78)

where ${\bss T}$ is a band matrix (the band width depending on the order of discretization of $y^{\prime\prime}$) and ${\bss K}$is a full matrix related to the kernel K(r,s), but also including factors arising from a quadrature rule. This eigenvalue problem could be solved for a selected eigenvalue using Rayleigh-Quotient iteration, but would then require the storage of the full matrix and ${\cal O}(M^3)$ floating point operations per iteration. An iterative procedure requiring only ${\cal O}(M)$ operations per iteration is used by the HF program.

Let {y^(k)(r), k = 0, 1, ,,} be a sequence of approximate solutions. Define

$$g^{(k)}(r) = \int_0^\infty K(r,s)\, y^{(k)}(r)\,ds$$ (79)

and

$$\left\{\frac{d^2}{dr^2} + f(r) -\varepsilon\right\}y^{(k+1)}(r) =g^{(k)}(r).$$ (80)

For the latter, the familiar Numerov method based on the first two terms of the discretization

$$\delta^2 y_j = h^2\left(1 +\frac{1}{12}\delta^2 -\frac{1}{240......+ \frac{31}{60480} \delta^6 - \cdots\right)y^{\prime\prime}_j$$ (81)

leads to a system of equations

$${\bss (T -\varepsilon S) y}^{(k+1)} = {\bss c}^{(k)}$$ (82)

where ${\bss T} = (t_{ij})$ and ${\bss S} = (s_{ij})$ are $(M-1)\times(M-1)$ tridiagonal matrices with

$$t_{ij} = \left\{ \begin{array}{ l l}-2 + (10/12)h^2f_i, & i...... 0, & \mbox{otherwise}\end{array} \right.f_j \equiv f(r_j)$$

and

$$s_{ij} = \left\{ \begin{array}{ l l}(10/12)h^2, & i=j\\(......t i-j\vert = 1\\0, & \mbox{otherwise}.\end{array} \right.$$

The column vector ${\bss c} = (c_1, c_2, \ldots, c_{M-1})^t$ has components

$$c_i = \frac{h^2}{12}(g_{i+1} + 10 g_i + g_{i-1}), \quad g_i\equiv g^{(k)}(r_i)$$ (83)

For simplicity the above assumes boundary conditions of zero at both ends.

Unlike eq:radial_integ, eq:radial_iter is no longer an eigenvalue problem and unique solutions to the discretized approximations of eq:radial_disc will exist for values of $\varepsilon$ if ${\bss (T -\varepsilon S)}$is non-singular. If the coefficient matrix is singular, unique solutions may still exist if either ${\bss c}^{(k)} =0$ or if ${\bss c}$ is orthogonal to the desired eigenvector. In either case, for the solutions of eq:radial_disc to converge to those of eq:radial_integ it is necessary to define $\varepsilon= \varepsilon^{(k)}$so that the sequence $\varepsilon^{(k)}, k= 0, 1, \ldots$ converges to an eigenvalue. A natural choice is the Rayleigh quotient for y^(k)(r). At the same time, we would like the iterations to converge to normalized eigenfunctions where normalization is defined in terms of numerical integration.

Define ${\bss (u,v)}$ to be the numerical approximation to $\int_0^\infty u(r)v(r)\,dr$ using only vectors ${\bss u}$ and${\bss v}$ representing values of the functions on the grid. Also define ${\vert\vert\bss u}\vert\vert^2 = ({\bss u,u })$. The method M1 implemented in the HF program combines these ideas as follows:

For $k=0, 1, 2, \dots$

1.

$\varepsilon^{(k)} = {\bss y}^{(k)}{^t}\left[{\bss T}{\bss y}^{(k)}-{\bss c}^{(k)}\right]/{\bss y}^{(k)}{^t}{\bss S}{\bss y}^{(k)}$

2.

$\left[ {\bss T}- \varepsilon^{(k)}{\bss S}\right]{\bss z}^{(k+1)}={\bss c}^{(k)}$

3.

${\bss y}^{(k+1)}= {\bss z}^{(k+1)}/\vert\vert{\bss z}^{(k+1)}\vert\vert$

This method works remarkably well in many cases but when the exchange effect is extremely small, the coefficient matrix is too nearly singular and a second method M2 is used. The latter, as ${\bss c}\rightarrow 0$ goes over to inverse iteration, a method for finding the eigenvector for the eigenvalue closest to $\varepsilon^{(k)}$.

Method M2 proceeds very much as M1, but also introduces an inverse iteration step

$$\left[ {\bss T}- \varepsilon^{(k)}{\bss S}\right]{\bss w}^{(k+1)}= {\bss S}{\bss y}^{(k)}.$$ (84)

(The latter is the variational equation for $w(r) = \partial y(r)/\partial \varepsilon$ under the assumption that f(r) and g(r) do not depend on the energy). Then the next estimate is formed as

$${\bss y}^{(k+1)}= {\bss z}^{(k+1)}+ \beta{\bss w}^{(k+1)}; \q......eta \mbox{such that }\vert\vert{\bss y}^{(k+1)}\vert\vert=1 .$$ (85)

The HF and MCHF programs solve the radial equation on the logarithmic grid not directly by linear algebra. Instead outward and inward integration are used, somewhat for historical reasons. Outward integration is stable in the oscillatory region where matrix methods would require pivoting. In the tail region the matrix is diagonally dominant and pivoting is not needed. In this region a system of equations is solved that matches the solution of outward integration, at the same time looking for a point where the boundary condition at large r can be applied. Thus the algorithm is ``adaptive'' and the user need not be concerned with the extent of the radial orbital. In both regions a particular solution and solution of the homogeneous equation are obtained. These are then combined so that the three-term recurrence at the join is satisfied.

In order to ensure that the iterative process converges to the desired eigenfunction, node counting is used and the notion of an ``acceptable'' solution is applied. Particularly in the early iterations, the solutions obtained by both M1 and M2 may be ``unusual''. In such cases, the Rayleigh quotient value of $\varepsilon$ must be replaced by some other value. A somewhat ad hoc search for an acceptable solution is initiated. In outward/inward integration, the three-term recurrence is satisfied at all points except the join. In some instances, it has proven useful to relate the residual of the equation at the join to the change in energy required to satisfy this condition.

When orbitals are multiply occupied, the function f(r) also depends on y(r) and the problem is non-linear. In such cases, the SCF iterations may oscillate and accelerating techniques may need to be used as described earlier.
 

2001-01-09