# GPU-Accelerated Fock Matrix Computation with Efficient Reduction
## 4 Proposed Method
### Distributed Atomic Reduction Through Replicated Fock Matrix Update  
Based on the GPU parallelization that computes only unique ERIs in Section 4.1, we efficiently update the Fock matrix while considering eight-fold symmetry. In our GPU implementation, each thread computes multiple primitive ERIs that contribute to different contracted ERIs. Consequently, we implement Equation (24) by directly accumulating the obtained primitive ERIs into the Fock matrix. Even with partial aggregation of primitive ERIs in thread-local buffers, \(O(M^{4})\) atomic additions to the \(O(M^{2})\) memory space of the Fock matrix are required. These exclusive additions are implemented using the atomicAdd operation in CUDA. The atomicAdd operation is a function that increments a variable stored at a specified address, executing the read and write as an indivisible process. The atomicAdd operations competing at the same address by multiple threads are processed sequentially, resulting in performance degradation in parallel computation. To avoid overhead caused by the memory contention, we distribute the atomicAdd operations into multiple replicated Fock matrices. Figure 3 illustrates the distributed atomic reduction method based on the replicated Fock matrices for an example of the \([\mathrm{ss}|\mathrm{sp}]\) sub-matrix. We decompose a single Fock matrix \(F\) into multiple replicas \(F_{k}(0\leq k\leq N_{\mathrm{Fock}}-1)\) and significantly reduce the memory contention by having each thread add values to different replicas. Each replica serves as an intermediate buffer of the same size as the original, and the desired Fock matrix is constructed by summing all the replicas. The index of replicas \(k\) contributed by each thread with index \(t\) is determined by  
\[k=t\bmod N_{\mathrm{Fock}}, \tag{29}\]  
ensuring that adjacent threads, which are prone to address conflicts, perform the atomicAdd operations on different replicas. Finally, the desired Fock matrix \(F\) is obtained by summing the \(N_{\mathrm{Fock}}\) replicas, each of which accumulates the product of the primitive ERIs and density matrix elements. We consider not only the eight-fold symmetry of ERIs but also the symmetry of the Fock matrix, ensuring that the updates according to Equation (24) are applied only to the upper triangular elements of \(F\). To detail our GPU implementation, Algorithm 1 presents pseudocode for the Fock matrix construction using the distributed atomic reduction. The distributed atomic reduction further reduces the memory contention that could not be avoided by existing thread-local reductions. Our replicated Fock matrices approach focuses not on reducing the number of atomicAdd operations, but rather on avoiding conflicts at the same address. Note that the \(N_{\mathrm{Fock}}\) replicas are allocated in global memory, just like the main Fock matrix itself, not in thread-local registers or shared memory. Ideally, the number of conflicts from atomicAdd operations by different threads is reduced by a factor of approximately \(1/N_{\mathrm{Fock}}\), but the optimal number of replicas needs to be determined experimentally. This replicated Fock matrix approach does not require additional registers or shared memory, thus not adversely affecting the occupancy of Streaming Multiprocessors. Furthermore, it does not require complex data structures or devised thread assignment, making it a simple yet powerful acceleration technique that can be easily integrated into existing GPU implementation of the Fock matrix computation. Although additional storage space of \(O(M^{2})\times N_{\mathrm{Fock}}\) in global memory is required, the space complexity increases only by a constant factor, having minimal impact on the size of molecules for which Direct-SCF can be executed on a GPU. Figure 4 shows the additional memory usage required by the proposed replicated Fock matrices. For reference, we consider the memory sizes available on NVIDIA A100 and H100 GPUs.