Two-electron integral algorithm
There are many things to consider and improve your code. Many components
in an ab initio electronic structure calculations need supervision to achieve a
good performance. I am no expert on those calculations, but I will point you to some ideas and references. On the theoretical side, you could calculate E_ab for all shell pairs first. For more clarification on the order of operations and their scaling look at
Figure 3 in Multi-electron integrals. The article is written by a well-recognized researcher in the field that wrote a monumental book where you can find all needed expressions.
On the computational side, I would say that maybe Numba is not able to optimize your code. Writing it with Numba is not a guarantee of speed. I just took a glimpse at the code and, without going deep in the code, I saw many if statements. Only a straightforward observation. As general advice, it is all about how data is moved
and the flow of operations. Avoid allocations, try to rewrite as much as possible.
Profile your code to find the bottlenecks. Try to write it as vector or matrix operations to be able to use the large registers of the processor. There are many tricks to do these operations fast, but sadly they are buried elsewhere.
I would like to end with a reference to the recent literature on the topic. I just came up with an article that I managed to find again now: Fast Evaluation of Two-Center Integrals over Gaussian Charge Distributions and Gaussian Orbitals with General Interaction Kernels.