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.