, and HCOO -• H + • -OOCH. They all have the shared proton sitting at the center, and their experimental gas phase infrared (IR) spectra have been available through cold-ion vibrational action spectroscopy. For comparison, we calculate their vibrational spectra using both DFT and CNEO-DFT methods. Static energy, gradient, and Hessian calculations are performed using an in-house version 43 of PySCF, and MD simulations are performed with Atomic Simulation Environment (ASE). In CNEO calculations, we only treat hydrogens quantum mechanically because heavy atoms generally have smaller NQEs, although we note that full-quantum calculations are possible. For all the calculations in this communication, if not specified, the aug-cc-pVDZ basis set 47 is adopted for electrons and the PB4-D basis is used for protons. Three electronic exchange-correlation functionals, B3LYP, 49,50 PBE0 51 and ωB97X, are tested in both DFT and CNEO-DFT calculations. Electron-nuclear correlations and nuclear-nuclear correlations are not included. For additional comparison, DFT-VPT2 calculations are performed using Gaussian 16. Since DFT-ωB97X has the best agreement with CCSD(T) in terms of harmonic frequencies for the proton transfer modes in H 5 O + 2 and H 3 O - 2 , we adopt ωB97X for MD simulations with the minimal error from the choice of the electronic functional. In principle, simulating IR spectra requires an ensemble average of a large number of N V E trajectories with starting configurations sampled from an equilibrated N V T trajectory. To reduce the calculation cost, we follow the practice established by existing literature and use a single 20 ps N V T trajectory at 10 K to simulate the IR spectra. We set the initial positions and velocities of each atom following the equipartition law at 10 K with equal amount of energy added to each vibrational mode obtained from Hessian analysis. In practical simulations, a time step of 0.5 fs is chosen and the Nose-Hoover thermostat with a large characteristic timescale being 1 ps is employed. In the end, the dipole autocorrelation function is calculated from each trajectory and then Fourier transformed to obtain the IR spectrum of each system. Results and discussions