  Nuclear quantum effects such as zero-point energy and hydrogen tunnelling
play a central role in many biological and chemical processes. The
nuclear-electronic orbital (NEO) approach captures these effects by treating
selected nuclei quantum mechanically on the same footing as electrons. On
classical computers, the resources required for an exact solution of NEO-based
models grow exponentially with system size. By contrast, quantum computers
offer a means of solving this problem with polynomial scaling. However, due to
the limitations of current quantum devices, NEO simulations are confined to the
smallest systems described by minimal basis sets whereas realistic simulations
beyond the Born--Oppenheimer approximation require more sophisticated basis
sets. For this purpose, we herein extend a hardware-efficient ADAPT-VQE method
to the NEO framework in the frozen natural orbital (FNO) basis. We demonstrate
on H$_2$ and D$_2$ molecules that the NEO-FNO-ADAPT-VQE method reduces the CNOT
count by a several orders of magnitude relative to the NEO unitary coupled
cluster method with singles and doubles (NEO-UCCSD) while maintaining the
desired accuracy. This extreme reduction in the CNOT gate count is sufficient
to permit practical computations employing the NEO method -- an important step
toward accurate simulations involving non-classical nuclei and
non--Born--Oppenheimer effects on near-term quantum devices. We further show
that the method can capture isotope effects and we demonstrate that inclusion
of correlation energy systematically improves the prediction of $\Delta$ZPE.
