Thanks for your reply. (Unfortunately) I need arbitrary order multipole integrals.
Based on your suggestion I tried the following:
wfn = psi4.core.Wavefunction.build(molecule, psi4.core.get_global_option(‘basis’))
mints = psi4.core.MintsHelper(wfn.basisset())
Then I can access the integral factory as mints.integral() and do something like mints_be.integral().ao_multipoles(2).compute_shell(0, 0), but how do I then extract the value?