# Initial guess for Unrestricted Hartree-Fock calculation

I am implementing an Unrestricted Hartree-Fock (UHF) calculation, following Ref. [1]. It works fine for a odd number of electrons ($N_\alpha\neq N_\beta$) since, despite the same initial guess for $\alpha$ and $\beta$ spins, the density matrices at the second step will be different and I get different orbital energies for different spin orbitals. The problem arises when I have an even number of electrons ($N_\alpha=N_\beta$) . Authors of Ref. [1] state that (p. 216)

When $N_\alpha=N_\beta$, a restricted solution of Roothaan equations is a solution to the unrestricted Pople-Nesbet equations. This restricted solution always exists and necessarily results if an initial guess $P_\alpha=P_\beta$ is used. [...] In seeking this second solution [the unrestricted one], it is imperative that an initial guess $P_\alpha \neq P_\beta$ be used or the equations will necessarily yield the restricted solution.

Unfortunately there is no clue on which initial guess should be taken. In fact in this book, for the restricted case, the Fock matrix is initially approximated by the core Hamiltonian (Guess=Core in Gaussian09), which consists in taking $P_\alpha=P_\beta=0$.

Since the core Hamiltonian guess cannot be differentiated for $\alpha$ and $\beta$ spins, I implemented Huckel guess (Guess=Huckel in Gaussian09) following Ref. [2] which consists in approximating the initial Fock matrix with the following: $$F_{ij} = KS_{ij}\frac{H_{ii}+H_{jj}}{2}$$ $H$ is the core Hamiltonian, $S$ is the overlap matrix and $K$ is a constant. I tried to slightly change the coefficient $K$ (the suggested value is 1.75) for different spins in order to get a different initial guess. However, the SCF procedure converges towards a value that is higher in energy than the restricted solution. This should not be the case, in particular for a dissociating $\ce{H2}$ molecule (i.e. distance $R=4.4$ Bohr).

How I can properly differentiate the initial guess for $\alpha$ and $\beta$ spins?

EDIT

The reference of Gaussian09 states the following:

Requests that the HOMO and LUMO be mixed so as to destroy α-β and spatial symmetries. This is useful in producing UHF wavefunctions for singlet states. [...]

I belive that this is the solution to my problem, but what does actually means HOMO and LUMO be mixed in practice (i.e. before starting the calculation)?

[1] Szabo and Ostlund, Modern Quantum Chemistry, Dover, 1996.

[2] M. Lee et al., Extended Hückel and Slater’s rule initial guess for real space grid-based density functional theory, Computational and Theoretical Chemistry, 2015.

• It only means spatial distribution of alpha and beta can be different, which is easiest to set up by mixing homo lumo levels – Greg Dec 13 '15 at 13:04
• @Greg I tried to differentiate the spatial distribution by choosing a different coefficient for the Hückel guess, but apparently it does not work. In order to mix the orbitals I need to perform at least one SCF step (i.e. I need to compute HOMO and LUMO orbitals for my initial guess of Fock matrix). But I didn't find anywhere how to mix these orbitals in practice afterwards. – user23061 Dec 13 '15 at 13:32
• Generally you use a new HOMO = constant (old HOMO + old LUMO); new LUMO = constant (old HOMO - old LUMO). The trick is that you actually should know what you are doing, and choose a related pair of orbitals (i.e. sometimes you have to mix HOMO-1 or HOMO-2 .. instead of HOMO...). Note, it does not necessarily work out of the box. Present toolbox of QChem is more about nice behaving closed shell molecules, open shell situations are often neglected as unimportant / difficult. – Greg Dec 13 '15 at 22:00

When $N_{\alpha} = N_{\beta}$, a restricted solution of Roothaan equations is a solution to the unrestricted Pople-Nesbet equations. This restricted solution always exists and necessarily results if an initial guess $P_{\alpha} = P_{\beta}$ is used. [...] In seeking this second solution [the unrestricted one], it is imperative that an initial guess $P_{\alpha} ≠ P_{\beta}$ be used or the equations will necessarily yield the restricted solution.

The quote you state from Szabo and Ostlund is not wrong, but isn't a practical quote in the sense that it helps to think about the different spin cases that one would want to handle. You already know how RHF/closed-shell calculations work in theory and in practice. I won't talk about restricted open-shell Hartree-Fock (ROHF), because it isn't relevant to this problem.

There are three UHF/open-shell scenarios, which will be considered in order:

• Odd number of electrons
• Even number of electrons, high spin
• Even number of electrons, low spin

## Odd number of electrons

For this case, take the copper(II) atom, which has a 3d$^{9}$ valence electron configuration; it is a doublet, so it has 1 more $\alpha$ electron than $\beta$ electron. There's no trouble here, because the indices for occupied orbitals run over different ranges; in the simplest case, the range for $\alpha$ spin will be 1 larger than for $\beta$ spin.

