Introduction

One thing that all the calculations in this thesis have in common is their reliance on the calculation of a number of two-electron repulsion integrals (ERIs) (for example, equation (1.24)) over contracted Gaussian-type basis functions. As mentioned in section (1.5.7), the number of ERIs grows as O(N²) and has become the computational bottleneck for most HF and DFT calculations. Because of this, the efficient calculation of ERIs has been the focus of much research [56,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176]. Even the new O(N) approaches for construction of the Coulomb matrix require the computation of ERIs for short-range interactions.

It should be noted that almost all ERIs calculated at present are over Gaussian type orbitals. This is entirely due to the Gaussian product rule: that the product of two Gaussian functions is another Gaussian, centered somewhere on the line connecting the first two centers. Thus, any four-center ERI involving GTOs can immediately be reduced to a two-center integral (over GTOs), a simplification not available with STOs. Thus the fundamental (contractionless and momentumless) integral
$$\begin{align}I=\iint e^{-\alpha\vert{\bf r}_{1}-{\bf A}\vert^{2}} e^{-\beta\vert... ...delta\vert{\bf r}_{2}-{\bf D}\vert^{2}} \, d{\bf r}_{1} d{\bf r}_{2} \end{align}$$
becomes
 $$\begin{align} I = U \iint e^{-\zeta\vert{\bf r}_{1}-{\bf P}\vert^{2}} \frac{1}{r... ...\eta\vert{\bf r}_{2}-{\bf Q}\vert^{2}} \, d{\bf r}_{1} d{\bf r}_{2}, \end{align}$$
where
$$\begin{align}U = \exp\left[\frac{-\alpha\beta}{\alpha + \beta}\vert{\bf A}-{\bf ... ...c{-\gamma\delta}{\gamma+\delta}\vert{\bf C}-{\bf D}\vert^{2}\right] \end{align}$$
and
$$\begin{align}\zeta &= \alpha + \beta& \eta &= \gamma + \delta \\ \mathbf{P} &... ...mathbf{Q} &= \frac{\gamma {\bf C}+ \delta {\bf D}}{\gamma + \delta}. \end{align}$$

How equation (5.2) is then computed depends on the integral method used. Each method has a number of subtle peculiarities, but all involve four basic steps (although the order of action varies from method to method).

The operator step, O, generates the momentumless two-center integrals over the two-electron operator. In the PRISM method [160] (which includes the McMurchie-Davidson [163], Obara-Saika [164], Head-Gordon-Pople [165], and Ten-no [166] as special cases) this forms the [0]^(m) integrals.

In the momentum step L, recursive identities are used to build the momentumless quantities into those of the required angular momenta. This is where the bulk of integral research has been focused over the last two decades [159,160,161,162,163,164,165,166,169,170,171,172,173,174,175].

The contraction step gathers together the primitive (uncontracted) components to form fully contracted contributions. The main advance of the PRISM method was the execution of the contraction step at the most optimal time. Previously, the contraction step has always been after the operator step (in the case of PRISM, forming the (0)^(m)s). The work presented in this chapter describes how contraction can be carried out first, thus providing the prospect of massive computational savings, and also forming the CO path to the COLD PRISM [176].

In the density step D, the contracted quantities are multiplied by the density matrix elements $P_{\mu \nu}$. This step is usually termed `digestion' of the integrals and was traditionally performed last. One of the discoveries of the COLD PRISM, however, is that large time savings can be achieved by performing this step earlier.