How to print natural occupation number

We need to know the method you want to print occupation numbers of.

For HF, the spinorbital occupation numbers are all 1 or 0, as a consequence of the wavefunction being a single determinant. It looks like you develop electronic structure methods. Am I misunderstanding the question, or did you just need a reminder of a fairly basic fact?

For CASSCF (or to doublecheck HF), you can access the density matrix as a numpy array by `wfn.Da_subset("MO").np`

(for a molecule with no symmetry) or as an element of the *list* `wfn.Da_subset("MO").nph`

(for a molecule with symmetry). From there, you determine the occupation numbers by diagonalization. This requires that the method you’re using set the density matrix properly. Both CASSCF and HF do.