Printing of molecular orbitals


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 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

Thanks that will get us going in the right direction. Will this only work effectively without symmetry?

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:

scf_wfn.Ca_subset("AO", "ACTIVE").print_out()

Ok, that might be an ok workaround. thanks!

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 =
    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

# 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 ''
                    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:

Wow, thanks Andy we’ll give it a shot.

best, -Kirk

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?


Never mind, I just found out that it can be easily printed by using
print_MOs true

Good deal, I was in the process of checking with my student Rulin to see if he had done anything differently.