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?


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


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



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


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.

  • $\begingroup$ If $\mathbf{B}$ is the "normal" coefficient matrix in the Roothan-Hall matrix equation, then the density matrix we use in our calculations is $\mathbf{C}\mathbf{C}^{\dagger}$Is $\mathbf{T}$ only over the occupied molecular orbitals. Is this what is referred to as the "reduced density matrix"? The $\mathbf{T}$ matrix in your answer is not what we get when we diagonalize the fock matrix, right? $\endgroup$ – Yoda Mar 6 '20 at 11:08
  • $\begingroup$ The density is just built from the occupied MO coefficients (see this tutorial). While you can solve the HF equation with C instead of T, the Fock matrix itself only depends on the occupied orbitals. $\endgroup$ – Tyberius Mar 6 '20 at 14:10
  • 1
    $\begingroup$ Related discussion on density matrices: chemistry.stackexchange.com/q/80351/41556 $\endgroup$ – Tyberius Mar 6 '20 at 14:11

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.

  • $\begingroup$ While your derivation makes sense, the choice of $X = S^{-1/2}$ is not the only possible orthogonalization, and I think that the final equation does not depend on the choice. $\endgroup$ – TAR86 Mar 6 '20 at 8:07
  • $\begingroup$ The transformation of the Fock matrix into an orthonormal AO basis is fine. But it is not clear to me why a very similar transformation of the density matrix changes the basis from MO to AO (and vice versa), and not just between AO and orthogonalized AO basis. $\endgroup$ – Yoda Mar 6 '20 at 17:04
  • $\begingroup$ Just noticed that you specifically asked about the MO basis rather any general orthogonalised basis. Hence why I started in the orthogonal AO basis. Will update to explain connection to MO basis shortly $\endgroup$ – user213305 Mar 6 '20 at 17:24
  • $\begingroup$ @user213305 Did you update your answer? $\endgroup$ – Yoda Mar 28 '20 at 11:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.