# Efficient GPU Implementation of the McMurchie-Davidson Method for Shell-Based ERI Computations
## 4 Proposed Method
### Parallel Batch Computation of the Recurrence Relation \(R\)  
This section explores the efficiency frontiers of the MD method through GPU acceleration schemes. Although the MD method offers an elegant formalism with a high degree of versatility to compute electron repulsion integrals over Gaussian functions (GTO-ERIs), several challenges prevent it from attaining practical performance for molecules of interest. One of the main drawbacks is the high computational cost associated with repetitive evaluations of the Boys function, and the redundancy of the identical Boys function values in iterative calls to the recursive function \(R\). To tackle this problem, our results presented in [23] address the issues mentioned earlier by expanding the set of recursive functions \(R\) and by dividing into independent computational units termed _batches_, thereby avoiding redundant calculations. Additionally, it is possible to precompute the required values of \(R_{t,u,v}^{0}\) by evaluating each element within a _batch_.  
In what follows, we provide an overview of the _batch_-based algorithm. First, let us define \(s_{x}=a_{x}+b_{x}+c_{x}+d_{x}\), \(s_{y}=a_{y}+b_{y}+c_{y}+d_{y}\), \(s_{z}=a_{z}+b_{z}+c_{z}+d_{z}\), and \(K=s_{x}+s_{y}+s_{z}\). To compute the GTO-ERI, all \(R_{t,u,v}^{0}\) values satisfying \(0\leq t\leq s_{x}\), \(0\leq u\leq s_{y}\), \(0\leq v\leq s_{z}\) are to be calculated. Among the expanded recurrence relations \((R_{t,u,v}^{n})\), we assign the values meeting the following condition to _batch_\(k\) (\(0\leq k\leq K\)) such that \(k=t+u+v\) and \(0\leq n\leq K-k\). The number of combinations \((t,u,v)\) that satisfy \(t+u+v=k\) is  
\[\frac{(k+2)!}{k!2!}=\frac{(k+1)(k+2)}{2}.\]  
Also, _batch_\(k\) includes the recurrence \(R_{t,u,v}^{n}\) for \(0\leq n\leq K-k\). Therefore, the number of values through recurrence \(R_{t,u,v}^{n}\) in _batch_\(k\) is given by  
\[\frac{(k+1)(k+2)(K-k+1)}{2}.\]  
Similarly, the total number of values computed through the recurrence \(R_{t,u,v}^{n}\) across all _batches_ is  
Figure 3: Basic idea behind the definition of shell-based ERIs. The term \((ss|\) implies that a _bra_ consists of two s-shells, and the term \(|sp\rangle\) implies that a _ket_ consists of one s-shell and one p-shell. The integral \((ss|sp\rangle\) consists of three GTO-ERIs: \([ss|sp_{x}]\), \([ss|sp_{y}]\), and \([ss|sp_{z}]\).  
\[\sum_{k=0}^{K}\frac{(k+1)(k+2)(K-k+1)}{2}=\frac{(K+1)(K+2)(K+3)(K+4)}{24}.\]  
Furthermore, by sequentially calculating from _batch_\(0\) to _batch_\(K\), all required values of the recurrence \(R_{i,u,v}^{0}\) are obtained. The above-mentioned observations become evident from the definitions of Equations (8) and (22). Moreover, as the values rendered from recurrence \(R_{i,u,v}^{n}\) within a _batch_ are independent according to Equation (22), all values within a _batch_ can be computed in parallel. For reader's reference, in order to portray a glimpse of the kind of dependencies in the computation paths involved, Figure 4 shows the computed values through the corresponding recurrences \(R\) through the use of _batches_ for \(K=4\).