Computing the commutator of the Fock and density matrices in AO basis (Hartree-Fock theory)

I am making a simple restricted HF code using the Python interface of Psi4. I am now evaluating convergence by tracking the change in the sum of orbital energies, but I want to do this in a better way. It is common to use the fact that at self-consistency, the Fock and density matrices commute


$$
[\mathbf{F}, \mathbf{D}] = \mathbf{FD} - \mathbf{DF} = \mathbf{0}
$$


However, the above expression is only valid in MO basis, while F and D in my code are computed in AO basis. So I need to derive an equivalent expression in AO basis. I am quite sure the correct expression is


$$
[\mathbf{F}, \mathbf{D}]^{\text{AO}} = \mathbf{FDS} - \mathbf{SDF}
$$


as this is equal to zero to within $1\times10^{-14}$. But how to derive this?


Derivation
----------


An arbitrary molecular orbital $\phi\_i$ is expanded in atomic orbital basis functions


$$
\phi\_i = \sum\_\alpha C\_{\alpha i} \chi \_{\alpha}
$$


Acting the commutator on $\phi\_i$ and expanding it to AO basis yields


$$
[\mathbf{F}, \mathbf{D}] = \mathbf{FD} \sum\_\alpha C\_{\alpha i} \chi\_{\alpha} - \mathbf{DF} \sum\_\alpha C\_{\alpha i} \chi\_{\alpha}
$$


Since we know the solution contains the overlap matrix $\mathbf{S}$, lets look at the definition


$$
\mathbf{S}\_{ij} = \langle \chi\_i(\mathbf{r}) \vert \chi\_j(\mathbf{r}) \rangle = \int d\mathbf{r} \chi^\*(\mathbf{r})\chi(\mathbf{r})
$$


Since this must be part of our expression, it seems to me a good approach is to multiply from the left by $\sum\_\beta C\_{\beta i}^\* \chi\_\beta^\*$ (dropping the $\mathbf{r}$ dependence from now on) and integrating over $\mathbf{r}$


$$
[\mathbf{F}, \mathbf{D}] = \int \sum\_\beta C\_{\beta i}^\* \chi\_\beta^\* \mathbf{FD} \sum\_\alpha C\_{\alpha i} \chi\_{\alpha} - \int \sum\_\beta C\_{\beta i}^\* \chi\_\beta^\* \mathbf{DF} \sum\_\alpha C\_{\alpha i} \chi\_{\alpha}
$$


from in braket notation becomes


$$
[\mathbf{F}, \mathbf{D}] = \langle \sum\_\beta C\_{\beta i} \chi\_\beta \vert \mathbf{FD} \vert \sum\_\alpha C\_{\alpha i} \chi\_{\alpha} \rangle - \langle \sum\_\beta C\_{\beta i} \chi\_\beta \vert \mathbf{DF} \vert \sum\_\alpha C\_{\alpha i} \chi\_{\alpha} \rangle
$$


At this point I am not sure what to do - or if am I even on the right track. I can see that we have the "pieces" that make up the overlap matrix, but I don't know how to put them together. Further, due to the orthonormality of the MOs, then I can imagine that the summation terms only survive when $\alpha = \beta$. But I'm not sure how to derive this properly.




Start from the AO Hartree-Fock equation and its adjoint 
$$\mathbf{F}^{AO}\mathbf{T}=\mathbf{S}\mathbf{T}\epsilon \text{ and } \mathbf{T^\dagger}\mathbf{F}^{AO}=\epsilon\mathbf{T^\dagger}\mathbf{S}$$


where $\mathbf{T}$ is an $N\times n$ matrix that is essentially the occupied block of $\mathbf{C}$ ($n$ is occupied, $N$ is total orbitals). We use this $\mathbf{T}$ matrix because it has the convenient property that 


$$\mathbf{T}\mathbf{T}^\dagger=\mathbf{D}^{AO}$$


