# Efficient GPU Implementation of the McMurchie-Davidson Method for Shell-Based ERI Computations
## 1 Introduction  
Quantum chemistry has been the cornerstone in offering the building blocks and the first principles to derive the molecular and physical properties arising from (sub)atomic interactions. As such, the field has been able to predict and explain the formation of electronic structures and chemical reactions in science, engineering, and medical applicationsubiquitously. One of the fundamental challenges in the realization of practical quantum chemistry lies in solving the Schrodinger equation [1] to compute the electronic orbitals representing molecular properties:  
\[\hat{H}\Psi=E\Psi, \tag{1}\]  
where \(\hat{H}\) represents the total energy of the molecular system (kinetic and potential energy, Hamiltonian), \(\Psi\) denotes the wave function encoding the quantum state of the electron within the molecule (electron kinematics), and \(E\) is the energy eigenvalue. The central aim of quantum chemistry is to find solutions to atomic and molecular Schrodinger equations. Yet, as molecules of practical interest become largely polyatomic, suggesting a large number of nuclei and electrons, solving the Schrodinger equations analytically becomes intractable due to the large number of (sub)atomic interactions (the many-body problem). Using reasonable assumptions on electron kinematics, approximate schemes, such as the Hartree-Fock method [2; 3], have become attractive to estimate energy states and the overall electronic structures [4]. The Hamiltonian operator of an \(N\) electron system is as follows:  
\[\hat{H}=\sum_{i=1}^{N}\hat{h}(i)+\sum_{i>j}^{N}\hat{g}(i,j) \tag{2}\]  
where \(i,j\) are electrons, \(\hat{h}(i)\) is a one-electron operator representing the kinetic energy of electrons, and electron-nucleus Coulomb attraction, and \(\hat{g}(i,j)\) is the two-electron operator representing the electron repulsion. The minimization of the expectation conditions of \(\hat{H}\) leads to the Hartree-Fock equations:  
\[\left[\hat{h}(i)+\sum_{k}\int\hat{g}(i,j)|\psi_{k}(j)|^{2}\,d\tau_{j}\right] \psi_{i}(i)-\sum_{k}\left[\int\hat{g}(i,j)\psi_{k}^{*}(j)\psi_{i}(j)\,d\tau_{j }\right]\psi_{k}(i)=\varepsilon_{i}\psi_{i}(i)\quad(i=1,2,\ldots,n) \tag{3}\]  
Solutions to the above (3) render the orbital functions \(\{\psi_{i}\}\) and the orbital energies \(\{\varepsilon_{i}\}\). To tackle the above, the Self-Consistent Field (SCF), introduced by [PERSON] [2; 5; 6], is often used to estimate the orbital functions \(\{\psi_{i}\}\) iteratively. By assuming an approximate initial guess of \(\{\psi_{i}\}\), replacing of \(\{\psi_{k}\}\) in (3) leads to the following:  
\[\hat{F}\psi=\varepsilon_{i}\psi \tag{4}\]  
in which the first solutions are computable from a normal eigenvalue problem. The first solutions \(\{\psi_{i}\}\) and \(\{\varepsilon_{i}\}\) are expected to be better than the initial guess. Afterward, a subsequent set of solutions is obtained after computing the operator \(\hat{F}\) using the previous solutions. In the above, initial guesses are used to compute the potential field felt by an electron. As such, the Hartree-Fock method assumes the single-electron configuration, in which electrons are subject to an average neighborhood electron potential. A relevant procedure is evaluating the potential encoded through a corresponding basis set and integrals. Often, solutions of the orbital functions \(\{\psi_{i}\}\) at each step of the SCF procedure are expressed as a linear expansion of \(M\) basis functions:  
\[\psi_{i}=\sum_{k}^{M}c_{ik}\phi_{k} \tag{5}\]  
where \(M\) is the number of basis functions, and \(c_{ik}\) is the expansion coefficient. Furthermore, the Hartree-Fock method solves the set of integrals, often termed the _two-electron repulsion integrals_ (Basis-ERIs), encoding the Coulombic repulsion interactions between electron pairs across different orbitals.