is there any convenient way to print the MOs from an SCF calculation that includes the basis function labels? For example, I’d like an easy way to see if a given MO is predominately d_z2 or d_x2-y2 or noting f orbital components in an MO that also includes p-type contributions.

Unfortunately, I do not believe there is anything like that baked in. If you are up for a bit of coding you may be able to extract it from the following:

molecule mol {
He
He 1 3.0
symmetry c1
}
scf_e, scf_wfn = energy('SCF/sto-3g', return_wfn=True)
scf_wfn.Ca().print_out() # Print orbital information
scf_wfn.basisset().print_detail_out() # Print basis set detailed information

Im not sure we ever generate the linear combinations in a human readable form. If you just want to compute in symmetry and then print the orbitals without symmetry, you can do the following:

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)

I reopen this question.
Is there now an easier way to print MOs with the contribution of AOs instead of using the very user-unfriendly scf_wfn.Ca().print_out() ? @Kirk: Can you share how you get it?