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, in this case, I have two Roothaan equations for two orthonormal orbitals: Fc⃗ i=Sc⃗ iϵ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)c⃗ 1=Sc⃗ 1ϵ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)c⃗ 2=Sc⃗ 2ϵ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 at each step, I have an oscillating result but the energy range obtained was not even closed to the right answer. My subprogram solving the generalized eigenvalue problem automatically generate the corresponding eigenvector and it is normalized with respect to my basis. I have posed this question in another forum but I was asked for more information. I am just hoping someone who has done it will realize what I have done wrong. Thanks in advance.