It looks like the Coulomb and Exchange energies are added together almost instantly upon computation. You’d either have to modify the compute_E functions in libscf_solver to store those contributions to variables or get the electron repulsion integral tensor and the orbital coefficients and manually compute the Coulomb and exchange components, Psi4Numpy style.

Thanks to both of you. It’s peculiar that our output file contains the “exchange gradient” but no one seems to have ever wanted to know the exchange energy. Of course I was wondering if this data was already printed or stored somewhere. I will implement something like you suggest. Thanks again.