Why does my Hartree Fock/SCF code work for small molecules (water, methane etc) but goes crazy for slightly larger ones (Benzene, Butadiene etc)? I made a Hartree Fock code, in the psi4 package, following the formulas given on this page:
http://sirius.chem.vt.edu/wiki/doku.php?id=crawdad:programming:project3
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. 
 A: There is nothing inherent in SCF that makes it lose efficacy as you increase the number of electrons. Certainly, coding errors are always possible so I would first check that the H2O case is absolutely correct. It seems to be so, however try tightening the SCF convergence tolerance in both your program and Psi4 and make sure that you reproduce each converged digit in both total HF energy and molecular orbital energies properly (i.e. convergence of 10^-6 on energy means that energy should match to 6 decimal places). Next, obviously, make sure your one- and two-body atomic orbital integrals are exactly the same what Psi4 uses. You should also check all the convergence parameters like DIIS size. Another thing to consider is how d-orbitals (and higher) are handled in your basis set. You should keep in mind that if your atomic orbital basis contains d-orbitals, there will be orbital degeneracies if you represent them using Cartesian Gaussians. This can mess up SCF when you try to orthonormalize orbitals to solve the pseudo-eigenvalue problem. In commercial codes, AO integrals are calculated using Cartesian Gaussians and transformed to real spherical Gaussians before any calculations are done (these are simple linear transformations can be found online or worked out using spherical basis vectors). Up to p-orbitals, Cartesian and spherical Gaussians are the same, so you may have missed this problem with just H2O. Finally, another culprit is that (I assume) you are not incorporating Abelian point group symmetry into your SCF solver. All commercial codes do this, not just for computational savings, but also because it can really clean up your orbital optimization space and prevent you from falling into certain local minima, which is likely what's happening given the oscillating output of your code. To check this, try running Psi4 with C1 symmetry and see where it lands. 
