End goal is irrelevant, but here is some code that behaves stochastically:
import psi4
mol = psi4.geometry("""
0 1
O
H 1 0.96
H 1 0.96 2 104.5
symmetry c1
""")
psi4.set_options({"basis" : "sto-3g"})
en, wfn = psi4.energy('scf', return_wfn=True)
basis = wfn.basisset()
basis.print_out()
factory = psi4.core.IntegralFactory(basis)
S = factory.ao_overlap()
S.basis1.print_out()
The output from basis.print_out() is sane. However, the output from S.basis1.print_out() is usually nonsense (and sometimes I only get a segfault). Looking on the C++ side, I don’t understand why the basis set isn’t propagating through correctly.