# Orthonormality of AOs and MOs in PySCF

I have a basic question about how to calculate molecular-orbitals from atomic-orbitals in PySCF using the self-consistent field approach. My understanding mathematically is as follows: we start with a set of non-orthonormal atomic-orbitals $$\{\phi_i\}$$ with overlap matrix $$S_{ij} = \langle \phi_i, \phi_j \rangle$$. Application of Hartree-Fock produces a new set of orbitals $$\{\psi_i\}$$ via a unitary orbital transformation $$U$$. The relationship between the old and new orbitals is: $$\psi_i = \sum_j U_{ij} \phi_j.$$

Since $$U$$ is unitary, the overlap matrix should be invariant under this transformation, implying that:

$$\langle \psi_i, \psi_j \rangle = \sum_{kl} U_{ik} U_{jl} S_{kl} = S_{ij}.$$

However, when I check this, I get that the molecular-orbitals are actually orthonormal, as their overlap matrix is the identity (see below for the code I used), and moreover, the transformation $$U$$ I obtain is not unitary. My question is: since the self-consistent field method should produce a unitary transformation of the orbitals, what is the approach used to obtain a non-unitary transformation. A follow-up question is: is there a benefit to orthonormalizing the MOs, rather than orthnormalizing the AOs first before obtaining the MOs?

import numpy as np
from pyscf import gto, scf

mol = gto.M(
atom='H 0 0 0 ; Li 0 0 1',
basis='sto-3g'
)

hf = mol.RHF()
hf.scf()
S = hf.get_ovlp()
U = hf.mo_coeff.T # this defines the transformation from AOs to MOs

S_MOs = np.einsum('ik,jl,kl->ij', C, C, S)

# S_MOs is an identity matrix
$$$$
`

In typical programs, this happens by the process of canonical orthonormalization due to Löwdin; see section 7 of the open access review Molecules 25, 1218 (2020). PySCF does not employ this approach by default; instead, it follows an alternative but equivalent approach of solving the generalized eigenvalue problem of the Roothaan equation $${\bf FC} = {\bf SCE}$$ with the eigh function of SciPy.
In either case, you will find out that the produced orbital coefficients satisfy $${\bf C}^\text{T} {\bf SC} = {\bf 1}$$, that is, the orbital coefficients are orthonormal in the overlap metric. Note that $${\bf C}^\text{T} {\bf C} \neq {\bf 1}$$ except in the case of an orthonormal basis, $${\bf S} = {\bf 1}$$,