Felix, thank you very much for the answer and data. It’s really useful.
When I posted the other day I had a large difference in energy. Finally I found out that the problem was not in the ERIs but on the way I was using the 8-fold permutational symmetry of the ERIs to do only the unique integrals. I constructed the ERI(a,b,c,d) matrix using this symmetry, and when I switched back to doing the full matrix, there were about 4 elements (in 17^4) that had some deviations. I still don’t know why is that, but started doing all the integrals, and then my energy was almost correct. I was about to post an answer here, but noticed that instead of -40.18… Ha for 6-31g methane obtained with Orca, I got -40.20… Ha. This small difference I still can’t figure out. So there is still some problem.
For the ERIs I’m using Obara-Saika derived methods, something similar to Head-Gordon Pople. I use the vertical recurrence to obtain the [L s s s] integrals (up to [g s s s] integrals with one algorithm and up to [L = 10 s s s] with another algorithm, slightly slower but allowing arbitrary angular momentum). I use the electron transfer recursion to obtain [La s Lc s] integrals, then I contract and then use horizontal recursion to obtain the full integral matrices for each shell. The code is vectorized, and works on all cartesian components simultaneously.
I do permutations to get all shells into the order I just described, and I think the permutations are correct.
Once again, thank you very much,
I’ll post when/if I solve the problem.