Hartree Fock iteration problem

Question

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.

Answer 1

If you choose to do it one orbital at a time, the second orbital needs to be solved for the lowest energy solution in the space orthogonal to the first orbital and for the same Fock operator. You need to solve ALL the orbitals before you update your Fock operator.

You were basically solving the same orbital over and over again, and took the change of that orbital between iterations to be the second orbital.

Instead of doing it one orbital at a time, it is actually easier to solve all orbitals at the same time. Basically you are solving $ \mathbf{F}\mathbf{C}=\mathbf{S}\mathbf{C}\mathbf{\epsilon} $ every iteration.

The standard way is to perform decomposition $ \mathbf{S}=\mathbf{U}^†\mathbf{\Lambda}\mathbf{U} $ and use $\mathbf{S}^{1/2}=\mathbf{U}^†\mathbf{\Lambda}^{1/2}\mathbf{U} $ to transform your basis and Fock operator.

This way you get a more standard eigenvalue problem $ \mathbf{\tilde{F}}\mathbf{\tilde{C}}=\mathbf{\tilde{C}}\mathbf{\epsilon} $ and you can simply pull out your favorite LAPACK implementation and get the orbitals fairly quickly.

Hartree-Fock procedure might oscillate or divergence. However, it really shouldn't do so for beryllium. The HOMO LUMO gap is fairly large and you will only need s orbitals in your basis.

Answer 2

You might find the Hartree-Fock tutorial on http://vergil.chemistry.gatech.edu/resources/programming/ useful.