In theoretical investigations, a significant challenge for the study of proton transfer is the incorporation of nuclear quantum effects (NQEs), given that protons are the lightest nuclei with strong NQEs. NQEs often include zero-point effects and tunneling effects, both of which are related to the quantum delocalized nature of nuclei. NQEs are closely related to many anomalous properties of water systems and some other hydrogen-bonded systems. For instance, the quantum delocalization of protons in the active site of the enzyme ketosteroid isomerase (KSI) significantly stabilizes the deprotonated form of an active-site tyrosine residue, resulting in its anomalous acidity. Additionally, NQEs can be manifested in vibrational spectra, where a good description of NQEs will naturally describe anharmonicity, which is significant in many bonds vibrations involving hydrogen atoms. There have been many methods developed to incorporate NQEs and anharmonicities in theoretical calculations. Based on high-dimensional potential energy surfaces (PESs), both static calculations using vibrational self-consistent field/vibrational configuration interaction (VSCF/VCI) and dynamic simulations using multi-configuration time-dependent Hartree (MCTDH) can give highly accurate results. More computationally affordable methods are vibrational perturbation theory (VPT2) and path-integral based simulation methods such as centroid molecular dynamics (CMD) and ring polymer molecular dynamics (RPMD). Compared to VSCF/VCI, VPT2 only needs local higher order derivative information of the PES to perform calculations. Due to the relatively low computational cost, it has been widely used to compute anharmonic corrections with high accuracy. However, VPT2 faces challenges in highly anharmonic systems, for example, the shared proton systems in this communication. Path-integral based methods have slightly higher costs than conventional classical molecular simulations. They have been applied to investigate proton transfer in larger systems such as between phenol and trimethylamine in polar solvents and through two-dimensional materials. In general, path-integral methods are excellent for investigating static properties such as structures and free energies but have limitations in predicting dynamic properties, especially in calculating vibrational spectra, although recent developments has made them more reliable for these calculations. Recently, constrained nuclear-electronic orbital density functional theory (CNEO-DFT) has been developed in our group. It simultaneously treats nuclei and electrons quantum mechanically but with a constraint on each quantum nuclear expectation position. The resulting CNEO-DFT energy surface is a function of quantum nuclear expectation positions as well as classical nuclear positions, and it inherently incorporates NQEs. The CNEO-DFT energy surface has many desirable properties. For instance, the harmonic vibrational frequencies obtained directly from Hessian calculations are comparable to the anharmonic frequencies obtained from DFT-VPT2. Additionally, molecular dynamics can be performed on CNEO-DFT energy surfaces (CNEO-MD), which incorporates NQEs in classical simulations with essentially the same computational cost as conventional DFT-based ab initio molecular dynamics (AIMD). CNEO-MD has been tested and found to be capable of accurately reproducing the key features in the vibrational spectra of a series of small molecules. Without any empirical scaling, CNEO-DFT and CNEO-MD are especially good at predicting the motion of free X-H bonds, where X is a heavy element. Although the low computational cost of CNEO-MD makes it directly applicable to the simulation of proton transfer processes in large-scale systems, the complexity of such systems and the interplay between different factors may make it hard to clearly assess the performance and understand the strengths and weaknesses of CNEO-MD. Fortunately, shared proton systems, in which a proton is shared between two neutral or charged moieties, have been widely investigated both experimentally and computationally. As prototypical patterns for the presence of protons in hydrogen-bonded networks, the shared proton systems can offer valuable insights into the nature of proton transfer processes in different chemical environment and are ideal systems for benchmarking. In water, there are two well-known structures for proton transfer processes: Eigen (e.g. H 9 O + 4 ) and Zundel (e.g. H 5 O + 2 ). In the Eigen structure, the proton is pocessed by one water only and it can be perceived as the starting and ending structure of the proton transfer process, whereas in the Zundel structure, the proton is shared equally by two adjacent water molecules, which is like a transition state for proton transfer. Out of the two, the Zundel structure is especially challenging for pure electronic structure calculations since the shared proton lies on a highly flat and highly anharmonic PES. As a result, not only the conventional harmonic approximation based on electronic structure calculations but also the popular VPT2 method fails to describe the vibration spectra of the Zundel structure, and more expensive methods are needed, such as diffusion Monte Carlo (DMC), MCTDH, and VSCF/VCI. In this communication, we will apply the CNEO method to investigate the vibrational spectra of shared proton systems, or more specifically, the challenging Zundel-like structures: