Column order, signs of MO coefficients, and density matrix construction in Hartree-Fock
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&lt;\vec{\mu}\right&gt; = 2 \sum_{\mu\nu} D_{\mu\nu} \left&lt; \phi_{\mu} | \hat{\mu} | \phi_{\nu} \right&gt;,
$$
  ...
  
  
  The factor 2 appearing above arises because the definition of the density used in this project differs from that used in Szabo &amp; 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.