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.