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

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?

• Welcome to Chemistry! Take the tour to get familiar with this site. Rather than posting screenshots of your terminal, what you can do is take the matrix blocks and add 4 spaces to the front of each line to get monospace text. – pentavalentcarbon Jun 12 '17 at 21:27
• Why do all these tutorials always provide the integrals? The most interesting (and expensive) part in programming a QC Program is the calculation of the integrals.... – Fl.pf. Jun 13 '17 at 6:58
• I agree in an academic sense, but practically, it's useful to have the values to make sure the rest of your code is working. When I'm done with the other parts, I'm going to make some actual integral-evaluating code. folk.uio.no/helgaker/talks/SostrupIntegrals_10.pdf This is something I have found helpful in evaluating the simpler one-electron integrals. For the two-electron integrals I may just link up to psi4 or a similar program, as that's an entire field in itself... – iammax Jun 13 '17 at 15:20
• I use both Psi4 and pyscf in a project of mine to generate 1- and 2-electron integrals. The problem is the other side of the same coin: you can work on methods without even knowing about integral evaluation and vice-versa. Implementing an integral engine is tedious for the naive case and obtuse/very difficult for the optimized case, and slow integrals are really slow. On Github, I have an Obara-Saika program and Joshua Goings has a McMurchie-Davidson program you can look at in Python. – pentavalentcarbon Jun 13 '17 at 17:26
• I'm likely to go with psi4 because my PI is (or rather, was) one of the developers so he has a lot of experience with it, but yes, actually writing integral-evaluating functions yourself will probably not be very good, although I am going to do it anyway, just for the learning experience of how the integrals are evaluated, which may give me some insight into how they are better optimized, even if I don't implement those optimizations myself. – iammax Jun 13 '17 at 20:03

## 2 Answers

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

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

• You're right, the factor of 2 was because of the notation difference. In the online guide, it was 2J -K, whereas in the book, it was J - K/2. I didn't even think to look at that before, but that accounts for the difference. – iammax Jun 12 '17 at 22:22

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.