Hartree Fock Algorithm too high values for Molecules

Question

I am working on a high school graduation paper on Hartree Fock and have gotten pretty far already. With help from this community, I managed to generate all the necessary integrals and implement the algorithm in python. I am currently checking whether the algorithm is actually doing what it's supposed to do. I compare my values with the ones provided by pyscf. See this example.

For noble gases I get correct values. If I try my implementation for molecules however I get too high values. I've tried $\ce{H2}$, $\ce{F2}$, $\ce{O2}$, and $\ce{HF}$. My algorithm works totally fine for atoms so I feel like there can't be something wrong with the implementation itself. I've also tried generating completely random initial guesses. I get the exact same results as well. I therefore feel like the algorithm is stopping on a local minimum or saddle point. Does anybody have an idea how to fix it or whether HF just returns such bad results?

This is my code if someone wants to check it out and give feedback. It only runs on MacOS and Linux distributions since pyscf can't run on Windows. Help is greatly appreciated since the due date for the paper is getting closer and closer ;)

TLDR: Restricted Hartree Fock SCF Algorithm returns too high orbital energies for molecules, but not for Atoms. How can I fix that?

Update: Here is a list of what I've tried:

1. I first checked my code by using a zero matrix as the initial guess. This I compared to the PySCF HF kernels result. The energies are too high.
2. I then checked whether the algorithm is doing anything anyway. I therefore used a completely random initial guess with values ranging between -100 and 100. The algorithm returned the same energies as with the zero matrix as the initial guess.
3. I then checked whether all the matrices are correct. I used PySCF get_hcore(), get_ovlp(), get_veff(), get_fock() for a given density matrix. All my matrices are the same as the ones generated by PySCF.

Here a quick explanation of my implementation of the algorithm:

1. Initial guess density Matrix $P$ is a zero matrix
2. Generate transformation matrix $X$ as $S^{-1/2}$.
3. Generate Fock matrix with given $P$.
4. Transform Fock matrix as $F^{'} = X^{\dagger}FX$
5. Diagonalize $F^{'}$
6. Transform coefficients $C^{'}$ back as $C = X C^{'}$
7. Use $C$ to generate new $P$ and repeat till energies meet convergence criteria. I assume that if the sum of orbital energies minus the sum of the last orbital energies is an absolute value of less than 0.0001 it converges.

There is no implementation of DIIS, level shifting or damping. I hope this makes it clearer.

Answer 1

Since you've done such an excellent job so far I thought I'd just give you some hints in the chat, but after banging your head on it for a week I think you deserve a proper answer. And with the new version of your code it's even easier to fix.

The Hartree-Fock orbitals are defined as: $$ \psi_i = \sum_j^N c_{ij}\phi_j $$ where $\phi_j$ are the basis functions, $N$ is the number of basis functions, and $c_{ij}$ are your expansion coefficients. From the orbitals, you construct the density matrix: $$ K_{ij} = \sum_k^N f_k c_{ik} c_{jk}^{*}. $$ Here the values $f_k$ have been introduced. $f_k$ is the occupation number of the orbital $\psi_k$.

Inside your code, you construct the density matrix with the following (abbreviated) routine:

for i in range(i_range):
  for k in range(k_range):
    row_i = old_matrix[i]
    row_k = old_matrix[k]
    product_array = np.multiply(row_i, row_k)
    matrix[i, k] = 2 * np.sum(product_array)

Having just reviewed the equations, you can see the missing ingredient: $f_k$. In fact, the multiplication you do by two essentially is assigning an occupation of $2$ to every orbital.

So why did it work for atoms and not molecules? It's just an artifact of the atoms you choose. Consider the Helium atom: it has one basis function and two electrons which doubly occupy the $1s$ subshell ($1s^2$). Neon was your second test, which has $10$ electrons and has configuration $1s^22s^22p^6$, which means the $1s$ is doubly occupied, as is the $2s$, as well as the $2p_x$, $2p_y$, and $2p_z$. So multiplying by $2$ worked.

Now consider the test I showed you: Carbon. It has $5$ basis functions, but just $6$ electrons ($1s^22s^22p^2$). And so the occupation numbers are $[2, 2, 2, 0, 0]$. This is why Carbon did not work, multiplication by $2$ was not sufficient. You can thus fix your density matrix routine by adding an additional argument which are the occupation numbers, and multiplying by them.