## Even number of electrons, high spin

For this case, let's take triplet oxygen. In practice, this is actually identical to the above case, except we have 2 more $\alpha$ electrons than $\beta$ electrons. Again, we have differently-sized occupied blocks for each set of MO coefficients which will lead to different densities.

## Even number of electrons, low spin

For this case, let's take singlet oxygen.

This is the difficult case, because the orbital degeneracy means for a proper qualitative description of the molecular orbitals, two determinants must be used, requiring a CAS(2,2) (at a minimum). Because this valence electron degeneracy isn't possible in a single determinant picture, the degenerate MOs will "split" into one lower in energy, the other higher. Requesting a singlet spin multiplicity means the two electrons will now occupy the same orbital; there's no reason for them to occupy different orbitals, at least for the initial guess. Hopefully it's (visually) clear that if the spatial components of $\alpha$ and $\beta$ spins are the same in this case, our guess is equivalent to the RHF guess. This means that for an open-shell singlet, one must break spatial symmetry in the initial guess, just as Szabo and Ostlund states. Since the RHF solution always exists as a local minimum on the "determinental" potential energy surface, whatever you choose as your SCF converger is most likely going to get stuck in this minimum.

This is a little hand-wavy, since we know that UHF calculations approximate two real determinants at the cost of spin contamination, but that isn't our goal.

# Breaking spatial symmetry

How do we break spatial symmetry? Consider the initial steps of a typical RHF calculation that diagonalizes the one-electron (core) Hamiltonian:

\begin{align} \mathbf{H}^{\text{core}} &= \mathbf{T}_{\text{e}} + \mathbf{V}_{\text{ne}} \\ \mathbf{F}^{0'} &= \tilde{\mathbf{S}}^{-1/2} \mathbf{H}^{\text{core}} \mathbf{S}^{-1/2} \\ \mathbf{F}^{0'} \mathbf{C}^{0'} &= \mathbf{\epsilon}^{0'} \mathbf{C}^{0'} \\ \mathbf{C}^{0} &= \mathbf{S}^{-1/2} \mathbf{C}^{0'} \\ P_{\mu\nu} &= \sum_{i}^{\text{occ}} C_{\mu i} C_{\nu i} \end{align}

where the tick mark signifies a quantity is in the orthogonal AO basis, as opposed to the original (non-orthogonal) AO basis. Now assume we're doing a UHF calculation:

