Alternative Program for Orbital Population Analysis 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.
ORCA can perform population analysis for each individual molecular orbital. This is an excerpt from the output file:

-----------------------------------
MULLIKEN ORBITAL POPULATIONS PER MO
-----------------------------------
THRESHOLD FOR PRINTING IS 0.1%
SPIN UP
                      0         1         2         3         4         5
                 -177.49372 -100.54186 -100.54186 -100.54184 -100.54179 -100.54174
                   1.00000   1.00000   1.00000   1.00000   1.00000   1.00000
                  --------  --------  --------  --------  --------  --------
  0Ti  1s             99.6       0.0       0.0       0.0       0.0       0.0
  0Ti  2s              0.4      -0.0      -0.0      -0.0      -0.0      -0.0
  1Cl  1s             -0.0       4.5       3.8       0.0      48.9       0.4
  1Cl  2s              0.0       0.0       0.0       0.0       0.3       0.0
  2Cl  1s             -0.0       4.9       4.1       0.0      50.5       0.4
  2Cl  2s              0.0       0.0       0.0       0.0       0.3       0.0
  3Cl  1s             -0.0      69.2       3.8       0.0       0.0      23.3
  3Cl  2s              0.0       0.5       0.0       0.0       0.0       0.2
  4Cl  1s             -0.0       1.4      71.5       0.0       0.0      23.4
  4Cl  2s              0.0       0.0       0.5       0.0       0.0       0.2
  5Cl  1s             -0.0       8.9       7.5      50.7       0.0      26.2
  5Cl  2s              0.0       0.1       0.1       0.3       0.0       0.2
  6Cl  1s             -0.0      10.4       8.7      48.6       0.0      25.6
  6Cl  2s              0.0       0.1       0.1       0.3       0.0       0.2


The program supports all options that you listed, and this input file header should get you started:

! B3LYP def2-SVP ECP{def2-SVP} COSMO(water) LargePrint
*xyz charge multiplicity
coordinates go here
*


You can find more detailed information in the manual.