Three- and four-electron integrals involving Gaussian geminals: fundamental integrals, upper bounds and recurrence relations

TL;DR

This paper develops methods for efficiently calculating three- and four-electron integrals involving Gaussian geminals, including fundamental integrals, recurrence relations, and upper bounds, to improve explicitly-correlated F12 quantum chemistry calculations.

Contribution

It introduces a comprehensive approach with fundamental integrals, recurrence relations, and upper bounds for three- and four-electron integrals involving Gaussian geminals, optimizing computational efficiency.

Findings

Derived straightforward methods for fundamental integrals.

Established recurrence relations for angular momentum buildup.

Provided upper bounds to identify significant integrals.

Abstract

We report the three main ingredients to calculate three- and four-electron integrals over Gaussian basis functions involving Gaussian geminal operators: fundamental integrals, upper bounds, and recurrence relations. In particular, we consider the three- and four-electron integrals that may arise in explicitly-correlated F12 methods. A straightforward method to obtain the fundamental integrals is given. We derive vertical, transfer and horizontal recurrence relations to build up angular momentum over the centers. Strong, simple and scaling-consistent upper bounds are also reported. This latest ingredient allows to compute only the $\order{N^2}$ significant three- and four-electron integrals, avoiding the computation of the very large number of negligible integrals.

Tables

Table 1. Table 1: Primitive bound sets [ I m ] delimited-[] subscript 𝐼 𝑚 [I_{m}] for three- and four-electron integrals and number N sig subscript 𝑁 sig N_{\text{sig}} of significant shell- m 𝑚 m tuplets in a large system with N 𝑁 N PGFs.

Integral	type	operator	shell-pair level ($m=2$)		shell-quartet level ($m=4$)
			$\left[I_2\right]/\left[\hat{1}\right]$	$N_{\text{sig}}$	$\left[I_4\right]$	$N_{\text{sig}}$
three-electron	chain	$C_{12}{G}_{13}$	$\{\left[\hat{1}\right]{\left[1\right]}^{12,13},{\left[\hat{1}\right]}^{13}{\left[1\right]}^{12},{\left[\hat{1}\right]}^{12}{\left[1\right]}^{13}\}$	$\left(N\right)$	$\{\left[\hat{1}\right]{\left[13\right]}^{12,13},{\left[\hat{1}\right]}^{12}{\left[13\right]}^{13}\}$	$\left(N\right)$
		$G_{13}{G}_{23}$	$\{{\left[\hat{1}\right]}^{23}{\left[1\right]}^{13},{\left[\hat{1}\right]}^{13}{\left[1\right]}^{23},\left[\hat{1}\right]{\left[1\right]}^{13,23}\}$	$\left(N\right)$	$\{{\left[\hat{1}\right]}^{23}{\left[13\right]}^{13},\left[\hat{1}\right]{\left[13\right]}^{13,23}\}$	$\left(N\right)$
	cyclic	$C_{12}{G}_{13}{G}_{23}$	$\{{\left[\hat{1}\right]}^{23}{\left[1\right]}^{12,13},{\left[\hat{1}\right]}^{13}{\left[1\right]}^{12,23},{\left[\hat{1}\right]}^{12}{\left[1\right]}^{13,23}\}$	$\left(N\right)$	$\{\left[\hat{1}\right]{\left[13\right]}^{12,13},{\left[\hat{1}\right]}^{12}{\left[13\right]}^{13}\}$	$\left(N\right)$
four-electron	chain	$C_{12}{G}_{14}{G}_{23}$	$\left[\hat{1}\right]\{{\left[\hat{1}\right]}^{23}{\left[1\right]}^{12,14},{\left[\hat{1}\right]}^{12,14}{\left[1\right]}^{23}\}$	$\left(N\right)$	$\left[\hat{1}\right]\{\left[\hat{1}\right]{\left[14\right]}^{12,14},{\left[\hat{1}\right]}^{12,23}{\left[14\right]}^{14},{\left[\hat{1}\right]}^{12}{\left[14\right]}^{14}\}$	$\left(N\right)$
		$C_{12}{G}_{13}{G}_{34}$	$\left[\hat{1}\right]\{{\left[\hat{1}\right]}^{34}{\left[1\right]}^{12,13},{\left[1\right]}^{34}{\left[\hat{1}\right]}^{12,13}\}$	$\left(N\right)$	$\{\left[\hat{1}\right]{\left[\hat{1}\right]}^{34}{\left[13\right]}^{12,13},{\left[\hat{1}\right]}^{12}{\left[\hat{1}\right]}^{34}{\left[13\right]}^{13}\}$	$\left(N\right)$
	trident	$C_{12}{G}_{13}{G}_{14}$	$\left[\hat{1}\right]\{{\left[1\right]}^{12,13,14}\left[\hat{1}\right],\left[1\right]{\left[\hat{1}\right]}^{12,13,14}\}$	$\left(N\right)$	$\{\left[\hat{1}\right]{\left[\hat{1}\right]}^{14}{\left[13\right]}^{12,13},{\left[\hat{1}\right]}^{12}{\left[\hat{1}\right]}^{14}{\left[13\right]}^{13}\}$	$\left(N\right)$
Integral	type	operator	shell-sextet level ($m=6$)		shell-octet level ($m=8$)
			$\left[I_6\right]$	$N_{\text{sig}}$	$\left[I_8\right]$	$N_{\text{sig}}$
three-electron	chain	$C_{12}{G}_{13}$	$\{{\left[13\right]}^{12,13}\left[2\right],{\left[13\right]}^{13}{\left[2\right]}^{12}\}$	$\left(N^2\right)$	—	—
		$G_{13}{G}_{23}$	$\left\{{\left[123\right]}^{13,23}\right\}$	$\left(N\right)$	—	—
	cyclic	$C_{12}{G}_{13}{G}_{23}$	$\left\{{\left[123\right]}^{13,23}\right\}$	$\left(N\right)$	—	—
