where each term in the Fock-like operator on the left side of the equations is obtained by taking the functional derivative of the corresponding energy term in Eq. 2 with respect to either the electron density or the proton density. In the Fock-like operator for the nuclear equations, an extra term f I • r shows up due to the constraints on the nuclear expectation positions. The equations 7 and 8 are coupled and can be solved iteratively together with the Lagrangian multipliers {f I }. The epc17 family of electron-proton correlation functions was developed analogously to the Lee-Yang-Parr (LYP) electron correlation functional, and more details can be found in Ref. 34. Based on an local density approximation, the electron-proton correlation functional takes the form