# Efficient GPU Implementation of the McMurchie-Davidson Method for Shell-Based ERI Computations
## 1 Introduction
### Related Works  
For systems with many electrons, the Hartree-Fock method becomes computationally expensive due to the challenge of solving a large number of integrals, e.g., Basis-ERIs require the computation of every quartet of basis functions representing the molecule of interest. In order to realize the efficient computation of electron repulsion integrals over Gaussian functions (GTO-ERIs), various algorithms have been proposed. Two key algorithms in the field are the well-known McMurchie-Davidson method (MD) [10] and the Obara-Saika method (OS) [11]. Inspired by the idea of using Cartesian Gaussian Functions [8], the MD method allows computing of the molecular integrals through Hermite Gaussians, expressed in compact form by auxiliary functions, and recursive relations that enable calculation of the higher angular momentum efficiently. Also, the MD method, due to its relatively simple mathematical structure, offers greater generality and extensibility compared to other approaches. The OS method uses the Gaussian Product theorem to render recurrence equations to simplify composite integrals. Several other efficient schemes exist, such as the Pople-Hehre (PH) method [12] using normal specification of basis functions sharing common computable information, allowing coordinate components to become zero or constant and decreasing computational load in the corresponding integrals, the Dupuis-Rys-King (DRK) method [13] using root search of (orthogonal) Rys polynomials, often termed the Rys quadrature method [14], yielding representations for basis functions of arbitrarily high angular momentum. Another approach is the Head-Gordon-Pople method (HGP) [15] that (1) extends the recurrence relations from OS to consider the shifting angular momentum, and (2) benefits from efficient evaluation of integrals by evaluating basis functions sharing the same information (coefficients and exponents), by evaluating similar integrals simultaneously, and by grouping similar atomic orbitals. The PH method is particularly suitable for calculations involving high-angular-momentum systems, whereas the HGP method is more efficient for systems with low-angular-momentum functions. Also, the PRISM method [16] capitalizes on the trichotomy between angular momentum and degree of contraction in ERIs to find tailored computation paths. As such, PRISM allows performing of contraction steps precisely where needed, and realizes a degree of flexibility to cover a wide range of ERI. The reader may find a further detailed description of the above in [17]. Computing the electronic structure configurations through Hartree-Fock methods is amenable to parallelization and acceleration schemes by Graphics Processing Units (GPUs). As such, the active developments of new GPU architectures and the concomitant Compute Unified Device Architecture (CUDA) have enabled several Hartree-Fock estimation schemes to capitalize on the merits of GPU hardware. The community has rendered several attempts on GPU acceleration of the MD method [18; 19; 20; 21; 22; 23]. In 2008, [PERSON] and [PERSON] [18] were among the first groups in accelerating the MD algorithm by tailored forms of stream processing on GPUs, and in evaluating specific mappings between threads and integral computation, e.g., one thread per contracted integral and one thread per primitive integral. The resulting algorithm achieved a 130\(\times\) speedup over CPU implementation.