Column order, signs of MO coefficients, and density matrix construction in Hartree-Fock

Question

I'm attempting to make a Hartree-Fock solver, using Python, following Daniel Crawford's guide.

I have been successful with the first four steps. However, I am having some trouble with this part in step 5:

$$ F^0_c C^{'}_0 = C^{'}_0 \epsilon_0 $$

I have the correct transformed Fock matrix (I can verify this as the site as step-by-step matrices that you can check against). My problem is that the MO coefficient matrix I get does not quite match the one on the site. I obtain it by multiplying the above equation by $C^{' -1}_0$ on the right:

$$ F^0_c = C^{'}_0 \epsilon_0 C^{' -1}_0 $$

Now it's in diagonalized form, which means $C^{'}_0$, the transformed MO coefficients, is the eigenvectors of $F^0_c$ and $\epsilon_0$ is a diagonal matrix of the eigenvalues. After I obtain $C^{'}_0$, I obtain the original untransformed MO coefficients $C_0$ using the next equation on the site:

$$ C_0 = S^{-1/2} C^{'}_0 $$

where $S^{-1/2}$ is the transformation matrix used to get the transformed Fock matrix from earlier. My transformation matrix, by the way, matches exactly the one on the site.

This is where my problem comes in. The $C_0$ matrix I get is this:

The one the website has as a reference is this:

If you look carefully, you will notice that mine and the reference one have the same columns; however, some are swapped around, and some have sign flips. This causes problems with the next steps, such as finding the density matrix.

I thought maybe the order should be such that $\epsilon_0$ is in order. The $\epsilon_0$ matrix I get is indeed out of order, and the same permutations and sign flips that will make it in order also makes my $C_0$ match the reference one in terms of where each column is, but NOT sign. Sadly the reference site does not have a $\epsilon_0$ matrix or $C{'}_0$ for me to compare with, only the final $C_0$, so I'm not sure if those are coming out of order also.

My first question is, will these sign flips matter? Is ordering the $C_0$ matrix in terms of getting $\epsilon_0$ the right thing to do?

I am carrying out these matrix operations with basic NumPy functions such as numpy.linalg.eig and numpy.dot. Nothing fancy there. It's worth noting that finding the transformation matrix, $S^{-1/2}$, requires diagonalizing in a similar way, but I didn't have any column switching problems there. It came out exactly as expected.

My second question relates to building the density matrix at the end of step 5. The website has this equation:

$$ D^0_{\mu\nu} = \sum_{m}^{occ} (C_0)^m_\mu (C_0)^m_\nu $$

It says sum over the occupied MOs. The molecule being used in the example calculation is water, so two H's, and one O, so there are 10 electrons, so 5 orbitals. If I sum using the correct $C_0$ matrix I get on the reference website, from m = 1 to 5, I get the correct $D^0_{\mu\nu}$ result. However my textbook (Szabo and Ostlund, Modern Quantum Chemistry) has this equation (changing some of the notation to match):

$$ D^0_{\mu\nu} = 2 \sum_{m}^{N/2} (C_0)^m_\mu (C_0)^m_\nu $$

Where N is the number of electrons. Well, 10/5 is 5, so I'm summing from 1 to 5 again, but now there's a 2 in front, so the sum will be doubled. Now, using the one on the website gets me the correct answer, but what could this 2 in the book mean?

Answer 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>, $$ ...

2. 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.

Answer 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.