  We present a novel theory and implementation for computing coupled electronic
and quantal nuclear subsystems on a single potential energy surface, moving
beyond the standard Born-Oppenheimer (BO) separation of nuclei and electrons.
We formulate an exact self-consistent nucleus-electron embedding potential from
the single product molecular wavefunction, and demonstrate that the fundamental
behavior of correlated nucleus-electron can be computed for mean-field
electrons that are responsive to a quantal anharmonic vibration of selected
nuclei in a discrete variable representation. Geometric gauge choices are
discussed and necessary for formulating energy invariant biorthogonal
electronic equations. 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 non-classical
nucleus-electron couplings for revising molecular bonding properties, but also
provides an alternative time-independent approach for deploying non-BO
molecular quantum chemistry.
