A Stochastic Formulation of the Resolution of Identity: Application to Second Order M{\o}ller-Plesset Perturbation Theory

TL;DR

This paper introduces a stochastic orbital approach to the resolution of identity for electron repulsion integrals, enabling efficient MP2 calculations with improved scaling and performance on water clusters.

Contribution

It presents a novel stochastic RI method with multiple orbitals that reduces MP2 computational scaling and outperforms traditional methods on water clusters.

Findings

Achieves $N^{2.4}$ scaling for water clusters

Outperforms MP2 for clusters with 21 water molecules

Demonstrates efficiency of stochastic RI-MP2 approach

Abstract

A stochastic orbital approach to the resolution of identity (RI) approximation for 4-index 2-electron electron repulsion integrals (ERIs) is presented. The stochastic RI-ERIs are then applied to M\o ller-Plesset perturbation theory (MP2) utilizing a \textit{multiple stochastic orbital approach}. The introduction of multiple stochastic orbitals results in an $N^3$ scaling for both the stochastic RI-ERIs and stochastic RI-MP2. We demonstrate that this method exhibits a small prefactor and an observed scaling of $N^{2.4}$ for a range of water clusters, already outperforming MP2 for clusters with as few as 21 water molecules.

Tables

Table 1. Table 1: MP2 and sRI-MP2 parameters and results for the water cluster test set. N e subscript 𝑁 𝑒 N_{e} = number of correlated electrons. MP2 and sRI-MP2 correlation energies per electron in Hartree. Error and standard error per electron in kcal/mol. Basis set: cc-pVDZ. Auxiliary basis set: cc-pVDZ-RI.

$N_e$	$N_{AO}$	$N_{aux}$	MP2	sRI-MP2	Error/$N_e$	Std Error/$N_e$	$N_{pairs}$
64	200	768	-0.0270	-0.0281	0.6750	0.8440	200
168	500	2016	-0.0268	-0.0261	0.3947	0.8422	200
256	800	3072	-0.0268	-0.0269	0.0577	0.6579	200
416	1300	4992	-0.0269	-0.0268	0.0426	1.0825	200
624	1950	7488	-0.0270	-0.0283	0.8304	1.1841	200
888	2775	10656		-0.0281		1.0755	200

Equations

$$\begin{array}{rl}(\alpha \beta |\gamma \delta )& =\iint d{r}_{1}d{r}_2\frac{{\chi }_{\alpha }(r_1){\chi }_{\beta }(r_1){\chi }_{\gamma }(r_2){\chi }_{\delta }(r_2)}{r_{12}}\end{array}$$

$$(\alpha \beta |A)=\iint d{r}_{1}d{r}_2\frac{{\chi }_{\alpha }(r_1){\chi }_{\beta }(r_1){\chi }_A(r_2)}{r_{12}}$$

$$V_{AB}=\iint d{r}_{1}d{r}_2\frac{{\chi }_A(r_1){\chi }_B(r_2)}{r_{12}}.$$

$$\begin{split}&(\alpha\beta|\gamma\delta)\approx\sum_{AB}^{N_{aux}}(\alpha\beta|A)[V^{-1}]_{AB}(B|\gamma\delta)\\ &=\sum_{Q}^{N_{aux}}\Big{[}\sum_{A}^{N_{aux}}(\alpha\beta|A)[V^{-\frac{1}{2}}]_{AQ}\Big{]}\Big{[}\sum_{B}^{N_{aux}}[V^{-\frac{1}{2}}]_{QB}(B|\gamma\delta)\Big{]}.\end{split}$$

$$K_{\alpha \beta }^Q\equiv \sum _A^{N_{aux}}(\alpha \beta |A)V_{AQ}^{-\frac{1}{2}},$$

$$(\alpha \beta |\gamma \delta )\approx \sum _Q^{N_{aux}}{K}_{\alpha \beta }^Q{K}_{\gamma \delta }^Q.$$

