We implemented the CNEO-DFT QM/MM in our locally-modified version of PySCF , which is available through our group GitHub page. Molecular dynamics simulations were carried out with GROMACS. In all the following calculations, the aug-cc-pVDZ electronic basis set was used for both CNEO-DFT and conventional DFT. For the glutamic acid-glutamate complex, density fitting with the aug-cc-pVDZ-RI auxilary basis set was used for electronic integrals. With CNEO-DFT, all hydrogen atoms in the QM region were treated as quantum nuclei with the PB4D protonic basis set. The B3LYP electronic exchangecorrelation functional was used and no electron-proton correlation (epc) functional was used. Note that our preliminary tests found that the currently developed epc functionals 100-103 tend to make negligible difference to molecular geometries and vibrational frequencies in the CNEO framework. Therefore, we chose to present results without epc functionals here and leave the detailed investigation of epc effects for future work. For MM-related calculations, we utilized a modified OPLS all-atom force field for the hydrogen-bonded complexes 104-106 and a modified TIP3P water. The polar hydrogens (hydrogen in N-H and O-H) of the hydrogenbonded molecules were assigned Lennard-Jones coefficients from Ref. 108. Additional computational details can be found in the Secion 1 of Supporting Information. Regarding computational cost, we observed that incorporating nuclear quantum effects in CNEO-DFT QM/MM results in only a minor increase in computational cost (typically 10% to 20% longer wall time) compared to conventional DFT-based QM/MM. This comparison is summarized in Section 2 of the Supporting Information.