four-electron	chain	$C_{12}{G}_{14}{G}_{23}$	$\left[\hat{1}\right]\{{\left[14\right]}^{12,14}{\left[2\right]}^{23},{\left[14\right]}^{14}{\left[2\right]}^{12,23}\}$	$\left(N^2\right)$	$\{{\left[14\right]}^{12,14}{\left[23\right]}^{23},{\left[14\right]}^{14}{\left[23\right]}^{12,23}\}$	$\left(N^2\right)$
		$C_{12}{G}_{13}{G}_{34}$	$\{\left[\hat{1}\right]{\left[134\right]}^{12,13,34},{\left[\hat{1}\right]}^{12}{\left[134\right]}^{13,34}\}$	$\left(N\right)$	$\{{\left[134\right]}^{12,13,34}\left[2\right],{\left[134\right]}^{13,34}{\left[2\right]}^{12}\}$	$\left(N^2\right)$
	trident	$C_{12}{G}_{13}{G}_{14}$	$\{\left[\hat{1}\right]{\left[134\right]}^{12,13,14},{\left[\hat{1}\right]}^{12}{\left[134\right]}^{13,14}\}$	$\left(N\right)$	$\{{\left[134\right]}^{12,13,14}\left[2\right],{\left[134\right]}^{13,14}{\left[2\right]}^{12}\}$	$\left(N^2\right)$

Table 2. Table 2: Bounds factors for three- and four-electron integrals. ζ ̊ i = ( 1 − σ i ) ​ ζ i subscript ̊ 𝜁 𝑖 1 subscript 𝜎 𝑖 subscript 𝜁 𝑖 \mathring{\zeta}_{i}=(1-\sigma_{i})\zeta_{i} , and h i subscript ℎ 𝑖 h_{i} , δ ̊ u subscript ̊ 𝛿 𝑢 \mathring{\delta}_{u} and Y ̊ ̊ 𝑌 \mathring{Y} are given by Eqs. ( 18a ), ( 18b ) and ( 18c ).

Bound factor	expression
${\left[1\right]}^{13,14}$	$h_1{\left(\frac{\pi }{{\mathring{\zeta }}_1+{\lambda }_{13}+{\lambda }_{14}}\right)}^{3/2}$
${\left[1\right]}^{12,13,14}$	$h_1\frac{2\pi }{{\mathring{\zeta }}_1+{\lambda }_{13}+{\lambda }_{14}}$
${\left[13\right]}^{13}$	$h_1{h}_3{\left(\frac{{\pi }^2}{{\mathring{\delta }}_0^{13}}\right)}^{3/2}\exp \left(-{\mathring{Y}}^{13}\right)$
${\left[13\right]}^{12,13}$	$\frac{2}{\sqrt{\pi }}\sqrt{\frac{{\mathring{\delta }}_0^{13}}{{\mathring{\delta }}_1^{13}-{\mathring{\delta }}_0^{13}}}{\left[13\right]}^{13}$
${\left[134\right]}^{13,14}$	$h_1{h}_3{h}_4{\left(\frac{{\pi }^3}{{\mathring{\delta }}_0^{13,14}}\right)}^{3/2}\exp \left(-{\mathring{Y}}^{13,14}\right)$
${\left[134\right]}^{12,13,14}$	$\frac{2}{\sqrt{\pi }}\sqrt{\frac{{\mathring{\delta }}_0^{13,14}}{{\mathring{\delta }}_1^{13,14}-{\mathring{\delta }}_0^{13,14}}}{\left[134\right]}^{13,14}$

Table 3. Table 3: 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]$	$\left[g\dots g\right]$
two-electron	chain	$G_{12}$	VV	(2,3)	(${𝐀}_2$,${𝐀}_1$)	3	6	10	15
			VT	(2,4)	(${𝐀}_2$,${𝐀}_1$)	4	9	16	25
		$C_{12}$	VV	(4,6)	(${𝐀}_2$,${𝐀}_1$)	4	13	25	48
			VT	(4,4)	(${𝐀}_2$,${𝐀}_1$)	7	19	37	61
