Sounds like a job for the JK class.
# If you want to run the job with symmetry but get J and K without, you'll need to build and recompute after the HF.
h2o = psi4.geometry("""
H 1 0.96
H 1 0.96 2 104.5
# If you don't include this command, Psi4 won't know to save the matrices!
eng, wfn = psi4.energy("scf/sto-3g", molecule=h2o, return_wfn=True)
# Calling J() returns a list for some reason I don't understand, but there's only one element in it,
# so that's the one we want. K probably will as well.