Solving Hartree Fock Iteration for Closed Shell Atoms with GTO Basis

[cheong]

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.

 

[cgk]

In SCF, you effectively do not have two different eigenvalue equations. You have just one, and take the two eigenfunctions from that which have the lowest energies. Note that a symmetric N x N matrix has not one, but N eigenvectors.

 

[cheong]

In SCF, you effectively do not have two different eigenvalue equations. You have just one, and take the two eigenfunctions from that which have the lowest energies. Note that a symmetric N x N matrix has not one, but N eigenvectors.

Thank you for the quick reply!
But it is not ture that each orbital itself has it own Fock operator? From my understanding, even in a SCF, the direct interaction terms and the exchange terms for each orbital is different, when I have more than 2 spin-orbital. And because the interaction terms are different, shouldn't I need two differen pseudo eigenvalue equations?

I have tried what you suggest too. I only got oscillating behavior. Perhaps my code is fundamentally wrong.