The nuclear-electronic orbital explicitly correlated Hartree-Fock (NEO-XCHF) approach for including electron-proton correlation in mixed nuclear-electronic wavefunctions is presented. A general ansatz for the nuclear-electronic wavefunction that includes explicit dependence on the nuclear-electronic distances with Gaussian-type geminal functions is proposed. Based on this ansatz, expressions are derived for the total energy and the electronic and nuclear Fock operators for multielectron systems. The explicit electron-proton correlation is incorporated directly into the self-consistent-field procedure for optimizing the nuclear-electronic wavefunction. This approach is computationally practical for many-electron systems because only a relatively small number of nuclei are treated quantum mechanically, and only electron-proton correlation is treated explicitly. Electron-electron correlation can be included by combining the NEO-XCHF approach with perturbation theory, density functional theory, and multiconfigurational methods. Previous nuclear-electronic orbital-based methods produce nuclear densities that are too localized, resulting in abnormally high stretching frequencies and inaccuracies in other properties relying on these densities. The application of the NEO-XCHF approach to the [He-H-He](+) model system illustrates that this approach includes the significant electron-proton correlation, thereby leading to an accurate description of the nuclear density. The agreement between the proton densities obtained with the NEO-XCHF and grid-based methods validates the underlying theory and the implementation of the NEO-XCHF method.