Chemistry - Column order, signs of MO coefficients, and density matrix construction in Hartree-Fock
Solution 1:
In the future, this should really be two separate questions, but I will try and tackle both.
MO coefficient column ordering and signs
The electronic energy is invariant to orbital swapping within subspaces. If you interchange one occupied orbital with another, the total energy of a molecular won't change. The same goes for virtual-virtual. Now, you know that the MO coefficients are by convention shaped as $C_{\mu i}$, where the rows are the AOs and the columns are the MOs. What this means is that for $\ce{H2O}$ in a STO-3G basis, where there are 7 AOs forming 7 MOs where 5 are doubly-occupied and 2 unoccupied, switching the order of the first 5 columns should not change the electronic energy. The same goes for the final 2 columns.
As you say, the $\epsilon_i$ follow the same ordering as the columns in $C_{\mu i}$; this is because they are the matching eigenvalues and eigenvectors from diagonalizing the Fock operator. The eigenvectors are not going to be affected by a sign change, which doesn't change how the space is spanned.
Regarding why the ordering is not what you expect, there is no guarantee made by NumPy or the underlying LAPACK routine DGEEV that the results will be sorted, just that the vector/value pairs will be in the same order. If you want to sort the results, you can use argsort:
idx = energies.argsort()
energies = energies[idx]
eigvecs = eigvecs[:, idx]
Make sure the ordering of the MO energies is maintained along with the ordering of the MO coefficient matrix columns.
Constructing the density matrix
This is a very irritating difference in convention among many sources. See later, for calculating the electronic part of the dipole moment:
$$ \left<\vec{\mu}\right> = 2 \sum_{\mu\nu} D_{\mu\nu} \left< \phi_{\mu} | \hat{\mu} | \phi_{\nu} \right>, $$ ...
- The factor 2 appearing above arises because the definition of the density used in this project differs from that used in Szabo & Ostlund.
This can also have an effect on whether or not there is a 1/2 in front of the $2J - K$ portion of the Fock operator.
Solution 2:
For your first question: since (if $A x = \lambda x$)
$$A (- x) = - A x = - (\lambda x) = \lambda (- x)\text{,}$$
an eigenvalue with sign flipped is still an eigenvalue. (In fact, in linear algebra, any multiplicative constant would be fine, but since we are interested in normalised states this reduces to $\pm 1$).
So no problem not matching signs.
For your second question: you're mixing stuff. Your first equation sums
$$\sum_m^{occ}$$
which means 10 terms in your case. Since every $\alpha$ spin-orbital contributes the same as another spatially-matching $\beta$ spin-orbital, you end up summing 10 terms in which each two are numerically the same.
Your second equation sums
$$2 \sum_m^{N/2}$$
which simply avoids this double counting: you sum 5 terms, each being twice the values in the last sum.