Nuclear potential integrals broken out by nucleus

I have been playing around with the Psi4 python interface (in particular, MintsHelper), assessing whether I can integrate it more smoothly into a new fragment-based method I have been working on.

It seems like what I need is to be able to get the ao_potential integrals broken out by nucleus. That is to say that, if the usual matrix contains elements
U_pq = \sum_N <p| Z_N(r_N)^-1 |q>
where r_N is the distance from the electron to nucleus N (and Z_N is its charge), then what I want are matrices with elements
U_Npq = <p| Z_N(r_N)^-1 |q>
such that these matrices would sum up to the usual nuclear potential integrals.

Is this possible?

Perhaps another way of phrasing the question in broader strokes, which might open additional paths to a solution would be: Is it possible to get the integrals corresponding to the attraction of electrons on one subsystem to the nuclei on a different subsystem (and excluding the attraction to its own nuclei)? If the capability described above exists, then clearly I can build this higher-level capability, but maybe something like this is already built in.

I have so far been using horrible hacks involving multiple calls with different atoms “ghosted.”

Thank you!