$$\begin{array}{rl}{K}_{p\gamma }^Q& =\sum _{\alpha }^{N_{AO}}{C}_{\alpha }^p{K}_{\alpha \gamma }^Q\\ K_{pq}^Q& =\sum _{\gamma }^{N_{AO}}{C}_{\gamma }^q{K}_{p\gamma }^Q.\end{array}$$

$$\begin{split}\Big{<}\theta\otimes\theta^{T}\Big{>}_{\xi}=I,\end{split}$$

$$\begin{split}\Big{<}\theta\otimes\theta^{T}\Big{>}_{\xi}&=\frac{1}{N_{s}}\sum_{\xi=1}^{N_{s}}\theta^{\xi}\otimes(\theta^{\xi})^{T}\equiv\begin{pmatrix}\left<\theta_{1}\theta_{1}\right>_{\xi}&\left<\theta_{1}\theta_{2}\right>_{\xi}\\ \left<\theta_{2}\theta_{1}\right>_{\xi}&\left<\theta_{2}\theta_{2}\right>_{\xi}\end{pmatrix}.\end{split}$$

$$\begin{split}&(\alpha\beta|\gamma\delta)\approx\sum_{PQ}^{N_{aux}}\sum_{AB}^{N_{aux}}(\alpha\beta|A)V^{-\frac{1}{2}}_{AP}I_{PQ}V^{-\frac{1}{2}}_{QB}(B|\gamma\delta)\\ &=\sum_{PQ}^{N_{aux}}\sum_{AB}^{N_{aux}}(\alpha\beta|A)V^{-\frac{1}{2}}_{AP}\left(\left<\theta\otimes\theta^{T}\right>_{\xi}\right)_{PQ}V^{-\frac{1}{2}}_{QB}(B|\gamma\delta)\\ &=\Big{<}\left[\sum_{A}^{N_{aux}}(\alpha\beta|A)\sum_{P}^{N_{aux}}V^{-\frac{1}{2}}_{AP}\theta_{P}\right]\\ &\hskip 60.00009pt\times\left[\sum_{B}^{N_{aux}}(B|\gamma\delta)\sum_{Q}^{N_{aux}}\theta_{Q}^{T}V^{-\frac{1}{2}}_{QB}\right]\Big{>}_{\xi},\\ \end{split}$$

$$R_{\alpha \beta }^{\xi }=\sum _A^{N_{aux}}(\alpha \beta |A)\left[\sum _P^{N_{aux}}[V_{AP}^{-\frac{1}{2}}{\theta }_P^{\xi }]\right]\equiv \sum _A^{N_{aux}}(\alpha \beta |A)L_A^{\xi }.$$

$$(\alpha \beta |\gamma \delta )\approx \frac{1}{N_s}\sum _{\xi }{R}_{\alpha \beta }^{\xi }{R}_{\gamma \delta }^{\xi }\equiv {⟨R_{\alpha \beta }{R}_{\gamma \delta }⟩}_{\xi }.$$

$$\begin{array}{rl}{R}_{p\beta }^{\xi }& =\sum _{\alpha }^{N_{AO}}{C}_{\alpha }^p{R}_{\alpha \beta }^{\xi }\\ R_{pq}^{\xi }& =\sum _{\beta }^{N_{AO}}{C}_{\beta }^q{R}_{p\beta }^{\xi },\end{array}$$

$$E_{MP2}=\sum _{abij}\frac{(ai|bj)[2(ai|bj)-(bi|aj)]}{{\epsilon }_i+{\epsilon }_j-{\epsilon }_a-{\epsilon }_b},$$

$$\begin{array}{rl}{E}_{sRI-MP2}& =\sum _{abij}\frac{{⟨R_{ai}^{\xi }{R}_{bj}^{\xi }⟩}_{\xi }\left[2{⟨R_{ai}^{\xi }{R}_{bj}^{\xi }⟩}_{\xi }-{⟨R_{aj}^{\xi }{R}_{bi}^{\xi }⟩}_{\xi }\right]}{{\epsilon }_i+{\epsilon }_j-{\epsilon }_a-{\epsilon }_b}.\end{array}$$

