Recurrence relations for four-electron integrals over Gaussian basis functions

TL;DR

This paper develops recurrence relations for efficient computation of four-electron integrals over Gaussian basis functions, extending previous three-electron algorithms to handle more complex integrals in quantum chemistry.

Contribution

It introduces a generalized recursive approach for four-electron integrals, applicable to various operators in explicitly-correlated F12 methods, enhancing computational efficiency.

Findings

Derived vertical, transfer, and horizontal recurrence relations

Applicable to a broad class of four-electron operators

Facilitates accurate and efficient integral calculations

Abstract

In the spirit of the Head-Gordon-Pople algorithm, we report vertical, transfer and horizontal recurrence relations for the efficient and accurate computation of four-electron integrals over Gaussian basis functions. Our recursive approach is a generalization of our algorithm for three-electron integrals [J.~Chem.~Theory Comput.~12, 1735 (2016)]. The RRs derived in the present study can be applied to a general class of multiplicative four-electron operators. In particular, we consider various types of four-electron integrals that may arise in explicitly-correlated F12 methods.

Tables

Table 1. Table 1: Kernels F ​ ( t ) 𝐹 𝑡 F(t) of the Gaussian integral representation for various f 12 subscript 𝑓 12 f_{12} operators. δ ​ ( x ) 𝛿 𝑥 \delta(x) and θ ​ ( x ) 𝜃 𝑥 \theta(x) are respectively the Dirac delta and Heaviside step functions, and erf ⁡ ( x ) error-function 𝑥 \erf(x) and erfc ⁡ ( x ) erfc 𝑥 \operatorname{erfc}(x) are the error function and its complement version, respectively. Olver  et al. ( 2010 )

$f_{12}$	$F\left(t\right)$
$1$	$\delta \left(t\right)$
$r_{12}^{-1}$	$2/\sqrt{\pi }$
$r_{12}^{-2}$	$2t$
${\left(r_{12}^2+{\lambda }^2\right)}^{-1/2}$	$\left(2/\sqrt{\pi }\right)\exp \left(-{\lambda }^2{t}^2\right)$
$\exp \left(-\lambda r_{12}\right)$	$\left(\lambda t^{-2}/\pi \right)\exp \left(-{\lambda }^2{t}^{-2}/4\right)$
$r_{12}^{-1}\exp \left(-\lambda r_{12}\right)$	$\left(2/\sqrt{\pi }\right)\exp \left(-{\lambda }^2{t}^{-2}/4\right)$
$\exp \left(-{\lambda }^2{r}_{12}^2\right)$	$\delta \left(t-\lambda \right)$
$r_{12}^{-1}\mathrm{erfc}\left(\lambda r_{12}\right)$	$\left(2/\sqrt{\pi }\right)\theta \left(t-\lambda \right)$
$r_{12}^{-1}\mathrm{erf}\left(\lambda r_{12}\right)$	$\left(2/\sqrt{\pi }\right)\left[1-\theta \left(t-\lambda \right)\right]$

Table 2. Table 2: Number of intermediates required to compute various integral classes for two-, three- and four-electron operators. The path generating the minimum number of intermediates is highlighted in bold. The number of terms in the RRs and the associated incremental center are also reported.

Integral	type	operator	path	number	centers	integral class
				of terms		$\left[p\dots p\right]$	$\left[d\dots d\right]$	$\left[f\dots f\right]$
two-electron		$f_{12}$	VV	(4,6)	(${𝐀}_2$,${𝐀}_1$)	4	13	25
			VT	(4,4)	(${𝐀}_2$,${𝐀}_1$)	7	19	37
