How to calculate the Fock matrix in the molecular orbital basis PySCF?

Question

I am interested in calculating the Fock matrix in the molecular orbital basis with PySCF, though I am not clear on the methodology behind this task.

---

In my attempt, I use the following script (for the example H$_{2}$ molecule):

from pyscf import gto, scf
geometry = '''
    H   0.000   0.000   0.000
    H   0.000   0.000   0.740
'''
mol = gto.Mole()
mol.atom = geometry
mol.basis = '3-21g'
mol.build()
​
mf = scf.RHF(mol)
mf.scf()
​
Fao = mf.get_fock()
Fmo = mf.mo_coeff.T @ Fao @ mf.mo_coeff
​
print('F_mo')
print(Fmo)

In this method, I first calculate the molecular mean-field. I then do matrix multiplication with the molecular coefficient transpose matrix (mf.mo_coeff.T), the Fock matrix in the atomic basis (Fao) and the molecular orbital coefficients (mf.mo_coeff).

---

The resulting off-diagonal matrix elements are essentially zero for the H$_{2}$ molecule and other larger systems taken to 10 decimal places (CH$_{4}$, NH$_{3}$, H$_{2}$O). This has confused me: I have seen other Fock matrices in the molecular orbital basis with off-diagonal elements present.

I am therefore looking for confirmation of my method, and if there is a better way of doing this task?

Answer 1

For canonical MOs, it is expected that the off-diagonal elements of the described product with the Fock matrix are zero (within reasonable accuracy) because they are the result of the diagonalization of the Fock matrix.

However, canonical MOs are not the only possible choice. One can "rotate" the MOs to minimize certain metrics in order to obtain, for instance, localized MOs (LMO). Such rotations will keep the total energy and the electron density intact, but will disallow the interpretation of the diagonal values of the Fock matrix (represented in the LMO basis) as orbital energies.