The indices µ and ⌫ refer in general basis functions, either electronic or nuclear. The fitting coe cients are determined through a socalled robust fitting described in Ref. [13]. The main objective of this work is to demonstrate how this computationally e cient implementation of NEO-DFT can be applied routinely to estimate anharmonic e↵ects in the zeropoint vibrational energy of hydrogen-bonded molecular clusters. Inherently including nuclear quantum e↵ects (NQEs) like anharmonic zero-point vibrational energy and delocalization within multicomponent DFT carries several advantages. These will be demonstrated on challenging examples whereby regular Born-Oppenheimer DFT is either providing inaccurate or mismatching results. For our study we chose the B3LYP-D3(BJ) electronic functional, the most commonly used functional to date for organic compounds. All of the experiments later mentioned in the text approach the verylow temperature limit and near-vaccuum conditions (or matrix when explicitly noted) so we will always be comparing the latter to 0 K computed enthalpies.