three-electron	chain	$f_{12}{h}_{23}$	VVV	(6,7,12)	(${𝐀}_3$,${𝐀}_1$,${𝐀}_2$)	34	230	881
			VVV	(6,10,9)	(${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	32	209	778
			VVV	(8,8,9)	(${𝐀}_2$,${𝐀}_3$,${𝐀}_1$)	32	212	801
			VVT	(6,7,6)	(${𝐀}_3$,${𝐀}_1$,${𝐀}_2$)	38	246	873
			VVT	(6,10,6)	(${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	43	314	1,256
			VVT	(8,8,6)	(${𝐀}_2$,${𝐀}_3$,${𝐀}_1$)	40	260	923
	cyclic	$f_{12}{g}_{13}{h}_{23}$	VVV	(8,10,12)	(${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	52	469	2,216
			VVT	(8,10,6)	(${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	61	539	2,426
four-electron	chain	$f_{12}{h}_{23}{i}_{34}$	VVVV	(10,11,20,24)	(${𝐀}_4$,${𝐀}_1$,${𝐀}_3$,${𝐀}_2$)	465	13,781	150,961
			VVVV	(10,18,20,17)	(${𝐀}_4$,${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	436	12,535	133,891
			VVVV	(14,14,20,17)	(${𝐀}_3$,${𝐀}_4$,${𝐀}_2$,${𝐀}_1$)	433	12,704	138,913
			VVVV	(14,18,16,17)	(${𝐀}_3$,${𝐀}_2$,${𝐀}_4$,${𝐀}_1$)	435	12,863	141,679
			VVVT	(10,11,20,8)	(${𝐀}_4$,${𝐀}_1$,${𝐀}_3$,${𝐀}_2$)	532	16,295	181,178
			VVVT	(10,18,20,8)	(${𝐀}_4$,${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	560	17,029	188,242
			VVVT	(14,14,20,8)	(${𝐀}_3$,${𝐀}_4$,${𝐀}_2$,${𝐀}_1$)	559	17,487	199,050
			VVVT	(14,18,16,8)	(${𝐀}_3$,${𝐀}_2$,${𝐀}_4$,${𝐀}_1$)	543	16,612	185,869
	trident	$g_{13}{h}_{23}{i}_{34}$	VVVV	(10,12,14,28)	(${𝐀}_4$,${𝐀}_2$,${𝐀}_1$,${𝐀}_3$)	445	13,139	143,619
			VVVV	(10,12,24,18)	(${𝐀}_4$,${𝐀}_2$,${𝐀}_3$,${𝐀}_1$)	447	13,381	148,911
			VVVV	(10,20,16,18)	(${𝐀}_4$,${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	418	12,447	133,853
			VVVV	(16,14,16,18)	(${𝐀}_3$,${𝐀}_4$,${𝐀}_2$,${𝐀}_1$)	418	12,054	129,322
			VVVT	(10,12,14,8)	(${𝐀}_4$,${𝐀}_2$,${𝐀}_1$,${𝐀}_3$)	470	13,306	136,584
			VVVT	(10,12,24,8)	(${𝐀}_4$,${𝐀}_2$,${𝐀}_3$,${𝐀}_1$)	546	16,917	191,171
			VVVT	(10,20,16,8)	(${𝐀}_4$,${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	521	15,515	168,958
			VVVT	(16,14,16,8)	(${𝐀}_3$,${𝐀}_4$,${𝐀}_2$,${𝐀}_1$)	499	13,969	142,264

Equations

$$f_{12}=\exp (-\lambda r_{12}^2),$$

$$f_{12}=\exp (-\lambda r_{12}).$$

$$\hat{I}\approx \sum _{\mu }^{N_{\text{RI}}}|{\chi }_{\mu }⟩⟨{\chi }_{\mu }|.$$

$${\phi }_{\mathbf{a}}^{\mathbf{A}}(\mathbf{r})=(x-A_x{)}^{a_x}(y-A_y{)}^{a_y}(z-A_z{)}^{a_z}{e}^{-\alpha {|\mathbf{r}-\mathbf{A}|}^2}.$$

$${\psi }_{\mathbf{a}}^{\mathbf{A}}(\mathbf{r})=\sum _{k=1}^{K_A}{D}_{\mathbf{a}k}(x-A_x{)}^{a_x}(y-A_y{)}^{a_y}(z-A_z{)}^{a_z}{e}^{-{\alpha }_k{|\mathbf{r}-\mathbf{A}|}^2},$$

$$\begin{array}{rl}\text{innerproduct}{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4{\mathbf{b}}_1{\mathbf{b}}_2{\mathbf{b}}_3{\mathbf{b}}_4& \equiv \text{matrixelement}{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4{f}_{1234}{\mathbf{b}}_1{\mathbf{b}}_2{\mathbf{b}}_3{\mathbf{b}}_4\\ & =\iiint {\psi }_{{\mathbf{a}}_1}^{{\mathbf{A}}_1}({\mathbf{r}}_1){\psi }_{{\mathbf{a}}_2}^{{\mathbf{A}}_2}({\mathbf{r}}_2){\psi }_{{\mathbf{a}}_3}^{{\mathbf{A}}_3}({\mathbf{r}}_3){\psi }_{{\mathbf{a}}_4}^{{\mathbf{A}}_4}({\mathbf{r}}_4)f_{1234}{\psi }_{{\mathbf{b}}_1}^{{\mathbf{B}}_1}({\mathbf{r}}_1){\psi }_{{\mathbf{b}}_2}^{{\mathbf{B}}_2}({\mathbf{r}}_2){\psi }_{{\mathbf{b}}_3}^{{\mathbf{B}}_3}({\mathbf{r}}_3){\psi }_{{\mathbf{b}}_4}^{{\mathbf{B}}_4}({\mathbf{r}}_4)d{\mathbf{r}}_{1}d{\mathbf{r}}_{2}d{\mathbf{r}}_{3}d{\mathbf{r}}_4.\end{array}$$

$$\begin{array}{rl}[{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4|{\mathbf{b}}_1{\mathbf{b}}_2{\mathbf{b}}_3{\mathbf{b}}_4]& \equiv [{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4|f_{1234}|{\mathbf{b}}_1{\mathbf{b}}_2{\mathbf{b}}_3{\mathbf{b}}_4]\\ & =\iiint {\phi }_{{\mathbf{a}}_1}^{{\mathbf{A}}_1}({\mathbf{r}}_1){\phi }_{{\mathbf{a}}_2}^{{\mathbf{A}}_2}({\mathbf{r}}_2){\phi }_{{\mathbf{a}}_3}^{{\mathbf{A}}_3}({\mathbf{r}}_3){\phi }_{{\mathbf{a}}_4}^{{\mathbf{A}}_4}({\mathbf{r}}_4)f_{1234}{\phi }_{{\mathbf{b}}_1}^{{\mathbf{B}}_1}({\mathbf{r}}_1){\phi }_{{\mathbf{b}}_2}^{{\mathbf{B}}_2}({\mathbf{r}}_2){\phi }_{{\mathbf{b}}_3}^{{\mathbf{B}}_3}({\mathbf{r}}_3){\phi }_{{\mathbf{b}}_4}^{{\mathbf{B}}_4}({\mathbf{r}}_4)d{\mathbf{r}}_{1}d{\mathbf{r}}_{2}d{\mathbf{r}}_{3}d{\mathbf{r}}_4.\end{array}$$

$$[0]=S_1{S}_2{S}_3{S}_4\iiint {\phi }_0^{{\mathbf{Z}}_1}({\mathbf{r}}_1){\phi }_0^{{\mathbf{Z}}_2}({\mathbf{r}}_2){\phi }_0^{{\mathbf{Z}}_3}({\mathbf{r}}_3){\phi }_0^{{\mathbf{Z}}_4}({\mathbf{r}}_4)f_{1234}d{\mathbf{r}}_{1}d{\mathbf{r}}_{2}d{\mathbf{r}}_{3}d{\mathbf{r}}_4,$$

$${\zeta }_i$$

$$f_{12}={\int }_0^{\infty }F(t_{12})\exp (-t_{12}^2{r}_{12}^2)d{t}_{12},$$

$$[0]=S_1{S}_2{S}_3{S}_4\text{iiiint}F(t_{12})G(t_{13})H(t_{23})I(t_{34})w_0(\mathbf{t})d\mathbf{t},$$

$$w_0(\mathbf{t})=\text{quantity}(\frac{{\pi }^4}{D(\mathbf{t})}{)}^{3/2}\exp [-\frac{N(\mathbf{t})}{D(\mathbf{t})}].$$

$$s_1$$

$$s_4$$

$$s_7$$

$$D(\mathbf{t})$$

$$+{\zeta }_4\left({\zeta }_1+{\zeta }_2+{\zeta }_3\right)s_5+\left({\zeta }_1+{\zeta }_2\right)\left({\zeta }_3+{\zeta }_4\right)s_6+{\zeta }_2\left({\zeta }_1+{\zeta }_3+{\zeta }_4\right)s_7+{\zeta }_1\left({\zeta }_2+{\zeta }_3+{\zeta }_4\right)s_8$$

$$+\left({\zeta }_1+{\zeta }_2+{\zeta }_3+{\zeta }_4\right)s_9,$$

$$N(\mathbf{t})$$

$$+\left[\left({\zeta }_3+{\zeta }_4\right){\kappa }_{12}+\left({\zeta }_1+{\zeta }_2\right){\kappa }_{34}\right]s_6+{\zeta }_2\left({\kappa }_{13}+{\kappa }_{14}+{\kappa }_{34}\right)s_7+{\zeta }_1\left({\kappa }_{23}+{\kappa }_{24}+{\kappa }_{34}\right)s_8$$

$$+\left({\kappa }_{12}+{\kappa }_{13}+{\kappa }_{14}+{\kappa }_{23}+{\kappa }_{24}+{\kappa }_{34}\right)s_9,$$

$$[0{]}^{\mathbf{m}}=S_1{S}_2{S}_3{S}_4\text{iiiint}F(t_{12})G(t_{13})H(t_{23})I(t_{34})w_{\mathbf{m}}(\mathbf{t})d\mathbf{t},$$

$$w_{\mathbf{m}}(\mathbf{t})=w_0(\mathbf{t})\prod _{k=1}^{\dim \mathbf{m}}\text{quantity}[\frac{s_k}{D(\mathbf{t})}{]}^{m_k}.$$

$$\begin{split}[\bm{a}_{1}^{+}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\bm{m}}&=\mathbf{Z}_{1}\mathbf{A}_{1}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\bm{m}}-\zeta_{2}\zeta_{3}\zeta_{4}\mathbf{Z}_{12}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{1\}}-\zeta_{2}\zeta_{3}\zeta_{4}\mathbf{Z}_{13}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{2\}}\\ &-\zeta_{4}(\zeta_{2}\mathbf{Z}_{12}+\zeta_{3}\mathbf{Z}_{13})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{5\}}-\zeta_{2}(\zeta_{3}+\zeta_{4})\mathbf{Z}_{12}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{6\}}\\ &-\zeta_{2}(\zeta_{3}\mathbf{Z}_{13}+\zeta_{4}\mathbf{Z}_{14})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{7\}}-(\zeta_{2}\mathbf{Z}_{12}+\zeta_{3}\mathbf{Z}_{13}+\zeta_{4}\mathbf{Z}_{14})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{9\}}\\ &+\frac{\bm{a}_{1}}{2\zeta_{1}}\Big{\{}[\bm{a}_{1}^{-}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\bm{m}}-\zeta_{2}\zeta_{3}\zeta_{4}\mathbf{Z}_{12}[\bm{a}_{1}^{-}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{1\}}-\zeta_{2}\zeta_{3}\zeta_{4}\mathbf{Z}_{13}[\bm{a}_{1}^{-}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{2\}}\\ &-\zeta_{4}(\zeta_{2}+\zeta_{3})[\bm{a}_{1}^{-}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{5\}}-\zeta_{2}(\zeta_{3}+\zeta_{4})[\bm{a}_{1}^{-}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{6\}}\\ &-\zeta_{2}(\zeta_{3}+\zeta_{4})[\bm{a}_{1}^{-}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{7\}}-(\zeta_{2}+\zeta_{3}+\zeta_{4})[\bm{a}_{1}^{-}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{9\}}\Big{\}}\\ &+\frac{\bm{a}_{2}}{2}\Big{\{}\zeta_{3}\zeta_{4}[\bm{a}_{1}\bm{a}_{2}^{-}\bm{a}_{3}\bm{a}_{4}]^{\{1\}}+\zeta_{4}[\bm{a}_{1}\bm{a}_{2}^{-}\bm{a}_{3}\bm{a}_{4}]^{\{5\}}+(\zeta_{3}+\zeta_{4})[\bm{a}_{1}\bm{a}_{2}^{-}\bm{a}_{3}\bm{a}_{4}]^{\{6\}}+[\bm{a}_{1}\bm{a}_{2}^{-}\bm{a}_{3}\bm{a}_{4}]^{\{9\}}\Big{\}}\\ &+\frac{\bm{a}_{3}}{2}\Big{\{}\zeta_{2}\zeta_{4}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{2\}}+\zeta_{4}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{5\}}+\zeta_{2}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{7\}}+[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{9\}}\Big{\}}\\ &+\frac{\bm{a}_{4}}{2}\Big{\{}\zeta_{2}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}^{-}]^{\{7\}}+[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}^{-}]^{\{9\}}\Big{\}},\end{split}$$

$$\begin{split}[\bm{a}_{1}\bm{a}_{2}^{+}\bm{a}_{3}\bm{a}_{4}]^{\bm{m}}&=\mathbf{Z}_{2}\mathbf{A}_{2}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\bm{m}}+\zeta_{1}\zeta_{3}\zeta_{4}\mathbf{Z}_{12}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{1\}}-\zeta_{1}\zeta_{3}\zeta_{4}\mathbf{Z}_{23}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{3\}}\\ &+\zeta_{4}(\zeta_{1}\mathbf{Z}_{12}-\zeta_{3}\mathbf{Z}_{23})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{5\}}+\zeta_{1}(\zeta_{3}+\zeta_{4})\mathbf{Z}_{12}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{6\}}\\ &-\zeta_{1}(\zeta_{3}\mathbf{Z}_{23}+\zeta_{4}\mathbf{Z}_{24})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{8\}}+(\zeta_{1}\mathbf{Z}_{12}-\zeta_{3}\mathbf{Z}_{23}-\zeta_{4}\mathbf{Z}_{24})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{9\}}\\ &+\frac{\bm{a}_{2}}{2\zeta_{2}}\Big{\{}[\bm{a}_{1}\bm{a}_{2}^{-}\bm{a}_{3}\bm{a}_{4}]^{\bm{m}}-\zeta_{1}\zeta_{3}\zeta_{4}[\bm{a}_{1}\bm{a}_{2}^{-}\bm{a}_{3}\bm{a}_{4}]^{\{1\}}-\zeta_{1}\zeta_{3}\zeta_{4}[\bm{a}_{1}\bm{a}_{2}^{-}\bm{a}_{3}\bm{a}_{4}]^{\{3\}}\\ &-\zeta_{4}(\zeta_{1}-\zeta_{3})[\bm{a}_{1}\bm{a}_{2}^{-}\bm{a}_{3}\bm{a}_{4}]^{\{5\}}+\zeta_{1}(\zeta_{3}+\zeta_{4})[\bm{a}_{1}\bm{a}_{2}^{-}\bm{a}_{3}\bm{a}_{4}]^{\{6\}}\\ &-\zeta_{1}(\zeta_{3}+\zeta_{4})[\bm{a}_{1}\bm{a}_{2}^{-}\bm{a}_{3}\bm{a}_{4}]^{\{8\}}-(\zeta_{1}+\zeta_{3}+\zeta_{4})[\bm{a}_{1}\bm{a}_{2}^{-}\bm{a}_{3}\bm{a}_{4}]^{\{9\}}\Big{\}}\\ &+\frac{\bm{a}_{1}}{2}\Big{\{}\zeta_{3}\zeta_{4}[\bm{a}_{1}^{-}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{1\}}+\zeta_{4}[\bm{a}_{1}^{-}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{5\}}+(\zeta_{3}+\zeta_{4})[\bm{a}_{1}^{-}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{6\}}+[\bm{a}_{1}^{-}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{9\}}\Big{\}}\\ &+\frac{\bm{a}_{3}}{2}\Big{\{}\zeta_{1}\zeta_{4}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{3\}}+\zeta_{4}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{5\}}+\zeta_{1}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{8\}}+[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{9\}}\Big{\}}\\ &+\frac{\bm{a}_{4}}{2}\Big{\{}\zeta_{1}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}^{-}]^{\{8\}}+[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}^{-}]^{\{9\}}\Big{\}},\end{split}$$