$$E_{sRI-MP2}={⟨\sum _{abij}\frac{R_{ai}^{\xi }{R}_{bj}^{\xi }[2R_{ai}^{{\xi }^{\prime }}{R}_{bj}^{{\xi }^{\prime }}-R_{aj}^{{\xi }^{\prime }}{R}_{bi}^{{\xi }^{\prime }}]}{{\epsilon }_i+{\epsilon }_j-{\epsilon }_a-{\epsilon }_b}⟩}_{\xi {\xi }^{\prime }}$$

$$\begin{split}E_{sRIMP2}&=\int_{0}^{\infty}\sum_{abij}\Big{<}\Big{[}2(R_{ai}^{\xi}R_{ai}^{\xi^{\prime}})(R_{bj}^{\xi}R_{bj}^{\xi^{\prime}})\\ &\hskip 10.00002pt-(R_{ai}^{\xi}R_{aj}^{\xi^{\prime}})(R_{bj}^{\xi}R_{bi}^{\xi^{\prime}})\Big{]}e^{-(\varepsilon_{i}+\varepsilon_{j}-\varepsilon_{a}-\varepsilon_{b})t}\Big{>}_{\xi\xi^{\prime}}dt\\ &=\int_{0}^{\infty}\left<2A(t)^{2}-Tr[E(t)^{2}]\right>_{\xi\xi^{\prime}}dt,\end{split}$$

$$\begin{array}{rl}A(t)& =\sum _i^{N_{occ}}\sum _a^{N_{virt}}{e}^{-({\epsilon }_i-{\epsilon }_a)t}{R}_{ai}^{\xi }{R}_{ai}^{{\xi }^{\prime }}\\ E(t{)}_{ij}& =\sum _a^{N_{virt}}{e}^{-({\epsilon }_i-{\epsilon }_a)t}{R}_{ai}^{\xi }{R}_{aj}^{{\xi }^{\prime }}.\end{array}$$

Full text

I Introduction

The vast majority of ab initio electronic structure methods require the calculation of 4-index electron repulsion integrals (ERIs). In fact, in some instances, when atom-centered gaussian basis sets are used the calculation of these integrals and their transformation from the atomic orbital (AO) to the molecular orbital (MO) basis is the computational bottleneck, e.g. Møller-Plesset perturbation theory (MP2). An appreciable reduction in the computational prefactor may be obtained through the resolution of identity (RI) approximation, also known as the density fitting approximation.Whitten (1973); Dunlap (1983); Dunlap et al. (1979); Vahtras et al. (1993); Feyereisen et al. (1993) The RI approximation expresses the 4-index ERIs in terms of 2-index and 3-index ERIs, the former being evaluated in an auxiliary basis and the latter as a combination of the AO and auxiliary basis sets. As only 2- and 3-index ERIs are needed, the RI approximation reduces the total number of integrals to be calculated and transformed. Today it has become common practice to apply the RI approximation to 4-index ERIs in order to lower the computational prefactor. However, in spite of these benefits, the assembly of the approximate ERIs scales as $O(N^{5})$ and therefore the scaling remains unaltered. Recent work focused on mitigating the high computational cost associated with the 4-index ERIs through the application of the tensor decomposition technique known as tensor hypercontractionHohenstein et al. (2012a); Parrish et al. (2012); Hohenstein et al. (2012b) has resulted in flexible factorization of the ERIs and reduced scaling.