\begin{align} \mathbf{H}^{\text{core}} &= \mathbf{T}_{\text{e}} + \mathbf{V}_{\text{ne}} \\ \mathbf{F}^{0'} &= \tilde{\mathbf{S}}^{-1/2} \mathbf{H}^{\text{core}} \mathbf{S}^{-1/2} \\ \mathbf{F}^{0'} \mathbf{C}^{\alpha, 0'} &= \mathbf{\epsilon}^{\alpha, 0'} \mathbf{C}^{\alpha, 0'} \\ \mathbf{F}^{0'} \mathbf{C}^{\beta, 0'} &= \mathbf{\epsilon}^{\beta, 0'} \mathbf{C}^{\beta, 0'} \\ \mathbf{C}^{\alpha, 0} &= \mathbf{S}^{-1/2} \mathbf{C}^{\alpha, 0'} \\ \mathbf{C}^{\beta, 0} &= \mathbf{S}^{-1/2} \mathbf{C}^{\beta, 0'} \\ P_{\mu\nu}^{\alpha} &= \sum_{i}^{\text{occ }\alpha} C_{\mu i}^{\alpha} C_{\nu i}^{\alpha} \\ P_{\mu\nu}^{\beta} &= \sum_{i}^{\text{occ }\beta} C_{\mu i}^{\beta} C_{\nu i}^{\beta} \end{align}

This is misleading! Because we only have one initial Fock matrix, it will only have one set of eigenvectors, so $\mathbf{C}^{\alpha, 0'} \equiv \mathbf{C}^{\beta, 0'}$! If the ranges for forming the two densities are different (as it would be when $N_{\alpha} \neq N_{\beta}$), then having an initial pair of identical MO coefficient matrices is not an issue. For your Huckel guess (called the Generalized Wolfsberg-Helmholtz or GWH guess in Q-Chem, scf_guess = gwh), the initial MO coefficient matrices would still be identical, because you modify a quantity they have in common (the core Hamiltonian).

So, we must modify the MO coefficients (coefficient matrices) themselves. There are three options that I know of:

• mixing
• rotation
• switching

These are variants of the same thing; rotation is the most general, since you can define switching or mixing of two orbitals via rotations.

## Mixing

The most straightforward approach of the three that requires no prior knowledge of your system is mixing. The usual choice is to "mix" some amount of the HOMO with the LUMO of the same spin:

\begin{align} C_{\mu,\text{HOMO}}^{\text{new}} &= \frac{1}{\sqrt{1 + k^{2}}} \left( C_{\mu,\text{HOMO}}^{\text{old}} + k C_{\mu,\text{LUMO}}^{\text{old}} \right) \\ C_{\mu,\text{LUMO}}^{\text{new}} &= \frac{1}{\sqrt{1 + k^{2}}} \left( -k C_{\mu,\text{HOMO}}^{\text{old}} + C_{\mu,\text{LUMO}}^{\text{old}} \right) \end{align}

where $k = [0, 1]$. If $k = 0$, nothing happens; if $k = 1$, the HOMO and LUMO become the symmetric and antisymmetric linear combinations of the two, respectively. In theory, this doesn't need to be the HOMO and LUMO; any two MO indices could be chosen, as long as they will result in different densities. This is also only for a single $\mu$; because we want to work with matrices, we would be working with whole rows of the MO coefficient matrix. For Q-Chem, this is $\alpha$ spin only (scf_guess_mix = k); the Gaussian manual doesn't say which spin it takes.

## Rotation

The most general form for operating on pairs of MOs is rotation. Recall the rotation matrix in two dimensions:

$$R(\theta) = \begin{bmatrix} \cos(\theta) & -\sin(\theta) \\ \sin(\theta) & \cos(\theta) \end{bmatrix}$$

This can be used to "rotate" two arbitrary molecular orbitals $i,j$ by an arbitrary angle:

$$\begin{bmatrix} C_{i}^{'} \\ C_{j}^{'} \end{bmatrix} = \begin{bmatrix} \cos(\theta) & -\sin(\theta) \\ \sin(\theta) & \cos(\theta) \end{bmatrix} \begin{bmatrix} C_{i} \\ C_{j} \end{bmatrix}$$

where again each $C$ is a row vector. This can be expended into a set of 2 equations:

\begin{align} C_{i}^{'} &= C_{i} \cos(\theta) - C_{j} \sin(\theta) \\ C_{j}^{'} &= C_{i} \sin(\theta) + C_{j} \cos(\theta) \end{align}

There are a few interesting cases.

### 180 degree rotation: change phase

$$R(180^{\circ}) = \begin{bmatrix} -1 & 0 \\ 0 & -1 \end{bmatrix}$$

\begin{align} C_{i}^{'} &= -C_{i} \\ C_{j}^{'} &= -C_{j} \end{align}

This is kind of pointless, since the energy should be invariant to a phase swap.

### 90 degree rotation: interchange

$$R(90^{\circ}) = \begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix}$$

\begin{align} C_{i}^{'} &= - C_{j} \\ C_{j}^{'} &= C_{i} \end{align}

This is the usual case when "rotate orbitals" is seen. This is equivalent to switching orbitals. I don't know about other programs, but if you do this in ORCA, the orbital energies will be swapped as well.

### 45 degree rotation: symmetric and antisymmetric linear combinations

$$R(45^{\circ}) = \begin{bmatrix} 1/\sqrt{2} & -1/\sqrt{2} \\ 1/\sqrt{2} & 1/\sqrt{2} \end{bmatrix}$$

\begin{align} C_{i}^{'} &= \frac{1}{\sqrt{2}} \left( C_{i} - C_{j} \right) \\ C_{j}^{'} &= \frac{1}{\sqrt{2}} \left( C_{i} + C_{j} \right) \end{align}

This last example is identical to the $k = 1$ (100%) case from the mixing example shown earlier.

## Switching

We could also alter occupation numbers, such as by placing the $\alpha$ (or $\beta$) HOMO electron into the same-spin LUMO, or some other occupied orbital, but this isn't an automated procedure like orbital rotation, because you need to know your occupied/virtual orbital indices prior to starting the calculation. This is equivalent to switching orbitals.

A short summary of program capabilities (to the best of my knowledge):

• Q-Chem: mix, switch
• Gaussian: mix, switch
• GAMESS: mix, switch
• DALTON: switch
• MOLPRO: rotate
• ORCA: rotate

An important note about mixing: it doesn't matter that choosing the HOMO and LUMO to permute may be incorrect. Programs that have baked-in support for symmetry (DALTON, MOLPRO, unsure about GAMESS) will only mix/switch/rotate orbitals within a symmetry unless symmetry is disabled, but the point is only to break symmetry in the initial guess so that the local RHF minimum is escaped. If you know beforehand that your system has particular symmetries that need to be preserved, then you would probably be doing something more complex than just mixing the HOMO and LUMO anyway.