# GPU-Accelerated Fock Matrix Computation with Efficient Reduction
## 4 Proposed Method
### GPU Parallelization of ERIs  
We begin by describing the GPU implementation of the ERI evaluation, which constitutes the first bottleneck in the Fock matrix construction. The ERI evaluation involves calculating an integral value for every combination of four basis functions, making it intuitive to parallelize by assigning each GPU thread to compute a single value of \((\mu v|\lambda\sigma)\). However, this approach is known to be inefficient due to redundant calculations of the Boys functions \(F_{j}(X)\). Therefore, we employ an efficient mapping where each thread evaluates a task, which is a set of primitive ERIs that share the same values of the Boys function. In this mapping, instead of grouping Gauss-type orbitals by basis functions, they are managed in a data structure called a primitive shell  denoted as \(S_{1},S_{2},S_{3},\ldots\) for each shell. Figure 1 shows an example for the mathematical expression of a p-orbital with the contraction degree \(K=3\). The p-orbital consists of three components: p\({}_{x}\)-, p\({}_{y}\)-, p\({}_{z}\)-orbitals, corresponding to angular momentum \((l,m,n)=(1,0,0)\), \((0,1,0)\), and \((0,0,1)\), respectively. We organize Gauss-type orbitals with common inputs for the Boys function (total angular momentum \(L\), orbital center \(A\), and exponent \(\rho_{a}\)) into a primitive shell data structure. Assigning each GPU thread to a quartet of these primitive shells minimizes the number of Boys function evaluations. A single shell quartet contains multiple primitive ERIs, each contributing to different contracted ERIs. The set of these primitive ERIs constitute a task, defined as combinations of ERIs over Gauss-type orbitals within four shells \(S_{a}\), \(S_{\beta}\), \(S_{\gamma}\), and \(S_{\delta}\) (\(0\leq\alpha,\beta,\gamma,\delta\leq N-1\)):  
\[[S_{a}S_{\beta}|S_{\gamma}S_{\delta}]=\big{\{}[ab|cd]\mid G_{a}\in S_{a},G_{b} \in S_{\beta},G_{c}\in S_{\gamma},G_{d}\in S_{\delta}\big{\}}, \tag{25}\]  
where \(N\) denotes the number of primitive shells, and the number of primitive ERIs in a task is determined by the type of the four shells. For example, we denote a task with three s-orbitals and one p-orbital as \([\mathrm{ss}|\mathrm{sp}]\) for convenience. Since the combination of the angular momentum for s-orbitals is only \((0,0,0)\), the \([\mathrm{ss}|\mathrm{sp}]\) task contains three ERIs over Using the shell data structure, the ERI calculations can be viewed as a task of evaluating primitive ERIs for all elements of an upper triangular matrix indexed by pairs of primitive shells. Figure 2 illustrates the task matrix for ERIs composed of s- and p-orbitals. We compute all unique ERI values by parallelizing this task matrix on GPUs. Each GPU thread is assigned a unique combination of four primitive shells and sequentially evaluates primitive ERIs within that task. Equation (14) can be implemented by appropriately accumulating the obtained primitive ERIs into the corresponding contracted ERIs. The task matrix can be classified into six integral types as task sub-matrices for s- and p-orbitals: \([\text{ss}|\text{ss}]\), \([\text{ss}|\text{sp}]\), \([\text{ss}|\text{pp}]\), \([\text{sp}|\text{sp}]\), \([\text{sp}|\text{pp}]\), and \([\text{pp}|\text{pp}]\), through sorting the primitive shell pairs \(S_{\alpha}S_{\beta}\) by their total angular momenta \((L_{\alpha},L_{\beta})\). According to the angular momenta, a single task in each sub-matrix has 1, 3, 9, 9, 27, and 81 primitive ERIs, respectively. In GPU parallelization, dedicated CUDA kernels can be developed for these sub-matrices to eliminate warp divergence caused by differences in the angular momenta \(l,m,n\) and contraction degrees \(K\) in Equation (14). Our GPU implementation employs the Schwarz screening of ERIs in accordance with the task matrix based on primitive shells. To determine whether to compute or skip entire primitive ERIs within a task, we use an ERI upper bound for each task \([S_{\alpha}S_{\beta}|S_{\gamma}S_{\delta}]\).