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.