The calculations were performed in the same fashion but with a few minor differences for the illustrative molecular systems (pyrene, protonated LeuEnk, and doubly protonated gramicidin S, see Figure ). For pyrene, in addition to harmonic spectra at all 27 chosen methods, the intensities of the overtones and combination bands at the r 2 SCAN-3c level of theory were computed with the "NearIR" procedure of ORCA. For protonated LeuEnk, only GGA and 3c-methods were tested. Nine representative conformers, given in Ref. , were optimized at the corresponding level of theory, and then their IR harmonic spectra were averaged with the weights provided in Ref. . For the doubly protonated gramicidin S, the given conformer from Ref. was used after re-optimization at several levels of theory from our list of methods. The IR spectra for these three systems were convoluted with the Gaussian functions with full width at half-maximum FWHM=30 cm -1 for pyrene, FWHM=110 cm -1 for protonated LeuEnk, and FWHM=10 cm -1 for gramicidin S, to match the widths observed in the experimental spectra.