# Two electron integral evaluation over primitive Gaussians without recurrence relations

I want to evaluate two-electron integrals over s, p, and d primitive Gaussian functions. The procedure is presented for s-type ones in "Modern Quantum Chemistry" by Szabo and Ostlund, but I need to evaluate these integrals for p and d ones. However, most books, such as "Molecular Electronic-Structure Theory" by Helgaker which I read, evaluate these integrals using recurrence relations. Is there any other method which can evaluate two-electron integrals, without using recurrence relations?

• You can use any numeric integration technique. It is, however, time-consuming, so people tried to reduce the time by exploiting properties of primitive gaussians. The relevant data can be found in either in papers on the subjects (see scholar.google.com and look for something like "primitive gaussian integral algorithm)", or in sources of relevant software (libint, gamess, and some other software). – permeakra May 12 '18 at 8:02
• @permeakra Thanks, but I look for an analytical solution. – Wisdom May 12 '18 at 8:34
• Sorry for being trivial, but... then use any method of analytical solution. The correct answer doesn't depend on the method, does it? – Ivan Neretin May 12 '18 at 10:51
• @Ivan Neretin The problem is that most of people use the recurrence relations in their codes and I have rarely seen analytical solution, specially for higher angular momentum than $s$ namely $p$ and $d$. Anyway I found a solution in a lecture note (still I'm not sure it's my desire answer) which I would post here, if no one post answer. – Wisdom May 12 '18 at 12:30
• Is there a specific reason why you do not want to use the recurrence relations? – Erik Kjellgren May 12 '18 at 21:05

As you mentioned, most implementations such as Obara–Saika (OS), McMurchie–Davidson (MD), Head-Gordon–Pople (HGP), and many others use a series of horizontal and vertical recurrence relations to build up angular momentum and transfer it among each of the four centers. These are analytic solutions, up until the need to solve the Boys function, however we could debate as to whether or not they are closed-form expressions by terminology. They are closed-form in the sense that each part of the recursive expansion does terminate properly with $(ss|ss)$, given in OS equation 44.

In general, any recursive function can be rewritten to be non-recursive using loops of sums and or products, though it may involve considerable effort. In that sense, the recurrence relationships are quite powerful and avoid lots of "boilerplate". If you are interested in a non-recursive expression, there is the earlier Taketa–Huzinaga–O-ohata (THO) paper that, while not very descriptive, includes explicit working expressions for 2-electron and 1-electron (overlap, kinetic energy, and nuclear-electron attraction) integrals. The most relevant expression is

\begin{align*} \text{ERI} &= \frac{2\pi^{2}}{\gamma_{1}\gamma_{2}} \left( \frac{\pi}{\gamma_{1} + \gamma_{2}} \right)^{1/2} \exp \left\{ - \frac{\alpha_{1}\alpha_{2}\overline{AB}^{2}}{\gamma_{1}} - \frac{\alpha_{3}\alpha_{4}\overline{CD}^{2}}{\gamma_{2}} \right\} \\ &\times \sum_{i_{1}i_{2}r_{1}r_{2}u} B_{i_{1}i_{2}r_{1}r_{2}u}(l_{1}, l_{2}, A_{x}, B_{x}, P_{x}, \gamma_{1} | l_{3}, l_{4}, C_{x}, D_{x}, Q_{x}, \gamma_{2}) \\ &\times \sum_{j_{1}j_{2}s_{1}s_{2}v} B_{j_{1}j_{2}s_{1}s_{2}v}(m_{1}, m_{2}, A_{y}, B_{y}, P_{y}, \gamma_{1} | m_{3}, m_{4}, C_{y}, D_{y}, Q_{y}, \gamma_{2}) \\ &\times \sum_{k_{1}k_{2}t_{1}t_{2}w} B_{k_{1}k_{2}t_{1}t_{2}w}(n_{1}, n_{2}, A_{z}, B_{z}, P_{z}, \gamma_{1} | n_{3}, n_{4}, C_{z}, D_{z}, Q_{z}, \gamma_{2}) F_{\nu}(\overline{PQ}^{2}/4\delta), \label{tho-2.22-1}\tag{THO eq. 2.22 part 1} \end{align*}

with

\begin{align*} &B_{i_{1}i_{2}r_{1}r_{2}u}(l_{1}, l_{2}, A_{x}, B_{x}, P_{x}, \gamma_{1} | l_{3}, l_{4}, C_{x}, D_{x}, Q_{x}, \gamma_{2}) \\ &\quad = (-)^{i_{2}} f_{i_{1}}(l_{1}, l_{2}, \overline{PA_{x}}, \overline{PB_{x}}) f_{i_{2}}(l_{3}, l_{4}, \overline{QC_{x}}, \overline{QD_{x}}) \\ &\qquad \times \frac{i_{1}!i_{2}!}{(4\gamma_{1})^{i_{1}}(4\gamma_{2})^{i_{2}}\delta^{i_{1} + i_{2}}} \cdot \frac{4\gamma_{1})^{r_{1}}(4\gamma_{2})^{r_{2}}\delta^{r_{1} + r_{2}}}{r_{1}!r_{2}!(i_{1}-2r_{1})!(i_{2}-2r_{2})!} \\ &\qquad \times [i_{1} + i_{2} - 2(r_{1} + r_{2})]! \frac{(-)^{u} p_{x}^{i_{1} + i_{2} - 2(r_{1} + r_{2}) - 2u} \delta^{u}}{u! \{ i_{1} + i_{2} - 2(r_{1} + r_{2} - 2u)\}!}, \label{tho-2.22-2}\tag{THO eq. 2.22 part 2} \end{align*}

