Hartree-Fock Method: Projection of density matrix onto a larger basis set

Question

I am implementing a Hartree-Fock SCF program with an educational purpose. Although the code works fine using a Core Hamiltonian guess for the SCF iterative process, I want to implement a more accurate initial guess, so I have programmed a simple Extended Hückel guess using a minimal basis set (MINI) and the parameters used by the GAMESS-US software as the diagonal elements of the matrix. This guess works fine for Hartree-Fock SCF calculations using minimal basis sets, but I want to generalize it for any type of basis set (my code support some split-valence basis sets, such as 3-21G). The manuals of some popular quantum chemistry packages state that the density matrix guess obtained using the Extended Hückel Method can be projected onto the desired basis set, but not explaining explicitly how it's done. After some search, I've found that the pyscf package implements a function called project_dm_nr2nr() which does exactly what I'm looking for. The guide of the program (https://pyscf.org/pyscf_api_docs/pyscf.scf.html) states that the projection is done using the equation: $$ D_{AO2}=PD_{AO1}P^T $$ with $$ P=S^{-1}_{AO2}\langle AO2|AO1\rangle $$ where $D_{AO2}$ is the electronic density matrix expressed in the new atomic orbital basis,$D_{AO1}$ is the density matrix expressed in the minimal basis set (in which the extended Hückel calculation has been done), $S^{-1}_{AO2}$ is the inverse of the overlap matrix of the new basis set and $\langle AO2|AO1\rangle$ is the overlap matrix between the two basis sets.

Although I've been searching information about the topic, I haven't found any information about where this formula(in textbooks about quantum chemistry nor the Internet) comes from and I wouldn't like to implement something that I don't understand in my program.I suspect that it is a basic question and I lack the mathematical knowledge for deriving the formula.Could you help me deriving this equation?

Answer 1

The spin-$\sigma$ density matrix is given by $$ \mathbf{P}^\sigma = \mathbf{C}_\text{occ}^\sigma (\mathbf{C}_\text{occ}^\sigma)^\text{T} $$ so if we know how the orbital coefficients $\mathbf{C}^\sigma$ transform, then we know how $\mathbf{P}^\sigma$ transforms.

Now, the orbitals are given by $$ |\psi_{i \sigma} \rangle = \sum_\alpha C_{\alpha i}^\sigma |\chi_\alpha\rangle $$ Inserting the resolution of the identity $$ \mathbf{1} = \sum_{AB} |A\rangle \langle A | B \rangle^{-1} \langle B| $$ from the left of the right-hand side we get $$ |\psi_{i \sigma} \rangle = \sum_\alpha C^\sigma_{\alpha i} |A\rangle \langle A | B \rangle^{-1} \langle B|\chi_\alpha\rangle \equiv \sum_A C^\sigma_{Ai} |A\rangle $$ from which we can read the transformed coefficients $$ C^\sigma_{Ai} = \sum_\alpha \langle A | B \rangle^{-1} \langle B|\chi_\alpha\rangle C^\sigma_{\alpha i}. $$

If I call the basis $|\chi_\alpha \rangle$ basis 1, and the basis $|A\rangle$ basis 2, this equation can be rewritten in compact matrix notation as $$ \mathbf{C}^\sigma_{2} = \mathbf{S}_{22}^{-1} \mathbf{S}_{21} \mathbf{C}^\sigma_{1}. $$

The equation for projecting the density matrix follows trivially: $$ \mathbf{P}^\sigma_{2} = \mathbf{S}_{22}^{-1} \mathbf{S}_{21} \mathbf{C}^\sigma_{1\text{, occ}} \left(\mathbf{S}_{22}^{-1} \mathbf{S}_{21} \mathbf{C}^\sigma_{1\text{, occ}}\right)^\text{T} = \mathbf{S}_{22}^{-1} \mathbf{S}_{21} \mathbf{C}^\sigma_{1\text{, occ}} (\mathbf{C}^\sigma_{1\text{, occ}})^\text{T} \mathbf{S}_{12}^\text{T} = \mathbf{S}_{22}^{-1} \mathbf{S}_{21} \mathbf{P}^\sigma_{11} \mathbf{S}_{12} \mathbf{S}_{22}^{-1} $$

Note that these equations don't make any assumptions of the size of the basis sets, or the level of theory used.