$$\begin{split}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{+}\bm{a}_{4}]^{\bm{m}}&=\mathbf{Z}_{3}\mathbf{A}_{3}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\bm{m}}+\zeta_{1}\zeta_{2}\zeta_{4}\mathbf{Z}_{13}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{2\}}+\zeta_{1}\zeta_{2}\zeta_{4}\mathbf{Z}_{23}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{3\}}\\ &+\zeta_{1}\zeta_{2}\zeta_{4}\mathbf{Z}_{34}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{4\}}+\zeta_{4}(\zeta_{1}\mathbf{Z}_{13}+\zeta_{2}\mathbf{Z}_{23})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{5\}}-\zeta_{4}(\zeta_{1}+\zeta_{2})\mathbf{Z}_{34}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{6\}}\\ &+\zeta_{2}(\zeta_{1}\mathbf{Z}_{13}-\zeta_{4}\mathbf{Z}_{34})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{7\}}+\zeta_{1}(\zeta_{2}\mathbf{Z}_{23}-\zeta_{4}\mathbf{Z}_{34})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{8\}}\\ &+(\zeta_{1}\mathbf{Z}_{13}+\zeta_{2}\mathbf{Z}_{23}-\zeta_{4}\mathbf{Z}_{34})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{9\}}\\ &+\frac{\bm{a}_{3}}{2\zeta_{3}}\Big{\{}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\bm{m}}-\zeta_{1}\zeta_{2}\zeta_{4}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{2\}}-\zeta_{1}\zeta_{2}\zeta_{4}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{3\}}-\zeta_{1}\zeta_{2}\zeta_{4}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{4\}}\\ &-\zeta_{4}(\zeta_{1}+\zeta_{2})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{5\}}-\zeta_{4}(\zeta_{1}+\zeta_{2})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{6\}}-\zeta_{2}(\zeta_{1}+\zeta_{4})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{7\}}\\ &-\zeta_{1}(\zeta_{2}+\zeta_{4})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{8\}}-(\zeta_{1}+\zeta_{2}+\zeta_{4})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{9\}}\Big{\}}\\ &+\frac{\bm{a}_{1}}{2}\Big{\{}\zeta_{2}\zeta_{4}[\bm{a}_{1}^{-}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{2\}}+\zeta_{4}[\bm{a}_{1}^{-}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{5\}}+\zeta_{2}[\bm{a}_{1}^{-}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{7\}}+[\bm{a}_{1}^{-}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{9\}}\Big{\}}\\ &+\frac{\bm{a}_{2}}{2}\Big{\{}\zeta_{1}\zeta_{4}[\bm{a}_{1}\bm{a}_{2}^{-}\bm{a}_{3}\bm{a}_{4}]^{\{3\}}+\zeta_{4}[\bm{a}_{1}\bm{a}_{2}^{-}\bm{a}_{3}\bm{a}_{4}]^{\{5\}}+\zeta_{1}[\bm{a}_{1}\bm{a}_{2}^{-}\bm{a}_{3}\bm{a}_{4}]^{\{8\}}+[\bm{a}_{1}\bm{a}_{2}^{-}\bm{a}_{3}\bm{a}_{4}]^{\{9\}}\Big{\}}\\ &+\frac{\bm{a}_{4}}{2}\Big{\{}\zeta_{1}\zeta_{2}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}^{-}]^{\{4\}}+(\zeta_{1}+\zeta_{2})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}^{-}]^{\{6\}}+\zeta_{2}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}^{-}]^{\{7\}}+\zeta_{1}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}^{-}]^{\{8\}}\\ &+[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}^{-}]^{\{9\}}\Big{\}}\end{split}$$

