# GPU-Accelerated Fock Matrix Computation with Efficient Reduction
## 2 Background
### Fock Matrix Computation
#### 2.2.2 Update of the Fock Matrix  
We restate the definition of the Fock matrix for \(M\) basis functions,  
\[F_{\mu\
u}=h_{\mu\
u}+\sum_{\lambda=1}^{M}\sum_{\sigma=1}^{M}D_{ \lambda\sigma}(2(\mu\
u|\lambda\sigma)-(\mu\lambda|\
u\sigma))=h_{\mu\
u}+2J_{ \mu\
u}-K_{\mu\
u}, \tag{21}\] \[J_{\mu\
u}=\sum_{\lambda=1}^{M}\sum_{\sigma=1}^{M}D_{\lambda \sigma}(\mu\
u|\lambda\sigma),\] (22) \[K_{\mu\
u}=\sum_{\lambda=1}^{M}\sum_{\sigma=1}^{M}D_{\lambda \sigma}(\mu\lambda|\
u\sigma). \tag{23}\]  
During the SCF procedure, while the molecular integrals \(h_{\mu\
u}\) and \((\mu\
u|\lambda\sigma)\) remain constant, the Fock matrix \(\mathbf{F}\) is updated in each iteration according to the density matrix \(\mathbf{D}\). Equation (21) is computationally demanding not only due to the evaluation of ERIs itself, but also because of the reduction of \(O(M^{4})\) ERIs into an \(O(M^{2})\) memory space of the Fock matrix. The efficient implementation of the ERI reduction, particularly using GPUs, is challenging and has been a subject of ongoing discussion in previous research over several years.  
In Direct-SCF, where ERIs are computed on-the-fly, an update equation of the Fock matrix is widely used rather than the definition itself. Equation (21) can be straightforwardly implemented in Stored-SCF, where all ERI values are stored in memory. However, in Direct-SCF, where ERIs are evaluated on-the-fly, redundant computations frequently occur, such as the overlap between \((\mu\
u|\lambda\sigma)\) for the Coulomb \(J\)-matrix and \((\mu\lambda|\
u\sigma)\) for the exchange \(K\)-matrix. Therefore, an update scheme for the Fock matrix that eliminates redundant integral calculations has been devised by considering the eight-fold symmetry of ERIs. Applying Equation (18) to Equation (21), the contribution of each unique ERI \((\mu\
u|\lambda\sigma)\) to the Fock matrix can be calculated as  
\[F_{\mu\
u}=F_{\mu\
u}+4D_{\lambda\sigma}(\mu\
u|\lambda\sigma), \tag{24}\] \[F_{\lambda\sigma}=F_{\lambda\sigma}+4D_{\mu\
u}(\mu\
u|\lambda \sigma),\] \[F_{\mu\lambda}=F_{\mu\lambda}-D_{\
u\
u}(\mu\
u|\lambda\sigma),\] \[F_{\
u\sigma}=F_{\
u\sigma}-D_{\mu\lambda}(\mu\
u|\lambda\sigma),\] \[F_{\mu\sigma}=F_{\mu\sigma}-D_{\
u\lambda}(\mu\
u|\lambda\sigma),\] \[F_{\
u\lambda}=F_{\
u\lambda}-D_{\mu\sigma}(\mu\
u|\lambda\sigma).\]  
Equation (24) can be used to implement the Coulomb and exchange terms of Equation (21) by accumulating each of the \(M(M+1)(M^{2}+M+2)/8\) unique ERIs into six elements in the Fock matrix. While this update scheme can eliminate redundant integral calculations, \(O(M^{4})\) ERIs are added directly to a small memory space of \(O(M^{2})\). Whenthis reduction operation is performed in parallel on a GPU, high memory contention arises from multiple threads simultaneously attempting to exclusively add integral values to the same elements of the Fock matrix, resulting in significant overhead. Consequently, even if it involves redundant ERI calculations, the strategy of constructing the Fock matrix by computing the \(\mathbf{J}\)- and \(\mathbf{K}\)-matrices separately is also being explored [12, 24, 27]. In this work, we provide an efficient GPU implementation that addresses the reduction issue of the memory contention in Equation (24), while fully exploiting the eight-fold symmetry of ERIs.