Unrestricted Hartree-Fock : Density matrices initialization

Question

I have some results about the energy of H2 as a function of the bond length using Restricted HF (RHF) and Unrestricted HF (UHF) methods. With a zero-initialization of the density matrix I have the same results which is expected. But When I use a random initialization (Uniform on [0,1]) for both density matrices, the scf loop converges to weird values. And when I only initialize one of the density matrices I have what I suppose to be the UHF solution.

How are density matrices initialized in Unrestricted calculations (UHF, Unrestricted Kohn-Sham)? Is there any schemes to have a more robust algorithm?

legend:

• uhf = UHF with zero initialization of both matrices
• uhf P_alpha random = P_alpha is initialized with uniform distribution on [0, 1] and P_beta as zeros
• uhf P_alpha and P_beta random = both matrices are initialized with uniform distribution

Edit:

I forgot the Coulomb interaction between the alpha and beta orbitals and I used rotation mentionned in Tyberius comment.

As initial guess I used the HOMO and LUMO orbitals of H2:

$$C_{\alpha} = C_{\beta} = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1\\ -1 & 1 \end{pmatrix}$$

Then I rotated the $\alpha$ orbitals, with $\theta = \frac{\pi}{2}$ : $$ C^{new}_{\alpha} = \begin{pmatrix} cos(\theta) & -sin(\theta)\\ sin(\theta) & cos(\theta) \end{pmatrix}C_{\alpha}$$ Here are the results:

Answer 1

When you dissociate the singlet state of H$_2$, you need the $\alpha$ and $\beta$ spatial orbitals to be different. If you use the same initial guess for both spin channels, as what happens when you set both initially to zero (this actually corresponds to the core guess i.e. diagonalization of the one-electron part of the Hamiltonian). In your random guess you clearly see that you have obtained at least two solutions to the SCF equations, which means that the higher-lying one has to be a saddle point solution. The trick is achieving this in a consistent mode. What typically works is starting from a large $R$, and reading in these orbitals as the initial guess for the calculation for the next largest value of $R$ (repeat ad nauseam).

The addendum with the rotated guess looks fine; at separation the energy goes to two times the energy of a single atom. You don't get $-1 E_h$ exactly, since you have basis set truncation error. If you plot the energy relative to that of two isolated atoms, you will see that this interaction energy goes to zero at large distance.