# Efficient GPU Implementation of the McMurchie-Davidson Method for Shell-Based ERI Computations
## 4 Proposed Method
### Parallel Algorithms  
The proposed parallel algorithm employs a two-stage parallelization strategy: the first stage computes the ERIs corresponding to the basis functions \((\mu v|\lambda\sigma)\) in parallel, and the second stage parallelizes the GTO-ERI algorithm \([ab|cd]\) by the MD method, which is conveniently computed in the first stage. Algorithm 1 parallelizes the computation of \(M\) basis functions by leveraging the independence of all combinations of four basis functions. Also, Algorithm 1 avoids redundant computations of the corresponding recurrences \(R\), and the results are preserved for further reuse. The inner loop of the algorithm is executed sequentially, while the GTO-ERI computation \([\mu_{e}v_{f}|\lambda_{g}\sigma_{h}]\) in line 4 utilizes the parallel GTO-ERI algorithm presented in Algorithm 2. Figure 4: Basic idea of the dependencies, denoted by arrows, behind computing the values of the corresponding recurrences \(R\) by using _batch_ concepts when \(K=4\). The values required by the MD method are highlighted in red. ```
0:\(N\) basis functions
0: Basic-ERIs for all combinations of four basis functions
1:for all \(\mu,\
u,\lambda,\sigma\) such that \(\mu\leq\
u\), \(\lambda\leq\sigma\), and \((\mu,\
u)\leq(\lambda,\sigma)\)do in parallel
2:\(t\gets 0\) // Equation (9)
3:for all \(e,f,g,h\)do
4:\(t\gets t+d_{\mu}d_{\
u\prime}d_{\lambda_{g}}d_{\alpha_{h}}\left[\mu_{e} \
u_{f}|\lambda_{g}\sigma_{h}\right]\) // Algorithm 2
5:endfor
6:\((\mu\
u|\lambda\sigma)\gets t\)
7:endfor
```  
**Algorithm 1** Parallel Basis-ERI algorithm. Algorithm 2 presents a parallel MD algorithm for computing the GTO-ERI. The algorithm first precomputes all necessary values of the recurrence \(R\). As explained earlier, a _batch_\(k\) can be computed from _batches_\(k-1\) and \(k-2\), and the _batches_ are computed sequentially from _batch_\(0\) to _batch_\(K\). Additionally, there are no dependencies among the values of \(R\) within each _batch_, enabling seamless parallel computation. In lines 4 through 8, the obtained values from the corresponding recurrences \(R\) are used to compute Equation (2). Since the overall procedure implies a six-fold loop computation, line 6 is computed in parallel, and the resulting sum is obtained using parallel addition techniques for shared memory. Since evaluating the recurrence formula \(E_{t}^{i,j}\) is computationally less intensive than evaluating the recurrence formula \(R_{t,\mu,\
u}^{n}\), the values derived from the recurrence \(E\) are obtained by computing the recursive expressions straightforwardly. Algorithm 3 shows a parallel algorithm for Basis-ERI computation based on shells. Whereas Algorithm 1 computes the Basis-ERIs, Algorithm 3 computes the shell-based ERIs. For a given combination of four shells, the algorithm computes the shell-based ERIs in parallel, as shown in Equation (28). As explained above, it is possible to reuse the Boys functions computed in the shells. Therefore, first, in line 2, all necessary Boys function values for the GTO-ERIs in the shells are computed independently and in parallel. Afterward, in lines 3-5, the GTO-ERIs in the shell-based ERI of Equation (8) are calculated sequentially. The result of each corresponding calculation is accumulated in the value of the basis function \((\mu\
u|\lambda\sigma)\). The reader may note that the key difference between Algorithms 1 and 3 lies in the order of calculations, yet the final computed results are equivalent.