# Efficient GPU Implementation of the McMurchie-Davidson Method for Shell-Based ERI Computations
## 3 McMurchie-Davidson Method  
Our approach is centered on the efficiency frontiers and the acceleration prospects of the McMurchie-Davidson (MD) method [10]. Thus, this section describes the preliminaries and key ideas behind the approach to efficiently compute electron repulsion integrals over Gaussian functions (GTO-ERIs) by the MD method. Generally speaking, the MD method uses Hermite Gaussian functions and auxiliary (tabulated) functions through an incremental approach, whose elegant formalism offers a high degree of versatility to parallelization schemes in GPUs. For simplicity and without loss of generality of our exposition, let us consider the GTO-ERI \([ab|cd]\) composed of the following four Gaussian-type orbitals (GTOs), each of which neglects normalization constants, as follows:  
\[G_{a}(\mathbf{r}) =(x-A_{x})^{a_{x}}(y-A_{y})^{a_{y}}(z-A_{z})^{a_{z}}\exp\Bigl{(}- \alpha|\mathbf{r}-\mathbf{A}|^{2}\Bigr{)}, \tag{10}\] \[G_{b}(\mathbf{r}) =(x-B_{x})^{b_{x}}(y-B_{y})^{b_{y}}(z-B_{z})^{b_{z}}\exp\Bigl{(}- \beta|\mathbf{r}-\mathbf{B}|^{2}\Bigr{)},\] (11) \[G_{c}(\mathbf{r}) =(x-C_{x})^{c_{x}}(y-C_{y})^{c_{y}}(z-C_{z})^{c_{z}}\exp\Bigl{(}- \gamma|\mathbf{r}-\mathbf{C}|^{2}\Bigr{)},\] (12) \[G_{d}(\mathbf{r}) =(x-D_{x})^{d_{x}}(y-D_{y})^{d_{y}}(z-D_{z})^{d_{z}}\exp\Bigl{(}- \delta|\mathbf{r}-\mathbf{D}|^{2}\Bigr{)}. \tag{13}\]  
Figure 2: Example of the relation between the Basis-ERIs and the GTO-ERIs. The row (column) directions represent the _bra_ (_ket_) Basis-ERIs. Each cell in the upper triangular matrix corresponds to a single Basis-ERI. Following the MD method [10], each of the GTO-ERIs is calculated as follows:  
\[[ab|cd]=\frac{2\pi^{5/2}}{pp^{\prime}\sqrt{p+p^{\prime}}}\sum_{t=0}^{a_{x}+b_{x}} \sum_{u=0}^{a_{y}+b_{y}}\sum_{z=0}^{a_{z}+b_{z}}\sum_{t^{\prime}=0}^{c_{z}+d_{x}} \sum_{u^{\prime}=0}^{c_{y}+d_{y}}\sum_{y^{\prime}=0}^{c_{z}+d_{z}}E_{r^{\prime},u^{\prime},r^{\prime}}^{t,u,v}R_{t+t^{\prime},u+u^{\prime},p+v^{\prime}}^{0}. \tag{14}\]  
Furthermore, in Equation (14):  
\[E_{t,u,v}^{t,u,v} =E_{t,u,v}^{a,b}E_{t^{\prime},u^{\prime},v^{\prime}}^{t,v^{\prime} }(-1)^{t^{\prime}+u^{\prime}+v^{\prime}}, \tag{15}\] \[E_{t,u,v}^{a,b} =E_{t}^{a_{y},b_{x}}E_{u}^{a_{y}b_{y}}E_{v}^{a_{z}b_{z}},\] (16) \[E_{t^{\prime},u^{\prime},v^{\prime}}^{c,d} =E_{t^{\prime}}^{c_{x},d_{x}}E_{u^{\prime}}^{c_{y}d_{x}}E_{v^{ \prime}}^{c_{z}d_{z}}, \tag{17}\]  
where \(p=\alpha+\beta,p^{\prime}=\gamma+\delta\). The recurrence formulas \(E_{t}^{a_{x},b_{x}},E_{u}^{a_{y}b_{y}}\), and \(E_{v}^{a_{z}b_{z}}\) are terms related to the GTOs \(G_{a}\) and \(G_{b}\). For example, \(E_{t}^{a_{x},b_{x}}\) is defined by the following recurrence equations:  
\[E_{t}^{i,j} =\frac{1}{2p}E_{t-1}^{i-1,j}-(A_{x}-B_{x})\frac{q}{\alpha}E_{t}^{i -1,j}+(t+1)E_{t+1}^{i-1,j}, \tag{18}\] \[G_{b}(\mathbf{r}) =(x-B_{x})^{b_{x}}(y-B_{y})^{b_{y}}(z-B_{z})^{b_{z}}\exp\Bigl{(}- \beta|\mathbf{r}-\mathbf{B}|^{2}\Bigr{)},\]0 (19) \[G_{b}(\mathbf{r}) =(x-B_{x})^{b_{x}}(y-B_{y})^{b_{y}}(z-B_{z})^{b_{z}}\exp\Bigl{(}- \beta|\mathbf{r}-\mathbf{B}|^{2}\Bigr{)},\]1 (20) \[G_{b}(\mathbf{r}) =(x-B_{x})^{b_{x}}(y-B_{y})^{b_{y}}(z-B_{z})^{b_{z}}\exp\Bigl{(}- \beta|\mathbf{r}-\mathbf{B}|^{2}\Bigr{)},\]2  
In the above relations,  
\[G_{b}(\mathbf{r}) =(x-B_{x})^{b_{x}}(y-B_{y})^{b_{y}}(z-B_{z})^{b_{z}}\exp\Bigl{(}- \beta|\mathbf{r}-\mathbf{B}|^{2}\Bigr{)},\]3  
Also, the terms \(E_{u}^{a_{y},b_{y}}\) and \(E_{v}^{a_{y}b_{z}}\) can be calculated using the coordinates \(y\) and \(z\), respectively, and the terms \(E_{t^{\prime}}^{c_{x},d_{x}},E_{u^{\prime}}^{c_{y},d_{y}}\), and \(E_{v^{\prime}}^{c_{y},d_{z}}\) can be calculated for \(G_{c}\) and \(G_{d}\) in the same fashion. Furthermore, although evaluating the recurrence formula \(E_{t}^{i,j}\) is computationally trivial and inexpensive in most cases, implementing the terms of \(E_{t}^{i,j}\) as a recursive function in a large number of basis sets is expected to contribute to the overall computational load and bottleneck due to the large number of computation paths.