I am writing a program to compute the ground state energy for any closed shell atom using Hartree Fock Roothaan method, with GTO basis. The code works for the simplest case, the helium, but it fails with beryllium (z=4).
I understand that I have two Roothaan equations for two orthonormal orbitals:
$$F\vec{c}_i=S\vec{c}_i\epsilon_{i}$$
As usual, I have initial guess for the coefficient matrix. Use it to generate the Fock matrix and find the coefficient for the first orbital by solving a generalized eigenvalue problem, choosing only the smallest energy.
$$F^{(0)}\vec{c}_1=S\vec{c}_1\epsilon_{1}$$
Than I use the coefficient of the new 1_st orbital (the 2_nd orbital is not changed yet) to generate a new Fock matrix to find the coefficient of the 2_nd orbital
$$F^{(1)}\vec{c}_2=S\vec{c}_2\epsilon_{2}$$
And use the new 2_nd orbital to find a new 1_st orbital again until energy converges.
My problem is that with that algorithm I always generate orbitals with the same coefficients. But different orbitals are supposedly orthonormal to each other.
If I impose orthonormal condition by gram-schmidting the coefficient vectors, I have an oscillating result but the energy range obtained was not even closed to the right answer.
My code to solve the generalized eigenvalue problem automatically generate the corresponding eigenvector and it is normalized with respect to my basis.
I really hope some knowledgeable people can point me to the right way to do it.