three-electron	chain	$G_{13}{G}_{23}$	VVV	(2,3,4)	(${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	5	13	26	45
			VVT	(2,3,6)	(${𝐀}_3$,${𝐀}_1$,${𝐀}_2$)	8	25	56	105
		$C_{12}{G}_{23}$	VVV	(4,5,7)	(${𝐀}_3$,${𝐀}_1$,${𝐀}_2$)	11	39	96	195
			VVV	(4,6,6)	(${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	10	39	96	196
			VVT	(4,5,6)	(${𝐀}_3$,${𝐀}_1$,${𝐀}_2$)	16	66	173	359
			VVT	(4,6,6)	(${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	15	65	171	357
	cyclic	$C_{12}{G}_{13}{G}_{23}$	VVV	(4,6,8)	(${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	12	46	119	250
			VVT	(4,6,6)	(${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	16	66	173	359
four-electron	chain	$C_{12}{G}_{14}{G}_{23}$	VVVV	(4,5,7,8)	(${𝐀}_4$,${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	21	108	344	847
			VVVV	(4,6,6,8)	(${𝐀}_4$,${𝐀}_1$,${𝐀}_3$,${𝐀}_2$)	19	88	260	607
			VVVT	(4,5,7,8)	(${𝐀}_4$,${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	33	208	736	1,926
			VVVT	(4,6,6,8)	(${𝐀}_4$,${𝐀}_1$,${𝐀}_3$,${𝐀}_2$)	33	204	716	1,866
		$C_{12}{G}_{13}{G}_{34}$	VVVV	(4,6,6,9)	(${𝐀}_4$,${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	22	113	360	888
			VVVV	(4,6,8,7)	(${𝐀}_4$,${𝐀}_3$,${𝐀}_1$,${𝐀}_2$)	20	98	302	726
			VVVT	(4,6,6,8)	(${𝐀}_4$,${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	33	204	716	1,866
			VVVT	(4,6,8,8)	(${𝐀}_4$,${𝐀}_3$,${𝐀}_1$,${𝐀}_2$)	34	214	756	1,976
	trident	$C_{12}{G}_{13}{G}_{14}$	VVVV	(4,6,6,9)	(${𝐀}_4$,${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	22	113	360	888
			VVVV	(4,6,8,7)	(${𝐀}_4$,${𝐀}_3$,${𝐀}_1$,${𝐀}_2$)	20	98	302	726
			VVVT	(4,6,6,8)	(${𝐀}_4$,${𝐀}_3$,${𝐀}_2$,${𝐀}_1$)	33	204	716	1,866
			VVVT	(4,6,8,8)	(${𝐀}_4$,${𝐀}_3$,${𝐀}_1$,${𝐀}_2$)	34	214	756	1,976

Equations

$$S_{12}=\exp (-{\lambda }_{12}{r}_{12}).$$

$$G_{12}=\exp (-{\lambda }_{12}{r}_{12}^2),$$

$$r_{12}\approx \sum _k{a}_k\text{quantity}[1-\exp (-{\lambda }_k{r}_{12}^2)],$$

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

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

$$|\mathbf{a}\mathbf{b}⟩\equiv {\psi }_{\mathbf{a}}^{\mathbf{A}}(\mathbf{r}){\psi }_{\mathbf{b}}^{\mathbf{B}}(\mathbf{r})=\sum _{i=1}^{K_a}\sum _{j=1}^{K_b}|\mathbf{a}\mathbf{b}{]}_{ij}$$

$$\begin{array}{rl}\text{innerproduct}{\mathbf{a}}_1\cdots {\mathbf{a}}_n{\mathbf{b}}_1\cdots {\mathbf{b}}_n& \equiv \text{matrixelement}{\mathbf{a}}_1\cdots {\mathbf{a}}_n{f}_{1\cdots n}{\mathbf{b}}_1\cdots {\mathbf{b}}_n\\ & =\text{idotsint}{\psi }_{{\mathbf{a}}_1}^{{\mathbf{A}}_1}({\mathbf{r}}_1)\cdots {\psi }_{{\mathbf{a}}_n}^{{\mathbf{A}}_n}({\mathbf{r}}_n)f_{1\cdots n}{\psi }_{{\mathbf{b}}_1}^{{\mathbf{B}}_1}({\mathbf{r}}_1)\cdots {\psi }_{{\mathbf{b}}_n}^{{\mathbf{B}}_n}({\mathbf{r}}_n)d{\mathbf{r}}_1\cdots d{\mathbf{r}}_n.\end{array}$$

$$[{\mathbf{a}}_1\cdots {\mathbf{a}}_n|{\mathbf{b}}_1\cdots {\mathbf{b}}_n]=\text{idotsint}{\phi }_{{\mathbf{a}}_1}^{{\mathbf{A}}_1}({\mathbf{r}}_1)\cdots {\phi }_{{\mathbf{a}}_n}^{{\mathbf{A}}_n}({\mathbf{r}}_n)f_{1\cdots n}{\phi }_{{\mathbf{b}}_1}^{{\mathbf{B}}_1}({\mathbf{r}}_1)\cdots {\phi }_{{\mathbf{b}}_n}^{{\mathbf{B}}_n}({\mathbf{r}}_n)d{\mathbf{r}}_1\cdots d{\mathbf{r}}_n.$$

$$[0]=\text{quantity}(\prod _{i=1}^n{S}_i)\text{idotsint}{\phi }_0^{{\mathbf{Z}}_1}({\mathbf{r}}_1)\cdots {\phi }_0^{{\mathbf{Z}}_n}({\mathbf{r}}_n)f_{1\cdots n}d{\mathbf{r}}_1\cdots d{\mathbf{r}}_n,$$

$${\zeta }_i$$

$${\mathbf{Z}}_i$$

$$S_i$$

$$C_{12}=\frac{2}{\sqrt{\pi }}{\int }_0^{\infty }\exp (-u^2{r}_{12}^2)du.$$

$$[0{]}^m=\frac{2}{\sqrt{\pi }}[0{]}_G\sqrt{\frac{{\delta }_0}{{\delta }_1-{\delta }_0}}\text{quantity}(\frac{{\delta }_1}{{\delta }_1-{\delta }_0}{)}^m{F}_m\text{quantity}[\frac{{\delta }_1\text{quantity}(Y_1-Y_0)}{{\delta }_1-{\delta }_0}],$$

$$[0{]}_G=\text{quantity}(\prod _{i=1}^4{S}_i)\text{quantity}(\frac{{\pi }^4}{{\delta }_0}{)}^{3/2}\exp (-Y_0)$$

$${\mathbf{\delta }}_u=\mathbf{\zeta }+{\mathbf{\lambda }}_u=\mathbf{\zeta }+\mathbf{G}+u^2\mathbf{C},$$

$$\mathbf{\zeta }$$

$$\mathbf{G}$$

$${\mathbf{\Delta }}_u$$

$${\delta }_u$$

$$[B_m]=\{\begin{array}{ll}\max [I_m],& m=2,\\ \min [I_m],& m>2,\end{array}$$

$$h=\text{absolutevalue}{D}_a{D}_b{\gamma }_a{\gamma }_b\sqrt{\frac{a^a{b}^b}{e^{a+b}}}\frac{\text{quantity}(4\alpha \beta {)}^{3/4}}{{\sigma }^{\frac{a+b}{2}}}{e}^{-\frac{(1-\sigma )\text{absolutevalue}{\mathbf{AB}}^2}{{\alpha }^{-1}+{\beta }^{-1}}},$$

$${\mathring{\delta }}_u^{13,14,23,34}=\det ({\mathring{\mathbf{\delta }}}_u^{13,14,23,34}),$$

$${\mathring{Y}}^{13,14,23,34}=\text{Tr}({\mathring{\mathbf{\Delta }}}_0^{13,14,23,34}\cdot {\mathbf{Y}}^2),$$

$${\gamma }_a=\prod _{i=0}^2\text{quantity}[\Gamma \text{quantity}(\lfloor \frac{a+i}{3}\rfloor +\frac{1}{2}){]}^{-1/2},$$

$${\mathring{\mathbf{\delta }}}_u^{13,14,23,34}=\mathring{\mathbf{\zeta }}+{\mathbf{\lambda }}_u=\left(\begin{array}{cccc}{\mathring{\zeta }}_1& 0& 0& 0\\ 0& {\mathring{\zeta }}_2& 0& 0\\ 0& 0& {\mathring{\zeta }}_3& 0\\ 0& 0& 0& {\mathring{\zeta }}_4\end{array}\right)+{\mathbf{\lambda }}_u,$$

$${\mathring{\mathbf{\Delta }}}_u^{13,14,23,34}=\mathring{\mathbf{\zeta }}\cdot {\left({\mathring{\mathbf{\delta }}}_u^{13,14,23,34}\right)}^{-1}\cdot \mathring{\mathbf{\zeta }},$$

$$\text{partialderivative}[1{]}^{13,14,23,34}{\sigma }_1=0.$$

$$\text{expectationvalue}{1}^{13,14}=\sum _{ij}^K[1{]}_{ij}^{13,14},$$

$$\text{expectationvalue}{13}^{12,13,23}$$

$$\text{expectationvalue}{123}^{12,13,23}$$

$$\text{expectationvalue}{134}^{12,13,14,34}$$

$${\check{Y}}^{13,14,23,34}=\text{Tr}({\check{\mathbf{\Delta }}}_0^{13,14,23,34}\cdot {\check{\mathbf{Y}}}^2),$$

$${\check{\mathbf{\delta }}}_u^{13,14,23,34}=\check{\mathbf{\zeta }}+{\mathbf{\lambda }}_u=\left(\begin{array}{cccc}{\check{\zeta }}_1& 0& 0& 0\\ 0& {\check{\zeta }}_2& 0& 0\\ 0& 0& {\check{\zeta }}_3& 0\\ 0& 0& 0& {\check{\zeta }}_4\end{array}\right)+{\mathbf{\lambda }}_u,$$

$${\check{\mathbf{\Delta }}}_u^{13,14,23,34}=\check{\mathbf{\zeta }}\cdot {\left({\check{\mathbf{\delta }}}_u^{13,14,23,34}\right)}^{-1}\cdot \check{\mathbf{\zeta }},$$

$${\check{\mathbf{Y}}}^k=\left(\begin{array}{cccc}0& {\check{Y}}_{12}^k& {\check{Y}}_{13}^k& {\check{Y}}_{14}^k\\ 0& 0& 0& {\check{Y}}_{24}^k\\ 0& 0& 0& {\check{Y}}_{34}^k\\ 0& 0& 0& 0\end{array}\right),$$

Full text

I

Introduction

It is well known that highly-accurate wave functions require the fulfilment (or near-fulfilment) of the electron-electron cusp conditions. Kato (1951, 1957); Pack and Byers Brown (1966); Morgan III and Kutzelnigg (1993); Myers et al. (1991); Tew (2008); Loos and Gill (2009, 2010, 2012); Kurokawa, Nakashima, and Nakatsuji (2013, 2014); Gruneis et al. (2017) For correlated wave functions expanded in terms of products of one-electron Gaussian basis functions, the energy converges as $\order{L^{-3}}$, where $L$ is the maximum angular momentum of the basis set. Kutzelnigg (1985) This slow convergence can be tracked down to the inadequacy of these products to properly model the Coulomb correlation hole. Hattig et al. (2012); Kong, Bischo, and Valeev (2012)

In the late 20’s, Hylleraas solved this issue for the helium atom by introducing explicitly the interelectronic distance $r_{12}=\absolutevalue{\bm{r}_{1}-\bm{r}_{2}}$ as an additional two-electron basis function. Hylleraas (1928, 1929) As Kutzelnigg later showed, this leads to a prominent improvement of the energy convergence from $\order{L^{-3}}$ to $\order{L^{-7}}$. Kutzelnigg (1985)

Around the same time, Slater, while studying the Rydberg series of helium, Slater (1928) suggested a new correlation factor known nowadays as a Slater geminal:

[TABLE]

Unfortunately, the increase in mathematical complexity brought by $r_{12}$ or $S_{12}$ has been found to be rapidly computationally overwhelming.

In 1960, Boys Boys (1960) and Singer Singer (1960) independently proposed to use the Gaussian geminal (GG) correlation factor

[TABLE]

as “there are explicit formulas for all of the necessary many-dimensional integrals”. Boys (1960) Interestingly, in the same article, a visionary Boys argued that, even if GGs do not fulfil the electron-electron cusp conditions, they could be used to fit $S_{12}$.

During the following years, variational calculations involving GGs flourished, giving birth to various methods, such as the exponentially-correlated Gaussian method. Rychlewski and Komasa (2003); Bubin, Cafiero, and Adamowicz (2005); Bubin, Stanke, and Adamowicz (2009); Szalewicz and Jeziorski (2010) However, this method was restricted to fairly small systems as it requires the optimization of a large number of non-linear parameters. In the 70’s, the first MP2 calculations including GGs appeared thanks to the work by Pan and King, Pan and King (1970, 1972) Adamowicz and Sadlej, Adamowicz and Sadlej (1977a, b, 1978) and later Szalewicz et al. Szalewicz et al. (1982, 1983) Even if these methods represented a substantial step forward in terms of computational burden, they still require the optimization of a large number of non-linear parameters.

In 1985, Kutzelnigg derived a first form of the MP2-R12 equations using $r_{12}$ as a correlation factor. Kutzelnigg (1985) Kutzelnigg’s idea, which was more formally described together with Klopper in 1987,Klopper and Kutzelnigg (1987) dredged up an old problem: in addition to two-electron integrals (traditional ones and new ones), three-electron and four-electron integrals were required. At that time, the only way to evaluate them would have been via an expensive one- or two-dimensional Gauss-Legendre quadrature. Preiskorn and Zurawski (1985); Clementi (1989) Nevertheless, the success of the R12 method was triggered by the decision to avoid the computation of these three- and four-electron integrals through the use of the resolution of the identity (RI) approximation. Kutzelnigg and Klopper (1991); Hattig et al. (2012); Werner, Adler, and Manby (2007) In this way, three- and four-electron integrals are approximated as linear combinations of products of two-electron integrals. Several key developments and improvements of the original MP2-R12 approach have been proposed in the last decade. Banchorz et al. (2011); Hattig et al. (2012); Werner, Adler, and Manby (2007); Kong, Bischo, and Valeev (2012) Of course, the accuracy of the RI approximation relies entirely on the assumption that the auxiliary basis set is sufficiently large, i.e. $N_{\text{RI}}\gg N$, where $N$ and $N_{\text{RI}}$ are the number of basis functions in the primary and auxiliary basis sets, respectively.

In 1996, Persson and Taylor killed two birds with one stone. Using a pre-optimized GG expansion fitting $r_{12}$

[TABLE]

they avoided the non-linear optimization, and eschewed the RI approximation thanks to the analytical integrability of three- and four-electron integrals over GGs. Persson and Taylor (1996) They were able to show that a six- or ten-term fit introduces a $0.5$ mE${}_{\text{h}}$ or $20$ $\mu$E${}_{\text{h}}$ error, respectively. Persson and Taylor (1996) Unfortunately, further improvements were unsuccessful due to the failure of $r_{12}$ in modelling the correct behaviour of the wave function for intermediate and large $r_{12}$. Tew and Klopper (2005, 2006) In fact, Ten-no showed that $S_{12}$ is near-optimal at describing the correlation hole, and that a 10-term GG fit of $S_{12}$ yields very similar results. This suggests that, albeit not catching the cusp per se, the Coulomb correlation hole can be accurately represented by GGs. Ten-no (2004a, 2007); May and Manby (2004); May et al. (2005)

Methods for evaluating many-electron integrals involving GGs have already been developed. As mentioned previously, Persson and Taylor Persson and Taylor (1997) derived recurrence relations based on Hermite Gaussians, analogously to the work of McMurchie and Davidson for two-electron integrals. McMurchie and Davidson (1978) These recurrence relations were implemented by Dahle. Dahle (2004); Dahle et al. (2007, 2008) Saito and Suzuki Saito and Suzuki (2001) also proposed an approach based on the work by Obara and Saika. Obara and Saika (1986, 1988) More recently, a general formulation using Rys polynomials King and Dupuis (1976); Dupuis, Rys, and King (1976); Rys, Dupuis, and King (1983) was published by Komornicki and King. Komornicki and King (2011) Even if limited to the three-center case, it is worth mentioning that May has also developed recurrence relations for two types of three-electron integrals. May (2006) These recurrence relations were implemented by Womack using automatically-generated code. Womack (2015) Recently, we have developed recurrence relations for three- and four-electron integrals for generic correlation factors. Barca, Loos, and Gill (2016); AQC17

A major limitation of all these approaches is that they do not include any integral screening.111Komornicki and King mentioned the crucial importance of an effective integral screening in Ref. Komornicki and King, 2011. Indeed, a remarkable consequence of the short-range nature of the Slater and Gaussian correlation factors is that, even if formally scaling as $\order{N^{6}}$ and $\order{N^{8}}$, there are only $\order{N^{2}}$ significant (i.e. greater than a given threshold) three- and four-electron integrals in a large system. Barca, Loos, and Gill (2016, tion) Therefore, it is paramount to devise rigorous upper bounds to avoid computing the large number of negligible integrals.

The present manuscript is organized as follows. In Sec. II, we discuss Gaussian basis functions, many-electron integrals and the structure of the three- and four-electron operators considered here. The next three sections contain the main ingredients for the efficient computation of three- and four-electron integrals involving GGs: i) fundamental integrals (FIs) in Sec. III, ii) upper bounds (UBs) in Sec. IV, and iii) recurrence relations (RRs) in Sec. V. In Sec. VI, we give an overall view of our algorithm which is an extension of the late-contraction path of PRISM (see Refs. Head-Gordon and Pople, 1988; Gill, 1994 and references therein). Note that the RRs developed in this study differ from the ones reported in our previous studies (Refs. Barca, Loos, and Gill, 2016; Barca and Loos, m ph) as they are specifically tailored for the unique factorization properties brought by the association of Gaussian basis functions and GGs. Atomic units are used throughout.

II

Generalities

II.1

Gaussian functions

A primitive Gaussian 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 function (CGF) is defined as a sum of PGFs

[TABLE]

where $K_{a}$ is the degree of contraction and the $D_{ak}$ are contraction coefficients. A CGF-pair

[TABLE]

is a two-fold sum of PGF-pairs $|\bm{a}\bm{b}]=\varphi_{\bm{a}}^{\mathbf{A}}(\bm{r})\varphi_{\bm{b}}^{\mathbf{B}}(\bm{r})$.

A primitive shell $|a]$ is a set of PGFs sharing the same total angular momentum $a$, exponent $\alpha$ and center $\mathbf{A}$. Similarly, a contracted shell $\ket{a}$ is a set of CGFs sharing the same PGFs and total angular momentum. A contracted shell-pair is the set of CGF-pairs obtained by the tensor product $\ket{ab}=\ket{a}\otimes\ket{b}$. Similarly, a primitive shell-pair $|ab]=|a]\otimes|b]$ is the set of PGF-pairs. Finally, primitive and contracted shell-quartets, -sextets and -octets are obtained in an analogous way. For example, $|a_{1}b_{1}a_{2}b_{2}]=|a_{1}b_{1}]\otimes|a_{2}b_{2}]$ and $\ket{a_{1}a_{2}b_{1}b_{2}}=\ket{a_{1}b_{1}}\otimes\ket{a_{2}b_{2}}$. Note that $|1]$ is a set of three $p$-type PGFs, a $|11]\equiv|pp]$ shell-pair is a set of nine PGF-pairs, and a $|2222]\equiv|dddd]$ shell-quartet is a set of $1,296$ PGF-quartets.

II.2

Many-electron integrals

Throughout this paper, we use physicists notations, and we write the integral over a $n$-electron operator $f_{1\cdots n}$ of CGFs as

[TABLE]

Additionally, square-bracketed integrals denote integrals over PGFs:

[TABLE]

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

[TABLE]

where

[TABLE]

and $\mathbf{A}_{i}\mathbf{B}_{i}=\mathbf{A}_{i}-\mathbf{B}_{i}$. We also define the quantity $\mathbf{Y}_{ij}=\mathbf{Z}_{i}-\mathbf{Z}_{j}$ which will be used later on.

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.3

Three- and four-electron operators

In this study, we are particularly interested in the “master” four-electron operator $C_{12}G_{13}G_{14}G_{23}G_{34}$ (where $C_{12}=r_{12}^{-1}$ is the Coulomb operator) because the three types of three-electron integrals and the three types of four-electron integrals that can be required in F12 calculations can be easily generated from it (see Fig. 1). These three types of three-electron integrals are composed by a single type of integrals over the cyclic operator $C_{12}G_{13}G_{23}$, and two types of integrals over the three-electron chain (or 3-chain) operators $C_{12}G_{23}$ and $G_{13}G_{23}$. F12 calculations may also require three types of four-electron integrals: two types of integrals over the 4-chain operators $C_{12}G_{14}G_{23}$ and $C_{12}G_{13}G_{34}$, as well as one type over the trident operator $C_{12}G_{13}G_{14}$. Explicitly-correlated methods also requires two-electron integrals. However, their computation has been thoroughly studied in the literature. 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); Ten-no and Noga (2012); Ten-no (2012); Reine, Helgaker, and Lind (2012); Kong, Bischo, and Valeev (2012); Hattig et al. (2012) Similarly, the nuclear attraction integrals can be easily obtained by taking the large-exponent limit of a $s$-type shell-pair.

Starting with the “master” operator $C_{12}G_{13}G_{14}G_{23}G_{34}$, we will show that one can easily obtain all the FIs as well as the RRs required to compute three- and four-electron integrals within F12 calculations. This is illustrated in Fig. 1 where we have used a diagrammatic representation of the operators. The number $N_{\text{sig}}$ of significant integrals in a large system with $N$ CGFs is also reported.

III

Fundamental integrals

Following Persson and Taylor, Persson and Taylor (1996) the $[\mathbf{0}]^{\bm{m}}$ are derived starting from the momentumless integral (9) using the following Gaussian integral representation for the Coulomb operator

[TABLE]

After a lengthy derivation, one can show that the closed-form expression of the FIs is

[TABLE]

where $m$ is an auxiliary index, $F_{m}(t)$ is the generalized Boys function, and

[TABLE]

is the FI of the “pure” GG operator $G_{13}G_{14}G_{23}G_{34}$ from which one can easily get the FI of the 3-chain operator $G_{13}G_{23}$ by setting $\lambda_{14}=\lambda_{34}=0$. While the FIs involving a Coulomb operator contain an auxiliary index $m$, the FIs over “pure” GG operators (like $G_{13}G_{23}$) do not, thanks to the factorization properties of GGs. Barca and Gill (2016)

The various quantities required to compute (12) are

[TABLE]

where

[TABLE]

and

[TABLE]

The generalized Boys function $F_{m}(t)$ in Eq. (12) can be computed efficiently using well-established algorithms. Gill, Johnson, and Pople (1991); Ishida (1996); Weiss and Ochsenfeld (2015)

IV

Upper bounds

In this section, we report UBs for primitive and contracted three- and four-electron integrals. Our UBs are required to be simple (i.e. significantly computationally cheaper that the true integral), strong (i.e. as close as possible to the true integral in the threshold region $10^{-14}-10^{-8}$), and scaling-consistent (i.e. the number of significant integrals $N_{\text{sig}}=\mathcal{O}(N_{\text{UB}})$, where $N_{\text{UB}}$ is the number of integrals estimated by the UB). We refer the interested reader to Refs. Gill, Johnson, and Pople, 1994; Barca and Gill, 2016 for additional information about UBs. A detailed study of these UBs (as well as their overall performances) will be reported in a forthcoming paper. Barca, Loos, and Gill (tion)

Our screening algorithms are based on primitive, $[B_{m}]$, and contracted, $\expectationvalue{B_{m}}$, shell-$m$tuplet bounds. These are based on shell-$m$tuplet information only: shell-pair ($m=2$), shell-quartet ($m=4$), shell-sextet ($m=6$) and shell-octet ($m=8$). Thus, for each category of three- and four-electron integrals, we will report from shell-pair to shell-sextet (or shell-octet) bounds.

Figure 2 is a schematic representation of the overall screening scheme for contracted four-electron integrals. First, we use a primitive shell-pair bound $[B_{2}]$ to create a list of significant primitive shell-pairs. For a given contracted shell-pair, if at least one of its primitive shell-pairs has survived, a contracted shell-pair bound $\expectationvalue{B_{2}}$ is used to decide whether or not this contracted shell-pair is worth keeping. The second step consists in using a shell-quartet bound $\expectationvalue{B_{4}}$ to create a list of significant contracted shell-quartets by pairing the contracted shell-pairs with themselves. Then, we combine the significant shell-quartets and shell-pairs, and a shell-sextet bound $\expectationvalue{B_{6}}$ identifies the significant contracted shell-sextets. Finally, the shell-sextets are paired with the shell-pairs. If the resulting shell-octet quantity is found to be significant, the contracted integral class $\innerproduct{a_{1}a_{2}a_{3}a_{4}}{b_{1}b_{2}b_{3}b_{4}}$ must be computed via RRs, as discussed in the next section. The number of significant shell-$m$tuplets generated at each step is given in Table 1. As one can see, the size of any shell-$m$tuplet list is, at worst, quadratic in a large system.

During the shell-pair screening, either a contracted or a primitive path is followed depending on the degree of contraction of the integral class $K_{\text{tot}}=\prod_{i=1}^{n}K_{a_{i}}K_{b_{i}}$. If $K_{\text{tot}}>1$, the contracted path is enforced, otherwise the primitive path is followed. This enables to adopt the more effective primitive bounds for primitive integral classes which are usually associated with medium and high angular momentum PGFs and, therefore, are more expensive to evaluate via RRs. The scheme for primitive four-electron integrals differs only by the use of primitive bounds instead of contracted ones. The three-electron integrals screening scheme can be easily deduced from Fig. 2.

Note that we bound an entire class of integrals with a single UB. This is a particularly desirable feature, especially when dealing with three- or four-electron integrals where the size of a class can be extremely large. For example, the simple $[ppp|ppp]$ and $[pppp|pppp]$ classes are made of 729 and 4,096 integrals!

IV.1

Primitive bounds

In this section we present UBs for primitive three- and four-electron integrals. Without loss of generality, we assume that the geminal operators are ordered by decreasing exponent, i.e. $\lambda_{13}\leq\lambda_{14}\leq\lambda_{23}\leq\lambda_{34}$.

All the required primitive bounds have the form

[TABLE]

where $m$ is the shell multiplicity. The bound sets $[I_{m}]$ are reported in Table 1. They also require the bound factors in Table 2, which are easily computed with the following quantities:

[TABLE]

where

[TABLE]

$\Gamma(x)$ and $\lfloor x\rfloor$ are the Gamma and floor functions respectively, Olver et al. (2010) and $\mathring{\zeta}_{i}=(1-\sigma_{i})\zeta_{i}$. We point out that bound factors $[1]^{12}$ or $[1]^{12,13}$ reported in Table 2 can be obtained from $[1]^{12,13,14}$ by setting $\lambda_{13}=\lambda_{14}=0$ or $\lambda_{14}=0$, respectively. The parameter $\sigma_{1}$ is obtained by solving the quadratic equation

[TABLE]

Factors of the kind $[\widehat{1}]^{{13},{14},{23},{34}}$ are also required for the bounds in Table 1. They are defined as the largest factor $[1]^{{13},{14},{23},{34}}$ within a given system and basis set, and can be pre-computed and stored with the remaining basis set information.

IV.2

Contracted bounds

Contracted integral bounds are straightforward variations of primitive ones. While contracting at the shell-pair level ($m=2$) only requires $\order{K^{2}}$ computational work, contracting at the shell-quartet, -sextet or -octet level would require $\order{K^{4}}$, $\order{K^{6}}$ or $\order{K^{8}}$ work, respectively.

Therefore, as sketched in Fig. 2, we use a primitive bound for a first screening of the shell-pairs, then contracted bounds are used to screen shell-pairs, -quartets, -sextets and -octets. Considering $K$-fold CGFs, the contraction step never exceeds $\mathcal{O}(K^{2})$ computational cost. Bound factors such as

[TABLE]

and its maximum within the basis set $\langle\widehat{1}\rangle^{13,14}$, are computed within the shell-pair loop. In Eq. (21), $i$ and $j$ refer to the PGFs $|a_{1}]_{i}$ and $|b_{1}]_{j}$ in the contracted shells $\ket{a_{1}}$ and $\ket{b_{1}}$, respectively (see Eq. (6)).

The expressions of the contracted bounds are identical to the primitive bounds, with the only exception that the contracted factors $\expectationvalue{13}^{12,13,23}$, $\expectationvalue{123}^{12,13,23}$ and $\expectationvalue{134}^{12,13,14,34}$ are bound by

[TABLE]

where

[TABLE]

can be evaluated with the following expressions

[TABLE]

where $\check{\zeta}_{i}$ is the smallest effective exponent $\mathring{\zeta}_{i}$ in the contracted shell-pair $\ket{a_{i}b_{i}}$, and

[TABLE]

is the distance between two spheres of diameters $\mathbf{A}_{i}\mathbf{B}_{i}$ and $\mathbf{A}_{j}\mathbf{B}_{j}$ (where $\mathbf{A}_{i}\mathbf{B}_{i}^{+}=\mathbf{A}_{i}+\mathbf{B}_{i}$).

V

Recurrence relations

V.1

Vertical recurrence relations

Following Obara and Saika, Obara and Saika (1986, 1988) vertical RRs (VRRs) are obtained by differentiation of Eq. (12) with respect to the centers coordinates. Ahlrichs (2006); Barca, Loos, and Gill (2016) For the integrals considered in this study, one can show that

[TABLE]

where $\delta_{ij}$ is the Kronecker delta, Olver et al. (2010)

[TABLE]

and

[TABLE]

with

[TABLE]

One can easily derive VRRs for other three- and four-electron operators following the simple rules given in Fig. 1. The number of terms for each of these VRRs is reported in Table 3 for various two-, three- and four-electron operators.

Note that for a pure GG operator, we have $m=0$ and $Y_{1}=Y_{0}$. Therefore, Eq. (26) reduces to a simpler expression:

[TABLE]

V.2

Transfer recurrence relations

Transfer RRs (TRRs) redistribute angular momentum between centers referring to different electrons. Barca, Loos, and Gill (2016) Using the translational invariance, one can derive

[TABLE]

Note that Eq. (31) can only be used to build up angular momentum on the last center. Moreover, to increase the momentum by one unit on this 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 (31)). Therefore, the TRR is computationally expensive for three- and four-electron integrals due to the large number of centers (see below). As mentioned by Ahlrichs, Ahlrichs (2006) the TRR can be beneficial for very high angular momentum two-electron integral classes.

