How to orthonormalize a set of Molecular orbitals?

Question

I am using PySCF, and checked that the molecular orbitals(MO) from the HF calculation are orthonormal. If S is the Overlap matrix, and V is the matrix of the MO coefficients, It can be seen from the equation, $ V^T * S * V = I $ when I is the Identity matrix. The S can be calculated from PySCF as,

S = mol.intor('int1e_ovlp')

Since $S$ is created from orthonormal orbitals, it has $1.00$ as diagonal elements. Now, the question is, I have a matrix of MO coefficients, which are not orthonormal. I was trying to orthonormalize them by using Lowdin symmetric orthogonalization. Which is actually, $ W' = S^{-1/2} W $ where, $W'$ is orthonormalized orbitals and S is calculated as before. However, the resulting orbitals are not orthonormal. Possible reasons:

1. The Lowdin method only does orthogonalization, so first, I need to normalize them
2. Use a different method, like the Gram-Schmidt method
3. Motivated from the [answer here][1], create $S$ from non-orthogonal orbitals, find $S^{-1/2}$ and do the transformation from non-orthonormal to orthonormal orbitals using $ W' = S^{-1/2} W $ . However, I am confused as to why Lowdin orthogonalization does not give the answers.

Are there other methods to do such orthonormalization? [1]: https://chemistry.stackexchange.com/questions/102144/is-the-lowdin-orthogonalization-used-in-diagonalizing-the-atomic-orbitals-really

Answer 1

Since S is created from orthonormal orbitals, it has 1.00 as diagonal elements.

If ${\bf S}$ is created from orthonormal orbitals, the overlap matrix is simply the unit matrix, ${\bf S}={\bf 1}$. When the basis functions are merely normalized, ${\bf S}$ has a unit diagonal, but typically also has non-zero off-diagonal elements.

Now, the question is, I have a matrix of MO coefficients, which are not orthonormal. I was trying to orthonormalize them by using Lowdin symmetric orthogonalization. Which is actually, $ {\bf W}′={\bf S}^{-1/2} {\bf W}$ where, ${\bf W}′$ is orthonormalized orbitals and ${\bf S}$ is calculated as before. However, the resulting orbitals are not orthonormal.

This is not surprising, since you have not orthonormalized the MO coefficients.

Let's follow convention and take ${\bf C}$ as the matrix of orbital coefficients. The MO overlap is then given by ${\bf S}^\text{MO} = {\bf C}^T {\bf SC}$. What you want to do is to make this matrix the unit matrix. This is achieved by transforming the orbital coefficients as ${\bf C} \to {\bf C}({\bf S}^{\rm MO})^{-1/2}$.