Potential Energy Curve from Hartree-Fock Calculations I am currently working on some quantum chemical calculations, in particular looking at dissociation energy curves generated from Hartree-Fock, CI, CCSD and MP2 methods for diatomic molecules. Now, I was always under the impression that in order to get a potential energy surface  (in the Born-Oppenheimer approximation)  one needs the energies for the electronic problem and the nuclear-nuclear interaction:
$$E(R) + V_{\mathrm{NN}}(R) = U(R)$$
Where $E(R)$ is the energy obtained from the solution of the electronic Schrödinger equation, $V_{\mathrm{NN}}(R)$ the internuclear coulomb potential and $U(R)$ the potential energy surface. I have now seen multiple publications and books that draw potential energy curves using only $E(R)$.
So my question: Are these publications/books wrong, or do I have a misunderstanding of how the PES is computed under the BOA?
Example publication: Dutta, A.; Sherrill, C. D. Full configuration interaction potential energy curves for breaking bonds to hydrogen: An assessment of single-reference correlation methods. J. Chem. Phys. 2003, 118 (4), 1610–1619. DOI: 10.1063/1.1531658.
The energy calculated by a quantum chemistry program is essential always the total potential energy $U(R)$.  For example Fig 1 and Fig 4 of Dutta and Sherrill referenced in the question both plot the total sum of electronic and nuclear energies $U(R)=E(R)+V_{NN}(R)$.
For an example, have a look at a PES I have generated for $\ce{H2}$, using Restricted Hartree-Fock and the STO-3G basis set using PySCF.  The characteristic PES, with a single minimum is actually the total sum of energies $U(R)$, labelled the &quot;Potential Energy Surface&quot;

N.B. The zero bond-length electronic energy is calculated using a helium nucleus as the program doesn't allow superposition of atoms.