memory 1024 mb
import numpy as np
import psi4
import math

molecule mol {
0 1
O 0.00000000 0.00000000 0.00000000
H -0.75301500 0.00000000 0.56774300
H 0.75301500 0.00000000 0.56774300

}

set scf_type df
set basis 6-31G(d_p)
set reference rhf
set MOLDEN_WRITE true
set WRITER_FILE_LABEL output
set print_mos true



e, wfn = energy('scf', return_wfn=True)

newline = """
"""
print_out("Molecular Orbitals with AO's Listed"+newline)
# Use Andy's handy little AO name helper function
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 = '{'+wfn.molecule().symbol(shell.ncenter)+center+'} '+amchar
    if amchar == 'd' and wfn.basisset().has_puream() == True:
        for i in range(1,6): ao_labels.append(basename+' '+str(i))
        continue
    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)

evals = list(np.array(wfn.epsilon_a_subset("AO", "OCC")).astype(str))
coeffs = list(np.array(wfn.Ca_subset("AO", "OCC")).astype(str))
length = int(math.ceil(len(evals)/6))
rem = len(evals)%6
start = 0
end = 0
for i in range(length+1):
    if rem == 0 and i == length: break
    if rem != 0 and i == length:
        start = end
        end = start + rem
    else:
        start = i*6
        end = (i+1)*6
    lin = "%12s" % (' ')
    for j in range(start,end): lin+= "%10d" % (j+1)
    print_out(lin+newline)
    lin = " "
    for num in range(start,end): lin+= "%10.3f" % (float(evals[num]))
    print_out("Eigenvalues: "+lin+newline)
    for j in range(len(coeffs)):
        lin = "%3d %5s" % (j+1,ao_labels[j].split()[0])
        lin += "%5s" % (" ".join(ao_labels[j].split()[1:]))
        for num in range(start,end): lin+= "%10.5f" % (float(coeffs[j][num]))
        print_out(lin+newline)
    print_out(newline)

oeprop(wfn, "MULLIKEN_CHARGES")
mol.print_out()
print_variables()