V.3

Horizontal recurrence relations

The so-called horizontal RRs (HRRs) enable to shift momentum between centers over the same electronic coordinate: Barca, Loos, and Gill (2016)

[TABLE]

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

VI

Algorithm

In this Section, we propose a recursive algorithm for the computation of a class of three- or four-electron integrals of arbitrary angular momentum. The present recursive algorithm is based on a late-contraction scheme inspired by the Head-Gordon-Pople algorithm Head-Gordon and Pople (1988) following a BOVVVVCCCCHHHH path. The general skeleton of the algorithm is shown in Fig. 3 for the trident operator $C_{12}G_{13}G_{14}$. We will use this example to illustrate each step.

Based on the shell data, the first step of the algorithm (step B) is to decide whether or not a given class of integrals is significant or negligible. If the integral class is found to be significant by the screening algorithm presented in Sec. IV and depicted in Fig. 2, an initial set of FIs is computed (step O) via the formulae gathered in Sec. III.

Starting with these FIs, angular momentum is then built up over the different bra centers $\mathbf{A}_{1}$, $\mathbf{A}_{2}$, $\mathbf{A}_{3}$ and $\mathbf{A}_{4}$ using the VRRs derived in Sec. V.1. To minimize the computational cost, one has to think carefully how to perform this step. Indeed, the cost depends on the order in which this increase in angular momentum is performed. This is illustrated in Fig. 4, where we have represented the various possible pathways for the 3-chain operator $C_{12}G_{23}$ (left) and the trident operator $C_{12}G_{13}G_{14}$ (right). The red path corresponds to the path generating the least intermediates (i.e. requiring the smallest number of classes in order to compute a given class). Different paths are compared in Table 3 for various two-, three- and four-electron operators, where we have reported the number of intermediates generated by each path for various integral classes.

