I have been trying to write a basic HF code that can compute the MOs and total energy of a molecule. I followed the following reference for overlap and kinetic integrals. During a conversation with a developer of Libint, I was told that I should try implementing spherical harmonics into the code. I found this paper for the conversion from cartesian to spherical harmonics. Any beginner or intermediate level pointers on how the integrals are computed in spherical harmonics would be much appreciated. I am not looking for specific implementation details, but for the general algorithm.
Molecular integrals aren't typically computed in spherical harmonics. What you do instead is to compute the integrals in the Cartesian representation, where they are separable, and then transform to the spherical representation as quickly as possible your algorithm allows, since there are $n_c=\lambda(\lambda+1)/2$ cartesians but only $n_s=2\lambda+1$ spherical harmonics for angular momentum $\lambda$.
What this means is that when computing the $(pq|rs)$ integrals, you convert an index to the spherical basis as soon as possible, so that the recursions for the remaining indices only deal with $n_s$ components for each transformed index instead of $n_c$ components for untransformed indices.
In practice, it is hard to say at which point the transformation should be done, since this is where all algorithms are different: McMurchie-Davidson, Obara-Saika / Head-Gordon-Pople, and Rys quadrature all have different recursion relations and different scaling with the total angular momentum.
For a concise introduction to molecular integrals, I recommend you consult "the purple Bible" i.e. Molecular Electronic-Structure Theory by Helgaker, Jørgensen, and Olsen, which discusses the various algorithms at length. Its chapter 9 is dedicated to molecular integral evaluation and is almost a 100 pages!