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

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?