wfn.molecule().nuclear_repulsion_energy() will give you back the nuclear energy in Hartree.

The other thing to note is that the orbitals are one iteration in advance of the current J/K matrices so your energy will differ slightly (usually within the current delta E).

Instead of np.trace(np.dot(D , H)) you can call np.vdot(D, H) for an inner product without trace requirements.

I ran the input provided by hokru above and the energies did not match exactly.
The printout is:
reference : -75.9607980647
HF energy : -75.9607971762
The output file matches the reference energy:
Total Energy = -75.9607980647016205

As noted above the J and K matrices are for the second to last iteration. We always diagonalize the Fock matrices and update the orbital during the last iteration, but do not recompute the J/K matrices. So there will always be a different in the example Hokru gave.

The only way to compute the exact same energy is to recompute the J/K matrices from the current Density matrix.

@jmisiewicz 4-orders of magnitude is somewhat concerning. The example by @bge6 is about 1.e-7 Hartree which is expected with default convergence metrics.