All NEO calculations have been carried out with Molpro, employing the B3LYP functional including the D3 dispersion correction with Becke-Johnson damping. The (local) density fitted version of the NEO-DFT module integrated in Molpro is an extension of our previously presented (local) density fitted NEO restricted HF implementation. Thereby, we employ (local) density fitting within the electronic subsystem and for the Coulomb coupling between quantum mechanical treated nuclei and electrons. The subsystem of the quantum nuclei is treated by an integral-direct implementation. Overall, the general implementation together with accuracy and performance of density fitting within multicomponent DFT was previously discussed by Mejía-Rodríguez et al. For the quantum mechanical nuclei the PB4-F2 nuclear basis set together with the even tempered 10s10p10d10f fitting basis set with exponents ranging from 2 p 2 to 64 are employed. Both NEO and regular DFT calculations are carried out employing the def2-TZVP, or def2-TZVPP in the case of the deprotonated formic acid trimer, basis sets with the def2-QZVPP-JKFIT density fitting basis set. A threshold of 10 7 a.u. for the energy di↵erence within the electronic and nuclear SCF computations, the di↵erence in the density between iterations and the gradient of the respective nuclear and electronic subiterations. The overall energy di↵erence in the NEO-DFT iterations was set to a threshold of 10 6 Hartree. All Molpro computations employ the direct inversion in the iterative subspace starting after the first iteration with a maximum of 10 Fock matrices as basis to extrapolate. In general the standard grid 3 is employed for the computations, whereas the formic acid clusters are computed with the standard grid 2. The electron-proton correlation is computed with the epc-17.2 functional. In order to assess the error of the density fitting, reference calculations for the formic acid clusters with regular NEO-DFT are carried out with Q-Chem 6.2. Thereby, the standard grid 2, a threshold of 10 8 Hartree for the energy and the geometric direct minimization algorithm were employed. For those systems also the threshold within Molpro was raised to 10 8 a.u. for the energy di↵erence within the NEO microiterations and 10 7 for the overall energy di↵erence. These tighter thresholds were also employed for the NEO-RHF reference wave function of the NEO(MP2)-PNO-LCCSD(T)-F12 method. Those calculations employ the cc-pVTZ-F12 basis set with the cc-pVQZ-JKFIT density fitting basis for the Fock and the exchange matrices as well as the complementary auxiliary basis set for the resolution of the identity and the cc-pVQZ-MP2FIT density fitting basis set. The F12b energies are obtained with the 3*A(LOC,FIX) ansatz. Moreover, the complementary auxiliary basis set singles correction together with the scaling of the perturbative triples are applied. All systems are optimized with Gaussian16 employing the def2-TZVP basis set with very tight SCF settings, tight optimization thresholds and a superfine grid utilizing B3LYP-D3(BJ). In the case of the deprotonated formic acid trimer the def2-TZVPP basis set is employed. Corresponding frequency calculations are carried out with the same settings. In order to obtain the zero-point vibrational energies of the systems without the contributions of the quantum mechanical treated protons, the isotope mass of the respective nuclear centers was set to 9.9•10 12 a.u., making those centers infinitely heavy. All nuclear densities shown are displayed at a 0.02 contour level generated with the PyMOL 2.5.2 program.