As an alternative to reduced scaling techniques focused on the ERIs, stochastic approaches to performing traditional electronic structure calculations have proven effective in reducing the high computational cost.Thom and Alavi (2007); Ohtsuka and Nagase (2008); Booth et al. (2009); Booth and Alavi (2010); Manni et al. (2016); Thom (2010); Spencer and Thom (2016); Willow et al. (2012, 2013); Willow and Hirata (2014); Neuhauser et al. (2013a, b); Baer et al. (2013); Neuhauser et al. (2014a, b); Ge et al. (2014); Rabani et al. (2015); Gao et al. (2015); Neuhauser et al. (2016) There are many successful stochastic techniques that can handle increasingly larger systems. We note, for example, that in certain situations the Full Configuration-Interaction Quantum Monte Carlo approach can handle systems with tens of electrons Thom and Alavi (2007); Ohtsuka and Nagase (2008); Booth et al. (2009); Booth and Alavi (2010) Likewise, Auxiliary-Field Monte which replaces the two-body interaction by an interaction with fluctuating densities and the fixed-node approximationZhang et al. (1995) when combined with the Shifted-Contour approachRom et al. (1997) give excellent results for systems with tens of electrons.Shee et al. For large systems containing hundreds or thousands of electrons several of the authors have developed stochastic methods for DFT and TDDFT Baer et al. (2013); Neuhauser et al. (2014c); Gao et al. (2015); Neuhauser et al. (2016), MP2 Neuhauser et al. (2013a); Ge et al. (2014), GF2Neuhauser et al. , GW Rabani et al. (2015); Vlček et al. (2016a, tted, b) and the Bethe-Salpeter equationRabani et al. (2015).

Given the success of the RI approximation and stochastic electronic structure methods it is therefore conceivable that methods that bring together the strengths of both approaches could prove extremely beneficial. In this letter, we present a hybrid approach, stochastic resolution of identity (sRI), that (i) lowers the computational scaling of the RI approximation to the 4-index ERIs and (ii) decouples pairs of indices within the 4-index ERI expression, a general feature capable of bringing about additional method-specific reductions in scaling. We apply the sRI approximation to the time-integrated MP2 expression obtaining an observed scaling of $O(N^{2.4})$.

II Theory

We use the usual notation, where the occupied, virtual and general set of MOs are represented by the indices $i,j,k,\dots$; $a,b,c,\dots$ and $p,q,r,\dots$ respectively. The AO Gaussian basis functions are represented by $\chi_{\alpha}(r)$ and greek indices $\alpha,\beta,\gamma,\delta,\dots$ while the auxiliary basis functions are represented by the indices $A,B,\dots$. Finally, the total number of AO basis functions, auxiliary basis functions, occupied MOs and virtual MOs are $N_{AO}$, $N_{aux}$, $N_{occ}$ and $N_{virt}$ respectively. Further, both $N_{aux}$ and $N_{AO}$ are proportional to the system size with $N_{aux}$ typically 3 to 6 times $N_{AO}$.

II.1 Deterministic Resolution of Identity

The 4-, 3- and 2-index ERIs are defined as:

[TABLE]

[TABLE]

[TABLE]

The approximate 4-index RI-ERIs are then expressed symmetrically in terms of the lower-rank integrals according to:

[TABLE]

Defining

[TABLE]

yields

[TABLE]

Summations over $A$ and $B$ (Eq. (2) and (3)) are usually performed beforehand and their contractions, $K_{\alpha\beta}^{Q}$ and $K_{\gamma\delta}^{Q}$, scale as $O(N_{AO}^{2}N_{aux})$ while the construction of $V^{-\frac{1}{2}}$ scales as $O(N_{aux}^{3})$. By expressing Eq. (2) in terms of $K_{\alpha\beta}^{Q}$ and $K_{\gamma\delta}^{Q}$ (Eq. (5)) the approximate ERIs now scale as $O(N_{AO}^{4}N_{aux})$

[TABLE]

ERIs are most often used in the MO basis and their transformation to the AO is done in a two step process with both the first and the second transformations (Eq. (6)) costing $O(N_{AO}^{3}N_{aux})$.

[TABLE]