Taking the 3-chain operator $C_{12}G_{23}$ as an example, one can see that, to compute a $[ppp]$ class, it is more advantageous to build momentum over center $\mathbf{A}_{3}$, then over centers $\mathbf{A}_{2}$, and finally over center $\mathbf{A}_{1}$ using VRRs with 4, 6 and 6 terms, respectively. The alternative path corresponding to building momentum over $\mathbf{A}_{3}$, $\mathbf{A}_{1}$, and then $\mathbf{A}_{2}$ with 4-, 5- and 7-term VRRs is slightly more expensive for a $[ppp]$ class but becomes affordable for high angular momentum classes. For both paths, using the TRR instead of the last VRR implies a large increase in the number of intermediates.

For the trident operator, we successively build angular momentum over $\mathbf{A}_{4}$, $\mathbf{A}_{3}$, $\mathbf{A}_{1}$ and $\mathbf{A}_{2}$ using VRRs with 4, 6, 8 and 7 terms. The pathway using VRRs with 4, 6, 6, and 9 terms is more expensive due to the large number of terms of the VRR building up momentum over the last center. Again, using the TRR instead of the last VRR significantly increases the number of intermediates.

The path involving the minimal number of intermediates is given in Table 3 for various two-, three- and four-electron operators. It is interesting to point out that it is never beneficial to use the TRR derived in Eq. (31) (see Sec. V.2).