$$\begin{split}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}^{+}]^{\bm{m}}&=\mathbf{Z}_{4}\mathbf{A}_{4}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\bm{m}}+\zeta_{1}\zeta_{2}\zeta_{3}\mathbf{Z}_{34}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{4\}}+\zeta_{3}(\zeta_{1}+\zeta_{2})\mathbf{Z}_{34}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{6\}}\\ &+\zeta_{2}(\zeta_{1}\mathbf{Z}_{14}+\zeta_{3}\mathbf{Z}_{34})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{7\}}+\zeta_{1}(\zeta_{2}\mathbf{Z}_{24}+\zeta_{3}\mathbf{Z}_{34})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{8\}}\\ &+(\zeta_{1}\mathbf{Z}_{14}+\zeta_{2}\mathbf{Z}_{24}+\zeta_{3}\mathbf{Z}_{34})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{9\}}\\ &+\frac{\bm{a}_{4}}{2\zeta_{4}}\Big{\{}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}^{-}]^{\bm{m}}-\zeta_{1}\zeta_{2}\zeta_{3}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}^{-}]^{\{4\}}-\zeta_{3}(\zeta_{1}+\zeta_{2})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}^{-}]^{\{6\}}\\ &-\zeta_{2}(\zeta_{1}+\zeta_{3})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}^{-}]^{\{7\}}-\zeta_{1}(\zeta_{2}+\zeta_{3})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}^{-}]^{\{8\}}-(\zeta_{1}+\zeta_{2}+\zeta_{3})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}^{-}]^{\{9\}}\Big{\}}\\ &+\frac{\bm{a}_{1}}{2}\Big{\{}\zeta_{2}[\bm{a}_{1}^{-}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{7\}}+[\bm{a}_{1}^{-}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}]^{\{9\}}\Big{\}}+\frac{\bm{a}_{2}}{2}\Big{\{}\zeta_{1}[\bm{a}_{1}\bm{a}_{2}^{-}\bm{a}_{3}\bm{a}_{4}]^{\{8\}}+[\bm{a}_{1}\bm{a}_{2}^{-}\bm{a}_{3}\bm{a}_{4}]^{\{9\}}\Big{\}}\\ &+\frac{\bm{a}_{3}}{2}\Big{\{}\zeta_{1}\zeta_{2}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{4\}}+(\zeta_{1}+\zeta_{2})[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{6\}}+\zeta_{2}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{7\}}+\zeta_{1}[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{8\}}\\ &+[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}^{-}\bm{a}_{4}]^{\{9\}}\Big{\}}\end{split}$$

