I made a Hartree Fock code, in the psi4 package, following the formulas given on this page:
When I use the same geometry they use (Water with a 104 bond angle and 1.1 Å bond length), I get the exact same results as on that website. Each iteration is the same. It converges to the same results.
Here is a screenshot of the iterations for Water, which match the results on the website:
You can also notice that it matches the reference energy, which is what psi4's built-in scf solver gives. I made my own hf/scf as an excersize, but psi4 also has one built in. It has agreement with mine, so I think I've coded scf correctly.
However if I try a larger molecule, like Butadiene or Benzene, I get some incorrect results. Here is a screenshot of the output of some iterations
Notice how in this case it bounces back and forth between two results. It doesn't converge. Why is this? The built in scf gives the correct answer; it should be in the mid to low -150s which is what I got from the built in one. But mine diverges.
My question is, does scf scale badly to larger molecules? Is there some trick I need to use? Or is this indicative that I made a coding error? I've been looking for one, but I can't find any.
I also coded scf all over again, in Python, because Python is very easy to work with. I get the same results for water: It matches exactly my psi4 code, and the website, for water. But for butadiene, I get a different failure. In this case, it DOES converge, but to the wrong answer. So I get the wrong think in both psi4 and python, but different wrong things... I'm a bit confused here. I also implemented DIIS/Pulay Mixing into my python version. It converges to the same wrong answer that the "basic SCF" does, but at least it does it faster...
I think SCF is capable of doing it, but maybe I have to do something different. The built-in SCF in psi4 (that is getting the correct energy for butadiene) uses the same initial guess that I do, which is to say, it guesses that the fock matrix is the core Hamiltonian initially.