# Unrestricted Hartree-Fock : Density matrices initialization

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:

• My understanding is that the UHF dissociation curve for $\ce{H2}$ should always be lower than the RHF curve, so that your putative UHF curve starts out higher in energy is already suspicious. I imagine random MO coefficients would be even worse than having them be all zero, since you could get some very odd contributions from the electron interaction terms. There is a related question/answer on Chem SE
– Tyberius
Commented Oct 15, 2023 at 21:47
• +1, the RHF and UHF curves look reasonable, and the points of the red line that do not fall onto the green line also look reasonable. But the green line appears unphysical to me. You may want to check whether your program has a bug, especially when the initial density matrix is not idempotent and when the spin symmetry is broken. Commented Oct 17, 2023 at 9:07

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.