$$\begin{array}{rl}[{\mathbf{a}}_1^{+}{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4]& =\frac{{\mathbf{a}}_1}{2{\zeta }_1}[{\mathbf{a}}_1^{-}{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4]+\frac{{\mathbf{a}}_2}{2{\zeta }_1}[{\mathbf{a}}_1{\mathbf{a}}_2^{-}{\mathbf{a}}_3{\mathbf{a}}_4]+\frac{{\mathbf{a}}_3}{2{\zeta }_1}[{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3^{-}{\mathbf{a}}_4]+\frac{{\mathbf{a}}_4}{2{\zeta }_1}[{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4^{-}]\\ & -\frac{{\zeta }_2}{{\zeta }_1}[{\mathbf{a}}_1{\mathbf{a}}_2^{+}{\mathbf{a}}_3{\mathbf{a}}_4]-\frac{{\zeta }_3}{{\zeta }_1}[{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3^{+}{\mathbf{a}}_4]-\frac{{\zeta }_4}{{\zeta }_1}[{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4^{+}]\\ & -\frac{{\beta }_1{\mathbf{A}}_1{\mathbf{B}}_1+{\beta }_2{\mathbf{A}}_2{\mathbf{B}}_2+{\beta }_3{\mathbf{A}}_3{\mathbf{B}}_3+{\beta }_4{\mathbf{A}}_4{\mathbf{B}}_4}{{\zeta }_1}[{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4].\end{array}$$

$$\text{innerproduct}{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4{\mathbf{b}}_4^{+}=\text{innerproduct}{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4^{+}{\mathbf{b}}_4+{\mathbf{A}}_4{\mathbf{B}}_4\text{innerproduct}{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4{\mathbf{b}}_4,$$

$$\text{innerproduct}{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4{\mathbf{b}}_3^{+}{\mathbf{b}}_4=\text{innerproduct}{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3^{+}{\mathbf{a}}_4{\mathbf{b}}_3{\mathbf{b}}_4+{\mathbf{A}}_3{\mathbf{B}}_3\text{innerproduct}{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4{\mathbf{b}}_3{\mathbf{b}}_4,$$

$$\text{innerproduct}{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4{\mathbf{b}}_2^{+}{\mathbf{b}}_3{\mathbf{b}}_4=\text{innerproduct}{\mathbf{a}}_1{\mathbf{a}}_2^{+}{\mathbf{a}}_3{\mathbf{a}}_4{\mathbf{b}}_2{\mathbf{b}}_3{\mathbf{b}}_4+{\mathbf{A}}_2{\mathbf{B}}_2\text{innerproduct}{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4{\mathbf{b}}_2{\mathbf{b}}_3{\mathbf{b}}_4,$$

$$\text{innerproduct}{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4{\mathbf{b}}_1^{+}{\mathbf{b}}_2{\mathbf{b}}_3{\mathbf{b}}_4=\text{innerproduct}{\mathbf{a}}_1^{+}{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4{\mathbf{b}}_1{\mathbf{b}}_2{\mathbf{b}}_3{\mathbf{b}}_4+{\mathbf{A}}_1{\mathbf{B}}_1\text{innerproduct}{\mathbf{a}}_1{\mathbf{a}}_2{\mathbf{a}}_3{\mathbf{a}}_4{\mathbf{b}}_1{\mathbf{b}}_2{\mathbf{b}}_3{\mathbf{b}}_4.$$

Full text

three-electron integral, four-electron integral, many-electron integral, Gaussian basis function, explicitly-correlated method, F12 method, recurrence relation

I

Introduction

In 1985, starting from the Hylleraas functional Hylleraas (1928, 1929); Sinagoǧlu (1962) and using the interelectronic distance $r_{12}=\absolutevalue{\bm{r}_{1}-\bm{r}_{2}}$ as a correlation factor, Kutzelnigg derived the seminal form of the MP2-R12 equations. Kutzelnigg (1985) This explicitly-correlated method (more formally stated together with Klopper in 1987 Klopper and Kutzelnigg (1987)) was later extended to higher level of theory and more accurate correlation factors $f_{12}=f(r_{12})$, Noga and Kutzelnigg (1994); Klopper and Kutzelnigg (1991); Klopper, Rohse, and Kutzelnigg (1991); Termath, Klopper, and Kutzelnigg (1991); Persson and Taylor (1996, 1997); May and Manby (2004); Ten-no (2004a); Tew and Klopper (2005); May et al. (2005) such as Gaussian geminals

[TABLE]

or Slater geminals

[TABLE]

The resulting “F12 methods” achieve chemical accuracy for small organic molecules with relatively small Gaussian basis sets Klopper et al. (2006); Hattig et al. (2012); Kong, Bischo, and Valeev (2012); Ten-no and Noga (2012); Ten-no (2012); Gruneis et al. (2017) and are quickly becoming the first-choice method for high accuracy. Ten-no and Noga (2012); Ten-no (2012)

However, the inclusion of the correlation factor $f_{12}$ dredged up an old problem: in addition to two-electron integrals (traditional ones and new ones), three-electron integrals over $f_{13}f_{23}$, $r_{12}^{-1}f_{13}$ and $r_{12}^{-1}f_{13}f_{23}$, as well as four-electron integrals over $r_{12}^{-1}f_{14}f_{23}$, $r_{12}^{-1}f_{13}f_{34}$ and $r_{12}^{-1}f_{13}f_{14}$ arise. Except when one uses Gaussian geminals, Boys (1960); Singer (1960) these integrals are not known analytically and, at that time, the only way to evaluate them would have been via expensive Gauss-Legendre quadratures. Preiskorn and Zurawski (1985); Clementi (1989) Additionally, citing Kutzelnigg and Klopper, Kutzelnigg and Klopper (1991) “even if fast procedures for the evaluation of these integrals were available, one would have to face the problem of the large number of these integrals; while that of two-electron integrals is $\sim N^{4}$, there are $\sim N^{6}$ three-electron and $\sim N^{8}$ four-electron integrals. The storing and manipulating of these integrals could be handled only for extremely small basis sets.”

Undoubtedly, in the late 80’s, the two-electron integrals technology was still in development. McMurchie and Davidson (1978); Pople and Hehre (1978); King and Dupuis (1976); Dupuis, Rys, and King (1976); Rys, Dupuis, and King (1983); Obara and Saika (1986); Head-Gordon and Pople (1988); Gill (1994) Nowadays, though still challenging, these integrals could be computed much more effectively via judicious recursive schemes, designing the quadrature only to the fundamental integrals. Barca, Loos, and Gill (2016) Another important remark is that the actual number of significant (i.e. greater than a given threshold) three- and four-electron integrals in a large system, is, at worst, $\order{N^{3}}$ or $\order{N^{4}}$. These kinds of scaling are achievable, for example, by exploiting robust density fitting Womack and Manby (2014) or upper bound-based screening methods. Barca, Loos, and Gill (tion)

Nevertheless, the success of the R12 method was due to the decision of avoiding three- and four-electron integrals entirely through the insertion of the resolution of the identity (RI) Kutzelnigg and Klopper (1991); Hattig et al. (2012); Werner, Adler, and Manby (2007)

[TABLE]

In this way, three- and four-electron integrals are approximated as linear combinations of products of more conventional two-electron integrals. Of course, the accuracy of the RI approximation (3) relies entirely on the assumption that the auxiliary basis set is sufficiently large ($N_{\text{RI}}\gg N$). Thus, in the general context of explicitly correlated methods, it is licit to ask: what is the suitable method to evaluate three- and four-electron integrals?

It not clear to us that such method is RI. In fact, eschewing the RI approximation would offer at least two advantages: i) smaller one-electron basis as the larger auxiliary basis set would not be required anymore; ii) the three- and four-electron integrals would be computed exactly. Moreover, one could avoid the commutator rearrangements involved in the computation of integrals over the kinetic energy operator. Rohse, Klopper, and Kutzelnigg (1993)

In a recent paper, Barca, Loos, and Gill (2016) we reported recurrence relations (RRs) to compute three-electron integrals over Gaussian basis functions for multiplicative general operators of the form $f_{12}g_{13}h_{23}$. Here, we generalize our previous study to four-electron integrals.

The present paper is organized as follows. In Sec. II, we introduce notations and we define various key quantities. Section III explains how to calculate fundamental integrals required to start the recursive scheme. In Sec. IV, we reports vertical, transfer and horizontal RRs for four-electron integrals. Finally, in Sec. V, we propose a recursive scheme based on these RRs to calculate classes of four-electron integrals. Atomic units are used throughout.

II

Four-electron integrals

A primitive Gaussian-type function (PGF) is specified by an orbital exponent $\alpha$, a center $\mathbf{A}=(A_{x},A_{y},A_{z})$, and angular momentum $\bm{a}=(a_{x},a_{y},a_{z})$:

[TABLE]

A contracted Gaussian-type function (CGF) is defined as a normalized sum of PGFs

[TABLE]

where $K_{A}$ is the degree of contraction and the $D_{\bm{a}k}$ are contraction coefficients. Throughout this paper, we use physicists notations, and we write the integral over a four-electron operator $f_{1234}$ of CGFs as

[TABLE]

Additionally, square-bracketed integrals denote integrals over PGFs:

[TABLE]

The fundamental integral (i.e. the integral in which all eight basis functions are $s$-type PGFs) is defined as $[\mathbf{0}]\equiv[\mathbf{0}\mathbf{0}\mathbf{0}\mathbf{0}|\mathbf{0}\mathbf{0}\mathbf{0}\mathbf{0}]$ with $\mathbf{0}=(0,0,0)$. The Gaussian product rule reduces it from eight to four centers:

[TABLE]

where

[TABLE]

with $\mathbf{A}_{i}\mathbf{B}_{i}=\mathbf{A}_{i}-\mathbf{B}_{i}$. For conciseness, we will adopt a notation in which missing indices represent $s$-type Gaussians. For example, $[\bm{a}_{2}\bm{a}_{3}]$ is a shorthand for $[\mathbf{0}\bm{a}_{2}\bm{a}_{3}\mathbf{0}|\mathbf{0}\mathbf{0}\mathbf{0}\mathbf{0}]$. We will also use unbold indices, e.g. $[a_{1}a_{2}a_{3}a_{4}|b_{1}b_{2}b_{3}b_{4}]$ to indicate a complete class of integrals from a shell-octet.

II.1

Four-electron operators

In the present study, we are particularly interested in the four-electron operators $g_{13}h_{23}i_{34}$ (trident) and $f_{12}h_{23}i_{34}$ (four-electron chain or 4-chain) because they can be required in explicitly-correlated methods such as F12 methods. Ten-no and Noga (2012); Ten-no (2012); Kong, Bischo, and Valeev (2012); Hattig et al. (2012) Explicitly-correlated calculations may also require three-electron integrals over the (3-chain) $f_{12}h_{23}$ and (cyclic) $f_{12}g_{13}h_{23}$ operators, as well as two-electron integrals over $f_{12}$. However, we will eschew the study of the two-electron integrals here as they have been extensively studied in the past 25 years. Kutzelnigg and Klopper (1991); Klopper and Rohse (1992); Persson and Taylor (1997); Klopper and Samson (2002); Manby (2003); Werner, Manby, and Knowles (2003); Klopper (2004); Ten-no (2004a, b); May et al. (2005); Manby et al. (2006); Ten-no (2007); Komornicki and King (2011); Reine, Helgaker, and Lind (2012) Note that the nuclear attraction integrals can be easily obtained by taking the large-exponent limit of a $s$-type shell-pair. We refer the interested reader to Refs. Gill and Pople, 1991; Gill, 1994 for more details about the computation of nuclear attraction integrals.

The structure of these operators is illustrated in Fig. 1, where we have adopted a diagrammatic representation. Starting with the “pacman” operator $f_{12}g_{13}h_{23}i_{34}$, we are going to show that one can easily derive all the RRs required to compute two-, three- and four-electron integrals following simple rules. Therefore, in the following, we will focus our analysis on this master “pacman” operator.

III

Fundamental integrals

The first step required to compute integrals of arbitrary angular momentum is the computation of the (momentumless) fundamental integrals $[\mathbf{0}]$. These are derived starting from Eq. (8) using the Gaussian integral representation of each two-electron operator. For instance, we have

[TABLE]

where $F(t_{12})$ is a Gaussian kernel. Table 1 contains kernels $F(t)$ for a variety of important two-electron operators $f_{12}$. From the formulas in Table 1, one can also easily deduce the kernels for related functions, such as $f_{12}^{2}$, $f_{12}/r_{12}$, and $\nabla^{2}f_{12}$, which are of interest in explicitly-correlated methods. More general kernels can be found in Ref. Barca, Loos, and Gill, 2016.

Next, the integration over $\bm{r}_{1}$, $\bm{r}_{2}$, $\bm{r}_{3}$ and $\bm{r}_{4}$ can be carried out, yielding

[TABLE]

where $\bm{t}=(t_{12},t_{13},t_{23},t_{34})$ and

[TABLE]

Defining the following polynomials

[TABLE]

we have

[TABLE]

where $\kappa_{ij}=\zeta_{i}\zeta_{j}\,\absolutevalue{\mathbf{Z}_{i}-\mathbf{Z}_{j}}^{2}$.

Following Obara and Saika, Obara and Saika (1986, 1988) vertical RRs (VRRs) are obtained by differentiation of Eq. (11) with respect to the center coordinates. Therefore, one can show that (11) has to be generalized to the following form:

[TABLE]

where $\bm{m}=(m_{1},m_{2},m_{3},m_{4},m_{5},m_{6},m_{7},m_{8},m_{9})$ is, for the pacman operator, a nine-dimensional auxiliary index, and

[TABLE]

As reported in Fig. 1, while the fundamental integrals of the pacman operator contains 9 auxiliary indices (i.e. $\dim\bm{m}=9$), the two interesting four-electron operators (trident and chain) contains 7. The cyclic and chain three-electron operators have only four and three, respectively, while the two-electron chain operator has a single $m$ component (i.e. $\dim\bm{m}=1$). These numbers are drastically reduced if one uses Gaussian geminals due to their factorization properties. Barca and Loos (tion)

IV

Recurrence relations

In this Section, we report vertical, transfer and horizontal RRs for the computation of four-electron integrals of arbitrary angular momentum. In particular, we refer the interested readers to the appendix of Ref. Barca, Loos, and Gill, 2016 for more details about how to derive these VRRs à la Ahlrichs. Ahlrichs (2006)

IV.1

Vertical recurrence relations

To build angular momentum over center $\mathbf{A}_{1}$, we have derived the following 24-term VRR:

[TABLE]

where the superscript $+$ or $-$ denotes an increment or decrement of one unit of Cartesian angular momentum. (Thus, $\bm{a}^{\pm}$ is analogous to $\bm{a}\pm\bm{1}_{i}$ in the notation of Obara and Saika.) The value in the curly superscript indicates which component of the auxiliary index vector $\bm{m}$ is incremented.

Because Eq. (17) builds angular momemtum over $\mathbf{A}_{1}$ and all four bra centers have non-zero angular momentum, we will call this expression VRR${}_{4}^{\mathbf{A}_{1}}$. The VRRs used to obtain $[\bm{a}_{1}^{+}\bm{a}_{2}\bm{a}_{3}]^{\bm{m}}$, $[\bm{a}_{1}^{+}\bm{a}_{2}]^{\bm{m}}$ and $[\bm{a}_{1}^{+}]^{\bm{m}}$ can be easily derived from Eq. (17) by setting successively $\bm{a}_{4}=\mathbf{0}$, $\bm{a}_{3}=\mathbf{0}$ and $\bm{a}_{2}=\mathbf{0}$. These are respectively named VRR${}_{3}^{\mathbf{A}_{1}}$, VRR${}_{2}^{\mathbf{A}_{1}}$ , VRR${}_{1}^{\mathbf{A}_{1}}$.

One can easily derive VRR4 for the trident and 4-chain operators following the simple rules given in Fig. 1. We obtain the VRRs for the trident operator by removing the terms $\{1\}$ and $\{6\}$. Similarly, 4-chain VRRs are obtained by removing the terms $\{2\}$ and $\{7\}$ . This yields a 18- and 17-term VRR${}_{4}^{\mathbf{A}_{1}}$ for the trident and 4-chain operators, respectively.

VRR${}_{4}^{\mathbf{A}_{2}}$, VRR${}_{4}^{\mathbf{A}_{3}}$ and VRR${}_{4}^{\mathbf{A}_{4}}$ are used to build angular momentum over $\mathbf{A}_{2}$, $\mathbf{A}_{3}$ and $\mathbf{A}_{4}$, respectively. Their expressions are

[TABLE]

[TABLE]

[TABLE]

Again, the corresponding expressions for VRR1, VRR2 , VRR3 can be easily derived from Eqs. (18), (19) and (20). The number of terms for each of these VRRs is reported in Fig. 2 for the 3-chain $f_{12}h_{23}$ (top left), cyclic $f_{12}g_{13}h_{23}$ (top right), 4-chain $f_{12}h_{23}i_{34}$ (bottom left) and trident $g_{13}h_{23}i_{34}$ (bottom right) operators.

IV.2

Transfer recurrence relations

TRRs redistribute angular momentum between centers hosting to different electrons. Using translational invariance, one can derive

[TABLE]

IV.3

Horizontal recurrence relations

The so-called HRRs enable to shift momentum between centers over the same electronic coordinate:

[TABLE]

Note that HRRs can be applied to contracted integrals because they are independent of the contraction coefficients and exponents.

V

Algorithm

In this Section, we describe a recursive scheme for the computation of three- and four-electron integrals based on a late contraction scheme à la Head-Gordon-Pople (HGP). Head-Gordon and Pople (1988) The general skeleton of the algorithm is shown in Fig. 3 for two representative examples: the 3-chain operator $f_{12}g_{23}$ (left) and the 4-chain operator $f_{12}h_{23}i_{34}$ (right). First, let us focus on the 3-chain operator.

To compute a class of three-electron integrals $\innerproduct{a_{1}a_{2}a_{3}}{b_{1}b_{2}b_{3}}$, starting from the fundamental integrals $[000]^{\bm{m}}$, we first build up angular momentum over center $\mathbf{A}_{3}$ with the 6-term VRR1 to obtain $[00a_{3}]$. Then, we use the 10-term VRR2 over $\mathbf{A}_{2}$ to obtain $[0a_{2}a_{3}]$. Finally, we build up momentum over the last bra center $\mathbf{A}_{1}$ using the 9-term VRR3 to get $[a_{1}a_{2}a_{3}]$.

Note that, in our previous paper, Barca, Loos, and Gill (2016) we claimed that it would be computationally cheaper to use the 6-term TRR instead of VRR3 because the number of terms in the TRR is much smaller than in VRR3. However, we have found here that the number of intermediates (i.e. the number of pre-computed classes needed to calculate a given class) required by the paths involving the TRR is much larger (see Table 2). This can be easily understood noticing that, to increase the momentum by one unit on the last center, one must increase the momentum by the same amount on all the other centers (as evidenced by the second term in the right-hand side of (21)). Hence, the TRR is computationally expensive for three- and four-electron integrals due to the large number of centers.

As illustrated in the top left graph of Fig. 2, other paths, corresponding to different VRRs, are possible. However, we have found that they do generate a larger number of intermediates, as reported in Table 2.

Similarly to the 3-chain operator, for the 4-chain operator, we get $[a_{1}a_{2}a_{3}a_{4}]$ by successively building up momentum over $\mathbf{A}_{4}$, $\mathbf{A}_{3}$, $\mathbf{A}_{2}$ and $\mathbf{A}_{1}$. The number of intermediates required by the other paths are gathered in Table 2. Again, the paths involving the 8-term TRR (reported in Eq. (21)) are much more expensive.

The last two steps of the algorithm are common to the three- and four-electron integral schemes. Following the HGP algorithm, Head-Gordon and Pople (1988) we contract the integrals: $[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}|\mathbf{0}\mathbf{0}\mathbf{0}\mathbf{0}]$ to form $\innerproduct{\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}\bm{a}_{4}}{\mathbf{0}\mathbf{0}\mathbf{0}\mathbf{0}}$ in the four-electron case, or $[\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}|\mathbf{0}\mathbf{0}\mathbf{0}]$ to form $\innerproduct{\bm{a}_{1}\bm{a}_{2}\bm{a}_{3}}{\mathbf{0}\mathbf{0}\mathbf{0}}$ in the three-electron case. More details about the contraction step can be found in Ref. Gill, 1994. The final step of the algorithm shifts momentum to the ket centers from the bra centers with the help of the 2-term HRRs reported in Sec. IV.3.

VI Concluding remarks

In this study, we have reported recurrence relations (RRs) for the efficient and accurate computation of four-electron integrals over Gaussian basis functions and a general class of multiplicative four-electron operators of the form $f_{12}g_{13}h_{23}i_{34}$. Starting from this master operator, one can easily derive the RRs for various operators arising in explicitly-correlated methods following simple diagrammatic rules (see Fig. 1).

Here, we have derived three types of RRs: i) starting from the fundamental integrals, vertical RRs (VRRs) allow to increase the angular momentum over the bra centers; ii) the transfer RR (TRR) redistributes angular momentum between centers hosting different electrons, and can be used instead of the VRR on the last bra center; iii) the horizontal RRs (HRRs) enable to shift momentum from the bra to the ket centers corresponding to the same electronic coordinate. Importantly, HRRs can be applied to contracted integrals.

Finally, after carefully studying the different paths one can follow to build up angular momentum (see Fig. 2), we have proposed a late-contraction recursive scheme which minimizes the number of intermediates to be computed (see Fig. 3). We believe our approach represents a major step towards a fast and accurate computational scheme for three- and four-electron integrals within explicitly-correlated methods. It also paves the way to contraction-effective methods for these types of integrals. Barca and Gill (2016) In particular, an early contraction scheme would have significant computational benefits.

Acknowledgements.

P.F.L. thanks the NCI National Facility for generous grants of supercomputer time, and the Australian Research Council for a Discovery Project grant (DP140104071).