One can easily show that, for operators involving the Coulomb operator, the number of intermediates required to compute a $n$-electron integral class $[a\ldots a]$ increases as $\order{a^{n+1}}$ for the VRR-only paths (see Table 3). This number is reduced to $\order{a^{n}}$ if one uses the TRR to build up angular momentum on the last center. However, the prefactor is much larger and the crossover happens for extremely high angular momentum for three- and four-electron integrals. For “pure” GG operators, such as $G_{12}$ or $G_{13}G_{23}$, the number of intermediates required to compute a class $[a\ldots a]$ increases as $\order{a^{n}}$ for any type of paths.

Finally, we note that the optimal path for the trident $C_{12}G_{13}G_{14}$ and the 4-chain $C_{12}G_{13}G_{34}$ is similar, thanks to their similar structure. Indeed, these two operators can be seen as two “linked” GGs ($G_{13}G_{14}$ or $G_{13}G_{34}$) interacting with the Coulomb operator $C_{12}$ (see Fig. 1), while the other 4-chain operator $C_{12}G_{14}G_{23}$ can be seen as two “unlinked” GGs ($G_{14}$ and $G_{23}$) interacting with the Coulomb operator.

When angular momentum has been built over all the bra centers, following the HGP algorithm, Head-Gordon and Pople (1988) we contract $[\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}}$ (step CCCC). We can perform the contraction at this point because all of the subsequent RRs are independent of the contraction coefficients and exponents. More details about this contraction step can be found in Ref. Gill, 1994.

