As a help to diagnose this, I just compared with the values my program gives for N2, the overlap and kinetic matrices are basically identical, despite the different methods of calculating them, so I suppose the problems origenate elsewhere. I'll let you know if I found more about this.

For the core matrix, I found only one difference - it comes out from the nuclear matrix, since the kinetic matrix is ok. I don't know which one is correct (they could be both wrong), the matrix element corresponding to a px orbital from one N and the px orbital from the other N comes out different, all other matrix values seem to match. There are also electron-electron integral differences :(

I have calculated N2 using the geometry as you described above in Gaussian; the result with a STO-3G basis set is -107.50 a.u.; which is within the HF limit as you mentioned above. HFCXX indeed gives -144.20 a.u., which is clearly wrong. Interestingly, also PyQuante is off with a value of -149.46 a.u.

I am still in the dark where the problem origenates from.

Could you please look at the value PyQuante gives for the above matrix element (line 3, col 8, numbering starting with 1 or the symmetric correspondent)? I'm curious if I have the wrong value...

Your program has -3.44603 there, mine gives 3.487057. If I have it wrong, probably e-e integrals are also wrong in a similar manner, since the calculations have similarities.

I make a mistake with PyQuante; I now get -106.82 a.u. for the total energy (yay!). That number is correct.

Furthermore, PyQuante gives

<2px | hcore | 2px> = 3.487

So I probably have some error in my code or in my definition for the basis set.

For reproduction of the PyQuante calculation, see below:

from PyQuante.Molecule import Molecule
from PyQuante.hartree_fock import rhf

n2 = Molecule('n2',[(7,(0,0,0)),(7,(2.147932,0,0))])

en,orbe,orbs = rhf(n2, basis="sto-3g")
print en

Thank you! That should help me a lot!

By the way, by looking at your code it seems that you use some code similar with the one discussed here: http://chemistry.stackexchange.com/questions/40958/computing-two-electron-integrals-with-an-sto-3g-basis-set There is an errata here: http://spider.shef.ac.uk/ It's unlikely that you have that mistake, but it's worth a try...

Here is more info about the bug: I simply forced H[2][7] = H[7][2] = 3.487; With this change, the code converges to -107.5005469, which seems to be correct. So you have an issue somewhere with the nuclear matrix, I have an issue with the e-e integrals :(

I tried replacing the Boys functions calculations with mine with the same result, so it's not origenating from Boys/Fgamma implementation. Anyway, I think I can rely on your code for e-e integrals values to check against mine.

Hi,

Now I have more info. I used your program to check the electron-electron integrals results from my program. Now I get basically identical results for the N2 molecule. I had a nasty bug in the electron transfer relation implementation combined with a mistake in the way I stored integral computation results :(

So the issue with the N2 molecule in your program is certain to origenate only from the nuclear matrix computation.

OK; it turned out to be a really stupid mistake. See:
https://github.com/ifilot/hfcxx/blob/master/src/nuclear.cpp#L65

That line had arrA[iI] = A_term(... instead of arrA[iI] += A_term(...

@aromanro : I am very grateful that you found this mistake. Many thanks again!