Orthogonalize first: X = S^(−1/2)
Learn
There's a classic maneuver for taming FC = SCε: manufacture a matrix X that makes the basis orthonormal, and the generalized problem becomes an ordinary one. The symmetric choice is the inverse square root of the overlap, built from S's own eigendecomposition:
Diagonalize F′ with a garden-variety symmetric eigensolver, then map back with C = XC′. For H₂, with its huge overlap S₁₂ = 0.659, the actual X is:
Those large off-diagonals are X actively un-mixing two heavily overlapping functions. (In practice, scipy.linalg.eigh(F, S) does this internally — but knowing what's inside the black box is the difference between using Roothaan-Hall and understanding it.)
In practice, how do you solve the generalized eigenvalue problem FC = SCε for a non-orthogonal Gaussian basis?
Either route, one diagonalization hands you all the orbitals at once — occupied and virtual.
Diagonalizing the Roothaan-Hall equations in an N-function basis returns N orbitals. The N/2 with the lowest energies are occupied. What are the rest?