# GPU-Accelerated Fock Matrix Computation with Efficient Reduction
## 2 Background
### Fock Matrix Computation
#### 2.2.1 Two-Electron Repulsion Integrals  
Two-electron repulsion integrals (ERIs) describe the Coulomb interactions between two electrons in a given chemical system and are defined as a double integral over basis functions:  
\[(\mu\
u|\lambda\sigma)=\iint\chi_{\mu}(\mathbf{r}_{1})\chi_{\
u}(\mathbf{r}_{1})\frac{ 1}{r_{12}}\chi_{\lambda}(\mathbf{r}_{2})\chi_{\sigma}(\mathbf{r}_{2})d\mathbf{r}_{1}d\mathbf{r} _{2}, \tag{11}\]  
where \(r_{12}\) denotes the Euclidean distance between electron 1 and 2. The basis functions \(\chi(\mathbf{r})\) are trial functions for atomic orbitals at the 3-dimensional coordinates \(\mathbf{r}=(x,y,z)\), representing electron orbitals around an atom. The chemical notation \((\mu\
u|\lambda\sigma)\) denotes an ERI over four basis functions \(\chi_{\mu}\), \(\chi_{\
u}\), \(\chi_{\lambda}\), and \(\chi_{\sigma}\) (\(1\leq\mu,\
u,\lambda,\sigma\leq M\)), corresponding to a single integral value. To simplify the calculations of molecular integrals, the basis functions are typically given as linear combinations of Gauss-type orbitals,  
\[\chi_{\mu}(\mathbf{r})=\sum_{a=1}^{K_{\mu}}w_{a}G_{lmn}(\mathbf{r},\mathbf{A},\rho_{a}), \tag{12}\]  
where \(K_{\mu}\) denotes the number of Gauss-type orbitals and \(w_{a}\) are the weight coefficients for the linear combination. The Gauss-type orbital is defined as  
\[G_{lmn}(\mathbf{r},\mathbf{A},\rho_{a})=c\cdot(x-A_{x})^{l}(y-A_{y})^{m}(z-A_{z})^{n} \text{exp}(-\rho_{a}|\mathbf{r}-\mathbf{A}|^{2}), \tag{13}\]  
where \(c\) is a normalization constant, \(\mathbf{A}=(A_{x},A_{y},A_{z})\) denotes the 3-dimensional coordinates of a nucleus, \(l,m,n\) are non-negative integers representing angular momenta, and \(\rho_{a}\) is the exponent that determines the radial extent of the orbital. The orbital shape is determined by the total angular momentum \(L=l+m+n\) and named such as s-, p-, and d-orbital, corresponding to \(L=0,1,2\), respectively. The set of parameters \(K\), \(w_{a}\), \(L\), and \(\rho_{a}\) is known as a basis set. Various basis sets have been proposed and are available in databases such as Basis Set Exchange [30]. In quantum chemical calculations, an appropriate basis set is selected based on the target system and application, with common examples including the Pople basis sets (e.g., STO-NG) and the Dunning basis sets (e.g., cc-pVXZ). By expanding Equation (11) with Equation (12), the ERI over basis functions can be calculated as a weighted sum of ERIs over Gauss-type orbitals,  
\[(\mu\
u|\lambda\sigma)=\sum_{a=1}^{K_{\mu}}\sum_{b=1}^{K_{\
u}}\sum_{c=1}^{K_ {\lambda}}\sum_{d=1}^{K_{\
u}}w_{a}w_{b}w_{c}w_{d}[ab|cd]. \tag{14}\]In contrast to the contracted ERI \((\mu\
u|\lambda\sigma)\), \([ab|cd]\) is referred to as a primitive ERI, similarly defined by four Gauss-type orbitals,  
\[[ab|cd]=\iint G_{a}(\mathbf{r}_{1})G_{b}(\mathbf{r}_{1})\frac{1}{r_{12}}G_{c}(\mathbf{r}_{2} )G_{d}(\mathbf{r}_{2})d\mathbf{r}_{1}d\mathbf{r}_{2}, \tag{15}\]  
where \(G_{a}(\mathbf{r})\) is a shortened form of the Gauss-type orbital \(G_{lmn}(\mathbf{r},A,\rho_{a})\). The primitive ERI \([ab|cd]\) is often calculated using recurrence relations such as the McMurchie-Davidson method [31], Obara-Saika method [32], and Head-Gordon-Pople method [33]. In either recursive function, Equation (15) is known to be expandable into a linear combination of a generalized mathematical function called the Boys function [34]:  
\[[ab|cd]=\sum_{j=0}^{L_{a}+L_{b}+L_{c}+L_{d}}c_{j}F_{j}(X), \tag{16}\]  
where the expansion coefficients \(c_{j}\) and \(X\) are determined by the orbital center \(\mathbf{A}\) and exponent \(\rho_{a}\), and the total angular momentum \(L\). The Boys function is given as an integral function as follows:  
\[F_{j}(X)=\int_{0}^{1}t^{2j}\text{exp}(-Xt^{2})dt, \tag{17}\]  
and several numerical evaluation methods have been devised for its efficient computation [21; 35]. In this work, we compute the primitive ERI \([ab|cd]\) using kernels expanded to the form of Equation (16) through the Head-Gordon-Pople method [25].