The functions Y^k(ab;r)

  By definition

$$\int_0^r \left(\frac{s}{r}\right)^k P(a;s)P(b;s)\,ds$$ Y^k(ab;r) = (61)

$$\quad + \int_r^\infty \left(\frac{r}{s}\right)^{k+1}P(a;s)P(b;s)\, ds .$$ (62)

Because integrals of this type occur frequently in Hartree-Fock calculations, it is desirable to determine them with maximum efficiency. Direct evaluation, using quadrature cannot be considered since an integral over the whole range would be required for each value of r. These functions can be determined much more rapidly as solutions of a pair of differential equations as suggested by Hartree (1957).

Define

$$Z^{k}(ab;r) = \int_0^r \left(\frac{s}{r}\right)^kP(a;s)P(b;s)\,ds .$$ (63)

Then it is easy to show that

$$\frac{d}{dr}Z^{k}(ab;r) = P(a;r)P(b;r) -\frac{k}{r}Z^{k}(ab;r)$$ (64)

and

$$\frac{d}{dr} Y^{k}({a}{b};r) = \frac{1}{r}\left[(k+1) Y^{k}({a}{b};r) -(2k+1)Z^{k}(ab;r)\right]$$ (65)

with the boundary conditions

$$Z^{k}(ab;0) = 0 \mbox{ and }Y^{k}({a}{b};r) \rightarrow Z^{k}(ab;r) \mbox{ as } r \rightarrow\infty .$$ (66)

With the $\rho$ variable, the equations for Z^k and Y^ktransform to

$$\frac{d}{d\rho}Z^k = r \overline{P}_a \overline{P}_b - k Z^k$$ (67)

$$\frac{d}{d\rho}Y^k = (k+1)Y^k - (2k+1)Z^k .$$ (68)

As was shown by Worsley (1958), both these equations have simple integrating factors, namely $e^{k\rho}$ and $e^{-(k+1)\rho}$, respectively, so that

$$d\left(e^{k\rho}Z^k\right) = r \overline{P}_a \overline{P}_be^{k\rho}d\rho$$ (69)

$$d\left(e^{-(k+1)\rho}Y^k\right) = - (2k+1)Z^ke^{-(k+1)\rho } d\rho .$$ (70)

On integrating these equations we get

$$Z^k_{j+m} = e^{-mhk}Z^k_j + \int_{\rho_j}^{\rho_{j+m}}r\overline{P}_a \overline{P}_b e^{k(\rho - \rho_{j+m}} d\rho$$ (71)

$$j = 1,2,\ldots$$ (72)

$$Y^k_j = e^{-mh(k+1)}Y^k_{j+m} + (2k+1)\int_{\rho_j}^{\rho_{j+m}}Z^k e^{-(k+1)(\rho - \rho_j)} d\rho$$ (73)

$$j = \ldots, 2, 1 .$$ (74)

Note that the calculation of Z^k proceeds by outward integration using known starting values whereas Y^k is obtained by inward integration. Here we have used the fact that $\rho_{j+m} -\rho_j = mh$. In the HF program m=2 is used; the first and last pair of intervals are integrated using Simpson's rule but for other pairs of intervals the more accurate ${\cal O}(h^7$) formula is used. For k >0 the method is stable in that errors at previous points are damped by an exponential factor, but for k=0, errors accumulate for the outward integration of Z^kthough the errors in Y^k are again damped.

The Y^k functions can also be thought of as solutions of a boundary value problem. Differentiating the first-order equation for Y^k and using the equation for Z^k to eliminate the latter, we obtain

$$\frac{d^2}{dr^2}Y^{k}({a}{b};r) = \frac{k(k+1)}{r^2} Y^{k}({a}{b};r)-\frac{(2k+1)}{r} P(a;r)P(b;r)$$ (75)

with the boundary conditions

$$Y^k(ab;0) = 0 \mbox{ and }\frac{d}{dr}Y^{k}({a}{b};r) = - \frac{k}{r}Y^{k}({a}{b};r) \mbox{ as }r \rightarrow \infty.$$ (76)

This equation is not stable for outward integration but has been used successfully in basis type methods (Brage and Froese Fischer 1994).

2001-01-09