# Efficient GPU Implementation of the McMurchie-Davidson Method for Shell-Based ERI Computations
## 1 Introduction
### Related Works
Also, a year later in 2009, [PERSON] and [PERSON] [19] extended their acceleration schemes by (1) incremental computations of the Fock matrix throughout successive iterations, allowing reduction of the number of computable elements of the Fock matrix, in particular, in late iterations when convergence is near, (2) double screening of integrals: in a first stage, _ket_ and _bra_ pairs become ordered by angular momentum, and at the second stage, such pairs become ordered within each angular momentum subject to contribution to the Schwartz upper bound, enabling effective parallelization and improved load balancing, and (3) computation of the Fock matrix components by blocks and by following the ordering of the integral grid while ignoring the integrals with the Schwartz bound below a threshold. The extended algorithm achieved a performance of up to 650\(\times\) when compared to a CPU-implementation. However, the implementation in [19] is restricted to s- and p-type Gaussian basis functions and assumes that the communication time between the GPU and CPU is negligible compared to the total time. In 2012, [PERSON] et al. [20] proposed automated code generation to accelerate the MD method by GPU and CPU kernels for integral evaluation. The generated kernels are expressed in symbolic form and are subject to algebraic transformations in the Mapple CAS system. The resulting framework was implemented in TeraChem, a popular tool in quantum chemistry, rendering kernels of about 2000 lines to support up to \(d\) angular momentum. Although the approach by [PERSON] et al. [20] has the potential to explore diverse kernel configurations for specific GPU hardware, the selection of the suitable kernel follows empirical testing resembling a trial-and-error approach. Furthermore, the challenge of fully automating kernel generation and optimization arises from the limitations faced by compilers, which are often constrained by both compilation time and the availability of source code. This makes it impractical to scale them to accommodate a large angular momentum. In 2014, [PERSON] and [PERSON] [21] proposed an acceleration scheme for MD through a calculation order for electron repulsion integrals (ERIs), and a method using multiple cooperating threads to reduce register usage. Computational experiments on CPUs and GPUs showed that the approach significantly accelerated the computation by a factor of four for transition metal clusters. Although the method proposed by [PERSON] and [PERSON] [21] enables finding calculation orders that minimize the intermediate ERIs that need to be simultaneously retained, scalability to higher angular momentum is still a challenge as minimizing register spills in ERI calculations is a nonpolynomial time search problem, allowing for only approximate solutions. Whereas compilers can effectively address simpler cases, the limited scalability to high-angular-momentum ERI calculations limits the overall effectiveness. In 2017, [PERSON] et al. [22] proposed the GPU acceleration of the MD method for arbitrary angular momentum. The approach is mainly oriented towards performing integral pre-screening (integral list generation) on the CPU side by the ORCA package, and CPU to GPU memory transfers of computed lists. As such, the computational load is distributed in a way that each GPU thread calculates a single integral batch at a time. The acceleration scheme achieved speedups of up to a factor of 30 relative to serial and parallel CPU implementations. Although the approach is presented for arbitrary angular momentum, aspects such as the complexity of high angular momentum inhibit the practical usage of GPU code to f-type integrals.