Computational quantum chemistry employs quantum mechanics to investigate the electronic structure of molecules and atoms. While the Schrödinger equation forms the foundation of this approach, its application to molecules with multiple nuclei and electrons poses significant computational challenges. To address this complexity, approximate methods like the Hartree-Fock method are commonly used. An integral part of Hartree-Fock calculations is the two-electron repulsion integral (ERI), computed for each combination of basis functions representing electron orbitals. Although several algorithms exist for ERI calculations, including the widely used McMurchie-Davidson (MD) algorithm, their implementation on GPUs is constrained by recursive formulas. This paper proposes an efficient parallel MD algorithm for GPUs and presents its implementation. Specifically, we divide the work into batches that can be processed concurrently and recursively, enabling efficient MD algorithm implementation on GPUs. By utilizing batches, parallel computation becomes feasible even with the small amount of shared memory of GPUs. Moreover, our approach allows computation with a single function without the need for separate functions for each combination of azimuthal quantum numbers representing electron orbitals. Experimental results show that the proposed GPU implementation can perform the ERI computation up to 101 and 24 times faster than the CPU implementation on an AMD EPYC 7702 CPU for monatomic and polyatomic molecules, respectively.