It’s not exactly what you want, but I cooked together some code that you can dump into an input file to find the composition of any symmetrized atomic orbitals. This may help with parsing the SO index of the wfn.Ca() array:
E,wfn = energy('hf/cc-pvtz', return_wfn=True)
#
# Cartesian AO labels
#
ao_labels = []
for s in range(wfn.basisset().nshell()):
shell = wfn.basisset().shell(s)
center = str(shell.ncenter+1)
am = shell.am
amchar = shell.amchar
basename = '{'+center+'}'+amchar
for j in range(0,am+1):
lx = am - j
for lz in range(0, j + 1):
ly = j - lz
ao_labels.append(basename+'x'*lx+'y'*ly+'z'*lz)
#
# SO labels
#
so_labels = []
sobasis = wfn.sobasisset()
petitelist = sobasis.petite_list()
U = petitelist.sotoao().nph
for irrep in range(wfn.nirrep()):
mat = U[irrep]
dims = mat.shape
for so in range(dims[0]):
soterms = ''
for ao in range(dims[1]):
aofunc = ao_labels[ao]
coef = mat[so,ao]
if(np.abs(coef) > 1E-10):
if not soterms:
join = '-' if coef < 0 else ''
else:
join = '- ' if coef < 0 else '+ '
soterms += " %s%.3f(%s)" % (join, np.abs(coef), aofunc)
so_labels.append(str(so)+'['+str(irrep)+']' + soterms)
#
# Print the results as
# number within irrep,[irrep], components
# where compnents are listed as
# coefficient({center} Cartesian component)
# N.B. ALL indices are 0 based, except atom center numbers
#
for sol in so_labels:
print(sol)