Traditional (0)^(m) Generation

Momentum is built up by recurrence relations which use (L+1) K_Bra K_Ket (where K_Bra is the number of Gaussians on center ${\bf A}$ times the number at center ${\bf B}$, and similarly for K_Ket with centers ${\bf C}$ and ${\bf D}$) [0]^(m)s (one m for each degree of momentum, L).

The simplest [0]^(m) (that is, [0]⁽⁰⁾) is the fundamental electron-repulsion integral, equation (5.2). The integration is performed by first replacing each of the three factors by its Fourier representation, giving (after some re-ordering)
$$\begin{align}[0]^{(0)} = &(2\pi)^{-3} U \left(\frac{\pi^{2}}{\zeta\eta}\right)^{... ...-{\bf k}_{2})}\,d{\bf r}_{2}\, d{\bf k}_{1}d{\bf k}_{2}d{\bf k}_{3}. \end{align}$$
The last two integrals above are Fourier representations of the three-dimensional Dirac delta function, allowing simplification to
$$\begin{align}[0]^{(0)} = \frac{U}{8(\zeta\eta)^{3/2}} \iiint e^{-k_{1}^{2}/4\zet... ...elta({\bf k}_{3}-{\bf k}_{2})\,d{\bf k}_{1}d{\bf k}_{2}d{\bf k}_{3}. \end{align}$$
Remembering that the Dirac $\delta$ function has the sampling property of
$$\begin{align}\int \delta({\bf r}_{1}-{\bf r}_{2})h({\bf r}_{1})\, d{\bf r}_{1} = h({\bf r}_{2}) \end{align}$$
for any function $h({\bf r})$ allows the triple integral to collapse to the single integral
$$\begin{align}[0]^{(0)} = \frac{2 \pi^{2} U}{(\zeta\eta)^{3/2}R} \int_{0}^{\infty} \frac{\sin u}{u} e^{-u^{2}/(4\mathrm{T})} \, du \end{align}$$
where, for convenience and computational efficiency, the following variables have been introduced:
$$\begin{align}{\bf R}&= {\bf Q}- {\bf P}\\ \mathrm{T} &= \upsilon^{2} R^{2} \\ \upsilon^{2} &= \frac{\zeta\eta}{\zeta + \eta}, \end{align}$$
and U is scaled by $(\frac{\pi^{2}}{\zeta\eta})^{3/2}$. The final integration is usually performed with the introduction of a special function G_m(T),
$$\begin{align}G_{m}(\mathrm{T}) = \sqrt{\frac{2}{\pi}} \int_{0}^{1} t^{2m} e^{-\mathrm{T}t^{2}} \, dt \end{align}$$
giving
 $$\begin{align}[0]^{(0)} = U \sqrt{2 \upsilon^{2}} G_{0}(\mathrm{T}). \end{align}$$

The higher momentum [0]^(m) formulae can be found via the relation
$$\begin{align}[0]^{(m+1)} = \frac{-d[0]^{(m)}/dR}{R} \end{align}$$
producing the general formula
$$\begin{align}[0]^{(m)} = U \left( 2 \upsilon^{2} \right)^{m+1/2} G_{m}(\mathrm{T}). \end{align}$$
If the value of T is less some critical value, the function G_L(T) is calculated using Chebyshev interpolation [167]. The G_m(T), $0 \le m < L$ are then calculated by downward recursion using
$$\begin{align}G_{m}(\mathrm{T}) = \frac{1}{2m+1} \left[ 2 \mathrm{T} G_{m+1}(\mathrm{T}) + e^{-\mathrm{T}} \right]. \end{align}$$
Note that this also requires an interpolation for e^-T.

If T is greater than a critical value, the distributions overlap negligibly, and G_m(T) can be approximated
$$\begin{align}G_{m}(\mathrm{T}) \approx \frac{(2m-1)!!}{(2\mathrm{T})^{m+1/2}} \end{align}$$
and the [0]^(m) reduces to the classical multipole formula
$$\begin{align}[0]^{(m)} = \left( \frac{U}{R} \right) \left( \frac{1}{R^{2}} \righ... ...( \frac{3}{R^{2}} \right) \ldots \left( \frac{2m-1}{R^{2}} \right ), \end{align}$$
which can be computed extremely fast via recursion.

The (0)^(m)s are formed if contraction is carried out immediately after forming the primitives. Contraction can occur early or late, depending on the PRISM path; however, it is advantageous to contract early if the total degree of contraction ( K_tot = K_BraK_Ket) is high. With the new algorithms mainly reducing the cost of high momentum, the point is fast being approached where ERI generation is dominated by making and contracting the [0]^(m)s. During contraction, scalings are introduced which are required for the later transformations, so a (0)^(m) is more accurately represented as:
$$\begin{align}{}_{a'b'p'}(0)_{c'd'q'}^{(m)} = \sum^{K_{Bra}}\sum^{K_{Ket}} \frac{... ...]^{(m)} \frac{(2\gamma)^{c'}(2\delta)^{d'}}{(2\gamma+2\delta)^{q'}}. \end{align}$$
The double sum above is the problem. The amount of work required to calculate (0)^(m) scales with K_BraK_Ket. Note that this means two interpolations and a divide are calculated K_BraK_Ket times! The rest of this chapter presents a method which scales independently of K_Bra and K_Ket.