According to Eq. (5) the cost of computing the RI-ERIs scale as $O(N_{AO}^{4}N_{aux})$; however, the total number of integrals that must be calculated grows only as $O(N_{AO}^{2}N_{aux})$. Since both $N_{AO}$ and $N_{aux}$ are dependent on the system size, the principle advantage of the RI approximation is therefore its ability to reduce the total number of integrals that must be calculated and stored while maintaining the same overall scaling.

II.2 Stochastic Resolution of Identity

The stochastic RI approximation we develop here utilizes the same set of 2- and 3-index ERIs while introducing an additional set of $N_{s}$ stochastic orbitals, $\{\theta^{\xi}\}$, $\xi=1,2,\cdots,N_{s}$. The stochastic orbitals are defined as arrays of length $N_{aux}$ with randomly selected elements $\theta_{A}^{\xi}=\pm 1$. The stochastic orbitals have the following property:

[TABLE]

where we have denoted the stochastic average over $N_{s}$ stochastic orbitals by \big{<}\big{>}_{\xi}. To better illustrate this, consider the case where the set $\{\theta^{\xi}\}$ contains $N_{s}$ elements, where each array $\theta^{\xi}$ is of length $N_{aux}=2$. The resulting stochastic average is then

[TABLE]

The individual matrix elements may be grouped as diagonal and off-diagonal elements. The stochastic element-by-element average of the diagonal elements, $\left<\theta_{A}\theta_{A}\right>_{\xi}$, is 1 and the stochastic average of the off-diagonal elements, $\left<\theta_{A}\theta_{B}\right>_{\xi}$, converges to 0 as $N_{s}\to\infty$, due to the random oscillations of $\theta_{A}^{\xi}\theta_{B}^{\xi}$ between $\pm 1$. The above example shows that the introduction of an identity matrix can be recast as the stochastic average over outer products of stochastic orbitals and is the underlying principle of the stochastic resolution of identity method.

The deterministic RI-ERIs in Eq. (2) are expressed symmetrically in terms of the 2-index and 3-index ERI matrix elements with the symmetric parts being coupled through a summation over the index $Q$. Inserting the stochastic identity matrix we obtain the expression for the sRI-ERIs:

[TABLE]

where $\left(\left<\theta\otimes\theta^{T}\right>_{\xi}\right)_{PQ}$ is the $PQ^{th}$ element of the stochastic identity matrix. We now define the $\xi^{th}$ elements of the stochastic average as

[TABLE]

With this definition, the ERI in the AO basis (Eq. (9)) is now given by a stochastic average, an $O(N_{s}N_{AO}^{4})$ step:

[TABLE]

Calculation of the $L^{\xi}_{A}$ terms in Eq. (10) scales as $O(N_{aux}^{2}N_{s})$ while the overall computational scaling of the $R^{\xi}$ matrices is $O(N_{s}N_{AO}^{2}N_{aux})$. This is similar to the deterministic RI components $K_{\alpha\beta}^{Q}$ and $K_{\gamma\delta}^{Q}$ but with an additional prefactor of $N_{s}$.

The transformation to the MO basis is given by

[TABLE]

and is a two step process with both transformation steps scaling as $O(N_{s}N_{AO}^{3})$ compared to the deterministic transformation that costs $O(N_{aux}N_{AO}^{3})$.

The stochastic error of the elements of the identity matrix and therefore the error of the ERIs is governed by the number of stochastic orbitals, $N_{s}$ as can be seen from Eq. (8). Since it is the length of stochastic arrays, $N_{aux}$, that increases with the system size rather than the number of stochastic orbitals, $N_{s}$ is expected to have little size dependence. We will show for a set of water clusters that $N_{s}$ remains approximately constant as a function of systems size for a fixed statistical error. Thus, the transformation from the AO to MO basis scales as $O(N_{AO}^{3})$, and the 4-index ERI assembly as $O(N_{AO}^{4})$ a factor of $N_{aux}/N_{s}$ less than deterministic RI.

II.3 Stochastic Resolution of Identity MP2

