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}
$$
leads to
$$\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}
$$
leads to
$$\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}
$$
leads to
$$\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.
