  We recently introduced [J. Chem. Phys. 152 2020, 204103] the
nuclear-electronic all-particle density matrix renormalization group method
(NEAP-DMRG) to solve the molecular Schr\"{o}dinger equation, based on a
stochastically optimized orbital basis, without invoking the Born-Oppenheimer
approximation. In this work, we combine the DMRG with nuclear-electronic
Hartree-Fock (NEHF-DMRG), treating nuclei and electrons on the same footing.
Inter- and intra-species correlations are described within the DMRG without
truncating the excitation degree of the full configuration interaction wave
function. We extend the concept of orbital entanglement and mutual information
to nuclear-electronic wave functions and demonstrate that they are reliable
metrics to detect strong correlation effects. We apply NEHF-DMRG to the
HeHHe$^+$ molecular ion, to obtain accurate proton densities, ground-state
total energies, and vibrational transition frequencies by comparison with
state-of-the-art data obtained with grid-based approaches and modern
configuration interaction methods. For HCN, we improve on the accuracy of the
latter approaches with respect to both ground-state absolute energy and proton
density which is a major challenge for multi-reference nuclear-electronic
state-of-the-art methods.