The last step of the algorithm (step HHHH) shifts momentum from the bra center $\mathbf{A}_{1}$, $\mathbf{A}_{2}$, $\mathbf{A}_{3}$ and $\mathbf{A}_{4}$ to the ket centers $\mathbf{B}_{1}$, $\mathbf{B}_{2}$, $\mathbf{B}_{3}$ and $\mathbf{B}_{4}$ using the two-term HRRs given by Eq. (32) in Sec. V.3.

VII Concluding remarks

We have presented the three main ingredients to compute three- and four-electron integrals involving GGs. Firstly, a straightforward method to compute the FIs is given. Secondly, scaling-consistent UBs are reported, as they allow to evaluate only the $\order{N^{2}}$ significant integrals in a large system. Finally, the significant integrals are computed via a recursive scheme based on vertical and horizontal RRs, which can be viewed as an extension of the PRISM late-contraction path to three- and four-electron integrals.

We believe our approach represents a major step towards an accurate and efficient computational scheme for three- and four-electron integrals. It also paves the way to contraction-effective methods for these types of integrals. In particular, an early-contraction scheme would have significant computational benefits.

Acknowledgements.

The authors thank Peter Gill for many stimulating discussions. P.F.L. thanks the Australian Research Council for a Discovery Project grant (DP140104071) and the NCI National Facility for generous grants of supercomputer time. Financial support from the Centre national de la recherche scientifique is also acknowledged.