As we have stated above in some instances the sRI approximation may lead to an additional decrease in scaling due to the decoupling of indices and we now demonstrate this for MP2. The MP2 energy expression for a closed shell system may be written as

[TABLE]

and implementing the sRI approximation we obtain a similar expression for sRI-MP2

[TABLE]

Although Eq. (13) is an $O(N_{occ}^{2}N_{virt}^{2})$ step, MP2 scales as $O(N_{occ}N_{AO}^{4})$ because of the 4-index ERI transformation, while RI-MP2 scales as $O(N_{occ}^{2}N_{virt}^{2}N_{aux})$ due to the reconstruction step in Eq. (5). Similarly, with the naive application of the sRI approximation in Eq. (14) one sees that sRI-MP2 is expected to scale as $O(N_{s}N_{occ}^{2}N_{virt}^{2})$. However, with the introduction of a second stochastic orbital in conjunction with Almlöf’sHäser and Almlöf (1992) time-integrated decomposition of the energy denominator, it is possible reduce the overall cost to that of the $R$ matrices (Eq. (10)). First the sRI-MP2 energy expression is written in terms of two rather than one stochastic orbital denoted by $\xi$ and $\xi^{\prime}$ in Eq. (15).

[TABLE]

The introduction of the second stochastic orbital doubles the number of $R^{\xi}$ matrices while leaving the number of elements in the stochastic average unchanged. The use of two stochastic orbitals is denoted by \big{<}\big{>}_{\xi\xi^{\prime}}. The modest increase in the computational prefactor and memory requirements of sRI-MP2 is extremely advantageous as it allows the stochastic average to be taken over the entire sRI-MP2 energy expression rather than individual integral pairs decoupling indices in the numerator. The numerator may now be rearranged in terms of products of the form $R^{\xi}_{ai}R^{\xi^{\prime}}_{ai}$ and $R^{\xi}_{ai}R^{\xi^{\prime}}_{aj}$ and the denominator rewritten as a time integral resulting in the time-integrated sRI-MP2 expression of Eq. (16).

[TABLE]

where

[TABLE]

The quantity $A(t)$ scales as $O(N_{occ}N_{virt})$ and the matrix $E(t)$ as $O(N_{occ}^{2}N_{virt})$. The overall scaling for the energy expression is $O(N_{s}N_{t}N_{occ}^{2}N_{virt})$, and in the case of small prefactors, $N_{s}$ and $N_{t}$, becomes $O(N_{occ}^{2}N_{virt})$.

III Results and Discussion

To study the observed scaling, stochastic errors and the impact of the prefactors, $N_{s}$ and $N_{t}$, on the sRI-MP2 method, we selected a test set of water clusters consisting of 8, 21, 32, 52, 78 and 111 water molecules. The sRI-ERI and time-integrated sRI-MP2 routines are implemented in a development version of the NWChem 6.6 package of computational chemistry tools.Valiev et al. Deterministic MP2 calculations were performed with the NWChem semi-direct MP2 module. Dunning’s correlation consistent basis sets of double zeta quality, cc-pVDZ,Dunning (1989) were used for all calculations and the corresponding, cc-pVDZ-RI, auxiliary basisWeigend et al. (2002); Hattig (2005) used in sRI-MP2 calculations. Schwarz integral screening was applied to all 4-, 3- and 2-index ERIs. All benchmark calculations were performed with the National Energy Research Scientific Computing Center resource Cori using a single Haswell compute node and 30 computational cores.

The results are listed in Table 1 where deterministic MP2 and sRI-MP2 correlation energies per electron are given in Hartree and the error in the correlation energy per electron and standard error of correlation energy per electron given in units of kcal/mol. As mentioned previously the computationally demanding step of the sRI approximation is the construction of the $R^{\xi}$ matrices which scales as $O(N_{s}N_{AO}^{2}N_{aux})$ while the sRI-MP2 energy expression is an $O(N_{s}N_{t}N_{occ}^{2}N_{virt})$ step. For the given test set ten quadrature points were found to be sufficient for the energy denominator decomposition. Therefore, the observed scaling of the method is dependent on $N_{s}$ remaining small with respect to the system size. The results listed in Table 1 show that using $N_{s}=200$ is sufficient to produce errors below 1 kcal/mol per electron for all systems within the test set.

