the results obtained with NEO-PNO-LMP2 to canonical DF-NEO-MP2 results. The absolute di↵erences are only 0.02, 0.01, 0.01 kcal mol 1 for the Eigen, cis-Zundel and trans-Zundel isomer with the Ring structure as reference. As mentioned before, the inclusion of explicit correlation will reduce this domain error further. Noticeably, NEO-PNO-LMP2-F12 and NEO(MP2)-PNO-LCCSD(T)-F12 are resulting in almost the same energy for the Eigen isomer only di↵ering by 0.05 kcal mol 1 . However, the di↵erence in the energetic ordering for the cis-Zundel and trans-Zundel are more pronounced with a di↵erence of 0.21 kcal mol 1 and 0.24 kcal mol 1 . As shown by the group of Hammes-Schi↵er this di↵erence can be reduced by employing spin-scaling and orbital optimization to the NEO-MP2 method. If compared to the anharmonic corrected energies obtained by the group of Hammes-Schi↵er, the energetic di↵erences to our NEO(MP2)-PNO-LCCSD(T)-F12 method with the Eigen isomer as reference are 0.48, 0.46, 0.27 kcal mol 1 for the Ring, cis-Zundel and trans-Zundel isomer. Compared to the NEO-CCSD method, the inclusion of explicit correlation and scaled perturbed triple excitations for the electronelectron correlation lowers the energies further by 0.52, 0.24, 0.23 kcal mol 1 of the aforementioned isomers. Moreover, we achieve a similar energetic degeneration of the cis-Zundel and trans-Zundel isomer. It should be noted that slight di↵erences in the energies could arise by using di↵erent optimization methods. We optimized the isomers utilizing B2PLYP-D3(BJ) whereas the study of Hammes-Schi↵er and coworkers utilized CCSD(T). However, one advantage of multicomponent approaches is the delocalization of nuclei which renders the exact placement of the basis functions and therefore the underlying optimization di↵erence in the for this system most important quantum mechanically treated nuclei less important. We will further discuss this aspect in other calculations. As another benchmark system we employ the proton-bound dimer of hydrogen sulfate and formate. This is a challenging case for electronic structure methods since it exhibits a shallow potential energy surface with zero-point energy in the order of the proton translocation. Therefore, accurate methods going beyond the Born-Oppenheimer approximation have to be employed to capture nuclear quantum e↵ects and restore a quantitative accurate description of the system. Recently, pathintegral molecular dynamic (PIMD) studies were carried out on the system within a temperature range from 100-300 K. In order to close the gap to the measurement at cryogenic temperature, multicomponent methods like our NEO(MP2)-PNO-LCCSD(T)-F12 approach can be utilized. The great advantage of such approaches is that nuclear delocalization and anharmonic zero-point energies of the quantum nuclei are inherently included. This is computationally much more e cient and robust. Thomas et al. were comparing di↵erent and commonly applied electronic structure methods and found a large spread between the obtained results. In Fig. tirely di↵erent. Under NEO, the global minimum is structure 3 followed by a slightly higher but almost energetically degenerate structure 2. Structure 1 is the most unstable, forming a shallow potential surface. This potential falls in agreement with the PIMD results reported in the literature. In the results from Udagawa et al. it was observed that at a high temperature of 300 K the picture is consistent with the double-hybrid DSDPBEP86-D3(BJ) potential surface (predominance of structure 1, followed by 3). However, by lowering the temperature down to 100 K the probability of observing structure 1 lowers while the populations of structures 2 and 3 increase, consistent with a shallow potential energy surface. This is in line with our NEO(MP2)-PNO-LCCSD(T)-F12 results, which are representative of 0 K. Taking a close look at the respective structures this trend can be further explained. While structure 1 only has one shared proton centered towards the formate, structure 2 exhibits a perfectly shared proton between the hydrogen sulfate and the formate and the energetically lowest structure 3 has two shared protons both centered towards the hydrogen sulfate. In future work it would be interesting to optimize these structures within a multicomponent framework and predict the actual vibrational frequencies accordingly. This would allow us to compare to the experimental values and benchmark the approach. 66 This could be for example achieved by using the NEO-DFT(V) method developed within the group of Hammes-Schi↵er. With this objective in mind, we are currently developing a computational e cient protocol of the NEO-DFT and NEO-TDDFT method within Molpro by employing local density fitting approximations. In the following, we aim to benchmark our NEO(MP2)-PNO-LCCSD(T)-F12 method with other commonly applied electronic structure methods and compare to experimental data. Therefore, we employ as test systems the hydrogen bound methanol complexes with furan (Fu), 2-methylfuran (MFu) and 2,5-dimethylfuran (DMFu) reference systems. Those systems were already thoroughly investigated by both theory and experiment as part of a blind challenge. We optimized the OH-O t , OH-O p and OH-⇡ structures of the respective dimers with DSDPBEP86-D3(BJ), B2PLYP-D3(BJ) and B3LYP-D3(BJ). The obtained structures are shown in Fig. and Fig. . The zero-point corrected energy differences between the most stable OH-O bound structures with respect to the OH-⇡ structure are shown in Fig. together with the respective NEO(MP2)-PNO-LCCSD(T)-F12 energies. Thereby, the hydrogen stretching frequency of the quantum mechanically treated hydrogen of the hydroxyl moiety was removed from the zero-point energy to avoid double counting since this contribution is intrinsically included in the NEO energy. The obtained energies are displayed in Tab. S2. We start by discussing the most noticeable feature of the NEO results shown in the graph, compared to the other methods. Although the individual structures di↵er quite substantially from each other by allowing the most important hydrogen to delocalize the di↵erence in energies narrows down. Therefore, NEO(MP2)-PNO-LCCSD(T)-F12 can be used to obtain reliable and robust results relatively independent of the underlying optimization method. Moreover, independently of the used optimization method the NEO(MP2)-PNO-LCCSD(T)-F12 method predicts, as the only method within the ref- erence set, that the OH-O t structure is energetically favorable over the OH-O p and OH-⇡ structures for the tested complexes. This was already suggested by experiment and later on explicitly confirmed for the methanol furan dimer by rotational spectroscopy. In contrast, every other method employed in this work fails to predict the physically correct trend and instead favours the OH-⇡ structure for the methanol furan dimer. One more point to note is that anharmonic corrections applied by multiple groups participating in the blind challenge also lead to inconsistent results regarding the stability of the conformers. This was, in fact, the most critical point of the entire challenge. This highlights the robustness of the NEO(MP2)-PNO-LCCSD(T)-F12 ansatz, given that all obtained NEO results are within the experimentally determined error bar.