Are the canonical orbitals of Hartree-Fock also the natural orbitals?

Question

My question stems from the comments on my answer to After a unitary transformation, is Koopmans' theorem still valid?. There was some confusion relating to differing terminology referring to different, mostly unrelated density matrices. As of now, I'm uncertain as to what density matrix one diagonalizes to obtain the natural orbitals. Is it the density matrix that appears in the HF calculation or some different density matrix?

Answer 1

There is a bit of a terminology problem in the field that makes things very confusing and I will try to clarify some of this here. Part of the problem arises from the fact that sometimes only one kind of density matrix is relevant and so people drop the adjectives when talking about specific approximations.

The Density Matrix

We're talking electronic structure here (density matrices are used in other areas of quantum and physical chemistry. Honestly, it's the same density matrix but people use different language and representations and normalizations to discuss them and everyone gets confused when non-EST and EST people talk to each other about density matrices) so I will restrict myself to that topic.

The ($n$-particle) density matrix is given by: $$\rho=|\Psi(x_1,x_2,\cdots x_n)\rangle\langle\Psi(x_1,x_2,\cdots x_n)|$$

where $x_i$ is the coordinates (space and spin) of the $i$th electron. That is, it gives the probability that the system is in a state with a given set of coordinates (space and spin) for all electrons in the system. The density matrix has many great properties, including being Hermitian and having the property that $\text{tr}{\rho A}=\langle A\rangle$ for any operator $A$. In other words, the $n$-particle density specifies everything about the wavefunction of your system. If you have the full wavefunction, you have the full $n$-particle density and vice versa.

Now, most operators are not $n$-electron operators. For example, the kinetic energy operator is a 1-electron operator (only the coordinates of a single electron are needed to calculate a kinetic energy), the standard electronic Hamiltonian is a 2-electron operator (because the Coulomb potential is only dependent on the positions of two electrons at a time). Hence, we don't need the full $n$-particle density matrix to get energies and properties.

The 1(or 2 or 3...)-particle (reduced) density matrix

Since we don't need all of those other coordinates we can integrate them out. Indeed, the trace operation that we described above would do just that. Say we want the expectation value of a 1-particle operator (like kinetic energy). We can define a 1-particle reduced density matrix (sometimes the word "particle" or "reduced" is omitted), often abbreviated as 1PDM, 1RDM, $P^{(1)}$, or $\gamma^{(1)}$ (also others).

$$P(x_1)=N\int \rho(x_1,x_2,\cdots x_n)dx_2\cdots dx_n$$

where $N$ is some normalization (commonly the number of electrons in the system ($n$) because, typically, people want the trace of the 1PDM to be the number of electrons in the system).

This is the matrix which is diagonalized to form the natural orbitals.

Relatedly, we can also define a 2-particle reduced density matrix (aka 2PDM, 2RDM, $\Gamma$, $P^{(2)}$, others) as

$$\Gamma(x_1,x_2)=N\int \rho(x_1,x_2,\cdots x_n)dx_3\cdots dx_n$$

where $N$ is a different normalization (usually, $n(n-1)$). Most operators that we care about (like the energy) are 2-particle operators so the problem of electronic structure is determining the 1 and 2PDMs.

The Hartree-Fock density

There are many ways to think about the Hartree-Fock approximation, but one way that is related to density matrices is that the HF approximation is that the 2-particle, 3-particle, etc. density matrices can be expressed as direct products of the 1PDM. Intuitively, since HF ignores electron correlation, it makes sense that the probability of finding electron 1 at position $x_1$ and electron 2 at position $x_2$ is related to just the product of the those individual probabilities (with appropriate fermionic antisymmetrization). In this way, the HF 1PDM specifies all density matrices in the HF approximation (which is why in HF you can compute $J$ and $K$ which are two-electron operators by contracting twice with the HF density).

$$P(x,x') = \sum_i^\text{occ} \phi_i(x)\phi_i^*(x')$$

where $\phi_i$ are the spin-orbitals. Note that this means for HF, the canonical orbitals are the natural orbitals. For a correlated method though, the 1PDM will generally be different. For example, CI methods like CAS express the 1PDM as the weighted sum of several HF-like 1PDMs, with coefficients determined by the CAS-CI optimization. Also, note that for Kohn-Sham DFT, the noninteracting reference system also looks like a system without electron correlation and so the 1PDM specifies terms that look like higher order density matrices. However, one does not obtain the 2PDM with KS DFT methods. One obtains the energy from the 1-particle density and an exchange-correlation functional although $i$-particle operators value can be extracted from this non-interacting $i$PDM.

---

That was a really long answer, but I sense there is a lot of confusion about these different density matrices, so I figured I would take the opportunity to give a very clear answer. In short, the answer to your question is that the HF (1-particle) density specifies all of the HF $i$-particle densities. Generally though, the $i$-particle densities are determined somewhat separately, as each are impacted by electron correlation differently. The orbitals that diagonalize the 1-particle density give the natural orbitals.