The observed MP2 and sRI-MP2 timings per core are plotted in Figure 1. For a system of eight water molecules the sRI-MP2 method is 3.5 times more expensive than the deterministic MP2. However, for systems above 161 correlated electrons (approximately 21 water molecules with $N_{e}$ = 168) the computational cost of sRI-MP2 drops below that of MP2 with an observed scaling of $O(N^{2.4})$.

If the extent of the sRI-MP2 capabilities were limited to converging the correlation energy per electron to within a given threshold of the deterministic results, sRI-MP2 would be of limited utility as in most practical applications to large systems, e.g materials, it is necessary to accurately calculate relative energies and forces. Specifically we must verified that a constant per-electron error leads to constant (small) error in the forces and relative energies. As an initial investigation we calculated the potential energy curve and numerical gradients of a system of two water molecules in a hydrogen-bonded configuration as a function of the internuclear distance along the hydrogen-bond coordinate. function of the intermolecular coordinate. The potential energy curve was generated on an equally spaced grid with $\Delta R$ = 0.2Å  and then fitted with a cubic spline to calculate the forces. We found that the most efficient sampling method to generate reasonably accurate potential energy curves was to average over $n$ sRI-MP2 calculations each performed with $N_{s}^{\prime}$ stochastic samples such that $N_{s}=nN_{s}^{\prime}$. The potential energy curves, MP2 and sRI-MP2 correlation energy and forces are plotted in Figure 2 for $N_{s}=$ 400 and 600 with $N_{s}^{\prime}=100$. This averaging approach resulted in faster convergence to the deterministic result with errors in the total relative energies of less than 1 kcal/mol for $N_{s}$ set at 400 and 600. From the correlation energies plotted in Figure 2 it is clear that the MP2 and sRI-MP2 correlation energies are significant, accounting for nearly half the total relative energy at the equilibrium distance. sRI-MP2 was able to reproduce the equilibrium geometry within 0.01 Å, while the forces were found to have errors of less than 1 (kcal/mol)/ Å  in the range -0.1 Å  to 2.0 Å  with respect to the equilibrium hydrogen bonded geometry. Errors in the potential energy curves and stochastic forces increased to a maximum of 3.3 (kcal/mol)/ Å  and 8.1 (kcal/mol)/ Å  respectively when the hydrogen bond distance was shorted by 0.4 Å  with respect to the equilibrium bond distance.

To conclude, we introduced a stochastic implementation of the resolution of identity approximation that reduced the scaling of the deterministic AO to MO transformation form $O(N^{5})$ (or $O(N^{4})$ for the deterministic RI approximation) to $O(N^{3})$ and overall memory requirements to $O(N^{2})$. It was then demonstrated that with the introduction of an additional stochastic orbital the stochastic averaging may take place over more complex expressions rather than individual 4-index ERIs leading to a decoupling of indices. This led to the time-integrated sRI-MP2 with a formal scaling of $O(N^{3})$ and an observed scaling of $O(N^{2.4})$ when applied to a set of 3 dimensional systems. Given that 4-index 2-electron ERI are ubiquitous in ab initio electronic structure methods we expect the sRI approximation to be widely applicable and readily interfaced with other reduced scaling techniques.

Acknowledgements.

This work was supported by the Laboratory Directed Research and Development Program of Lawrence Berkeley National Laboratory under U.S. Department of Energy Contract No. DE-AC02-05CH11231. D. Neuhauser and R. Baer are grateful for support by the National Science Foundation Division of Materials Research and Binational Science Foundation, grant numbers 1611382 and 2015687.