Gaussian has Orbital-by-Orbital Population Analysis.
From the manual:
# UHF/6-311+G(d) Pop=Orbitals=3
Here is the resulting output from a calculation on FeO+ quartet:
Atomic contributions to Alpha molecular orbitals:
Alpha occ 16 OE=-0.923 is Fe1-d=1.00
Alpha occ 17 OE=-0.699 is O2-p=0.88
Alpha occ 18 OE=-0.690 is O2-p=0.68 Fe1-s=0.21
Alpha vir 19 OE=-0.253 is Fe1-s=0.70 Fe1-p=0.27
Alpha vir 20 OE=-0.188 is Fe1-p=0.71 O2-p=0.29
Alpha vir 21 OE=-0.133 is Fe1-p=1.04
Atomic contributions to Beta molecular orbitals:
Beta occ 13 OE=-0.801 is O2-p=0.79
Beta occ 14 OE=-0.783 is Fe1-d=1.00
Beta occ 15 OE=-0.758 is O2-p=0.89
Beta vir 16 OE=-0.241 is Fe1-s=0.81 Fe1-p=0.17
Beta vir 17 OE=-0.139 is Fe1-p=0.91 Fe1-d=0.14
This can be used to provide a more quantitative characterization of the nature of canonical molecular orbitals.
Due to certain restrictions on the use of Gaussian (http://www.bannedbygaussian.org/), I am unable to use Gaussian. Are there other programs that can perform such a calculation? NBO only seems to provide such a breakdown in terms of NBOs, and I can't find the feature available in any other quantum chemistry program. I realize that I can do this myself by extracting the C matrix, figuring which coefficients go with which orbital, and summing up myself. However, this would probably take me a day to program.
Note: My calculations use B3LYP, ECPs, and the COSMO solvent model. So the program would need to be able to at least produce single points with such a method.