Now, we can multiply the HF equation by $\mathbf{T^\dagger}\mathbf{S}$ on the right and multiply its adjoint by $\mathbf{S}\mathbf{T}$ on the left, which gives


$$\mathbf{F}^{AO}\mathbf{T}\mathbf{T^\dagger}\mathbf{S}=\mathbf{S}\mathbf{T}\epsilon\mathbf{T^\dagger}\mathbf{S}$$


$$\mathbf{S}\mathbf{T}\mathbf{T^\dagger}\mathbf{F}^{AO}=\mathbf{S}\mathbf{T}\epsilon\mathbf{T^\dagger}\mathbf{S}$$


Subtracting the first equation from the second and yields the desired commutator relationship.


$$\mathbf{F}^{AO}\mathbf{D}^{AO}\mathbf{S}-\mathbf{S}\mathbf{D}^{AO}\mathbf{F}^{AO}=\mathbf{0}$$


The notation I use here is based on a similar derivation given in Chapter 6 of McWeeny's Methods of Molecular Quantum Mechanics, 2nd edition.


 6

Starting from statement that the Fock matrix and the density matrix commute in an orthonormal basis.
$$
[\mathbf{F}, \mathbf{D}] = \mathbf{FD} - \mathbf{DF} = \mathbf{0}
$$
The orthonormal basis matrices can be substituted for their equivalents in an atomic orbital basis
\begin{align}
\mathbf{F} = {} & \mathbf{X}^\dagger \mathbf{F}^{AO} \mathbf{X} \\
\implies \mathbf{F}^{AO} = {} & \left[\mathbf{X}^\dagger\right]^{-1} \mathbf{F} \mathbf{X}^{-1}\\
\mathbf{D}^{AO} = {} & \mathbf{X} \mathbf{D} \mathbf{X}^\dagger \\
\implies \mathbf{D} = {} & \mathbf{X}^{-1} \mathbf{D}^{AO} \left[\mathbf{X}^\dagger\right]^{-1} 
\end{align}


where $\mathbf{X}$ is an orthogonalisation matrix. :
\begin{align}
[\mathbf{F}, \mathbf{D}] = {} & \mathbf{FD} - \mathbf{DF} = \mathbf{0} \\
= {} & \mathbf{X}^\dagger \mathbf{F}^{AO} \mathbf{X}
\mathbf{X}^{-1} \mathbf{D}^{AO} \left[\mathbf{X}^\dagger\right]^{-1}
- \mathbf{X}^{-1} \mathbf{D}^{AO} \left[\mathbf{X}^\dagger\right]^{-1} 
\mathbf{X}^\dagger \mathbf{F}^{AO} \mathbf{X} \\
= {} & \mathbf{X}^\dagger \mathbf{F}^{AO} \mathbf{D}^{AO} \left[\mathbf{X}^\dagger\right]^{-1}
- \mathbf{X}^{-1} \mathbf{D}^{AO} \mathbf{F}^{AO} \mathbf{X}
\end{align}


When $\mathbf{X} = \mathbf{S}^{-\frac{1}{2}}$, pre and postmultiplying by $\mathbf{X}^{-1} = \mathbf{S}^{\frac{1}{2}}$:


\begin{align}
\mathbf{X}^{-1} \mathbf{0} \mathbf{X}^{-1} = {} &
\mathbf{X}^{-1} \mathbf{X}^\dagger \mathbf{F}^{AO} \mathbf{D}^{AO} \left[\mathbf{X}^\dagger\right]^{-1} \mathbf{X}^{-1}
- \mathbf{X}^{-1}
\mathbf{X}^{-1} \mathbf{D}^{AO} \mathbf{F}^{AO} \mathbf{X}
\mathbf{X}^{-1} \\
\mathbf{0} = {} & \mathbf{F}^{AO} \mathbf{D}^{AO} \mathbf{S} - \mathbf{S} \mathbf{D}^{AO} \mathbf{F}^{AO}
\end{align}


as you have suggested.


 5