  We present a novel theory and implementation for computing coupled electronic
and nuclear subsystems by formulating an exact self-consistent nucleus-electron
embedding potential from the single product molecular wavefunction. Moving
beyond the standard Born-Oppenheimer (BO) separation of nuclei and electrons,
we demonstrate that the fundamental significance of non-classical
nucleus-electron correlation can be well captured for mean-field electrons that
are coupled with a quantal anharmonic vibration of selected nuclei. Our method
is further applied to characterize vibrationally averaged molecular bonding
properties of molecular energetics, bond length, protonic and electron density.
Moreover, post-Hartree-Fock electron correlation can be conveniently computed
on the basis of nucleus-electron coupled molecular orbitals, as demonstrated to
correlated models of second-order M{\o}llet-Plesset perturbation and full
configuration interaction theories. Our approach not only accurately quantifies
nucleus-electron coupling fingerprints for molecular bonding properties, but
also provides an alternative time-independent approach for deploying non-BO
molecular quantum chemistry.
