# GPU-Accelerated Fock Matrix Computation with Efficient Reduction
## 4 Proposed Method
### GPU Parallelization of ERIs
This upper bound is defined as  
\[\forall[ab|cd]\in[S_{\alpha}S_{\beta}|S_{\gamma}S_{\delta}],\quad[ ab|cd]\leq Q_{\alpha\beta}Q_{\gamma\delta}, \tag{26}\] \[Q_{\alpha\beta}=\max_{G_{\alpha}\in S_{\alpha},G_{\beta}\in S_{ \beta}}\sqrt{[ab|ab]},\] (27) \[Q_{\gamma\delta}=\max_{G_{\alpha}\in S_{\gamma},G_{\alpha}\in S _{\delta}}\sqrt{[cd|cd]}. \tag{28}\]  
where \(Q_{\alpha\beta}\) denotes the maximum value of primitive ERIs within each task corresponding to diagonal elements of the task matrix in Figure 2. The upper bound factors \(Q_{\alpha\beta}\) are calculated in advance with a computational complexity of \(O(N^{2})\). For each shell quartet \([S_{\alpha}S_{\beta}|S_{\gamma}S_{\delta}]\), if its upper bound in Equation (26) is less than the predefined cutoff threshold \(\theta\), we omit the computation of all primitive ERIs within a task. However, divergent screening results among threads within a warp lead to performance degradation in parallelization, as threads that could skip tasks become idle while waiting for the other threads in the same warp to complete their ERI computations. To avoid this situation, we sort the pairs of primitive shells \([S_{\alpha}S_{\beta}|\) within each segment, \([\text{ss}]\), \([\text{sp}]\), and \([\text{pp}]\) according to the upper bound factors \(Q_{\alpha\beta}\) ensuring the uniformed screening result within warps. Furthermore, even if a task cannot be skipped, additional screening based on Equation (20) is applied to individual primitive ERIs to omit as many ERI computations as possible. Figure 1: Data structure of primitive shells. Efficient construction of the Fock matrix in Equation (24) based on the primitive shell data structure presents a nontrivial challenge for GPU implementation. The Fock matrix is allocated in global memory, which involves high latency to access, necessitating minimization of the access frequency. GPU computing typically employs hierarchical reduction strategies to reduce the global memory accesses, utilizing registers within each thread and shared memory within each CUDA block. Based on this concept, it is ideal to implement aggregation of integral values using shared memory in addition to thread-local reduction in registers. However, in the primitive shell-based GPU implementation, utilizing shared memory as an intermediate buffer for the Fock matrix is challenging. The contributions of multiple primitive integrals computed by each CUDA block to contracted ERIs and ultimately to Fock matrix elements are input-dependent and cannot be statically determined. Consequently, an appropriate mapping mechanism of the write destination for each ERIs is required to manage the limited shared memory space, which cannot store the entire Fock matrix. Additionally, sorting the pairs of primitive shells to avoid warp divergence in Schwarz screening further exacerbates the irregular memory access pattern. Figure 2: ERI task matrix based on the primitive shell data structure for GPU parallelization. The following proposed method using replicated Fock matrices offers a simple yet effective approach that is readily applicable to such complex shell-structured GPU implementations.