and intermediate quantities are

$$\nu = i_{1} + i_{2} + j_{1} + j_{2} + k_{1} + k_{2} - 2 (r_{1} + r_{2} + s_{1} + s_{2} + t_{1} + t_{2}) - u - v - w,$$

\begin{align*} \gamma_{1} &= \alpha_{1} + \alpha_{2},~\gamma_{2}\text{ is similar}\\ \delta &= \frac{1}{4\gamma_{1}} + \frac{1}{4\gamma_{2}} \\ \mathbf{P} &= \frac{\alpha_{1} \mathbf{A} + \alpha_{2} \mathbf{B}}{\alpha_{1} + \alpha_{2}},~\mathbf{Q}\text{ is similar} \\ \mathbf{p} &= \mathbf{Q} - \mathbf{P} \\ \overline{AB} &= A - B \end{align*}

$$F_{\nu}(t) = \int_{0}^{1} u^{2\nu} \exp(-tu^{2}) \mathop{du},~\text{the Boys function.} \tag{THO eq. 2.11}$$

$(l, m, n)$ are the usual polynomial exponents leading to s-, p-, and d-type functions, $(\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D})$ are the usual basis function positions, the $\alpha$ are the basis function exponents, and $f_{j}(l, m, a, b)$ is the coefficient of $x^{j}$ in the expansion of $(x + a)^{l} (x + b)^{m}$.

In \eqref{tho-2.22-1} the summations with respect to the indices $i_{1}, i_{2}, r_{1}, r_{2}, u$ extend from 0 to $l_{1} + l_{2}$, $l_{3} + l_{4}$, $[i_{1}/2]$, $[i_{2}/2]$ and $[(i_{1} + i_{2})/2 - r_{1} - r_{2}]$ respectively. The ranges of $(j_{1}, j_{2}, s_{1}, s_{2}, v)$ or $(k_{1}, k_{2}, t_{1}, t_{2}, w)$ may be found in the same way.

Regarding the indices used in the summations:

1. If the index $i_{1}$ runs from 0 to $l_{1} + l_{2}$, and $i_{2}$ runs from 0 to $l_{3} + l_{4}$ we can deduce that $i_{1}$ must cover all possible values of the total angular momentum along the x-component for the bra, and $i_{2}$ goes up to the total angular momentum along the x-component for the ket. The indices $j,k$ match with $m,n$, meaning the y- and z-component.

2. $r,s,t$ are not to be confused with distances themselves, like how we usually define $r_{A} = |\mathbf{r}_{A}| = |\mathbf{r} - \mathbf{A}|$. Both $r,s,t$ and $u,v,w$ are also related to the three Cartesian coordinates, but in a more complicated manner than $i,j,k$; they appear to originate from the closed-form expressions for the Hermite polynomials.

Their implementation is slightly modified, with the necessary expressions in equations 3.3, 3.4, and 3.5, though the above can be implemented directly as well. A quick search revealed this code that implements all the THO integrals in C++. It is worth commenting that in real packages, the Boys function is precomputed from interpolation tables, rather than from the gamma and incomplete gamma functions.

• Thanks a lot. If I got it right, there is no way out of the Recurrence Relations, but the THO formula provides an escape! However I didn't understand what is the upper value of $i$ in summation of the expansion? – Wisdom May 13 '18 at 8:00
• I added a block about $i,j,k$; essentially they are related to the angular momentum of the x, y, and z components of each integral. Let me know if anything is still unclear. – pentavalentcarbon May 13 '18 at 15:04
• Thank you so much dear @pentavalentcarbon . In the case which I'm studying, all primitive Gaussian functions have been located at the origin. Can you tell me what happens to the THO formula with such assumption? – Wisdom May 14 '18 at 4:19
• Without writing it out on paper, I think $\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D}$ would all be zero (they are absolute positions, not differences), leading to $\mathbf{P}, \mathbf{Q}, \mathbf{p}$ being zero. That should get you started. – pentavalentcarbon May 14 '18 at 15:16