A commonly used approximation in the study of vibrational spectroscopies is the double-harmonic approximation,20 where molecular vibrations are described as harmonic oscillators and where the fundamental properties describing the spectroscopic intensities are determined by the first-order geometric derivative of the polarization property governing the spectroscopic phenomenon under study. This corresponds to a description which obeys the well-known selection rules for e.g. infrared (IR) and Raman spectroscopies,21 where it is the magnitude of the first-order geometric derivative of the molecular electric dipole moment and polarizability, respectively, that determines the spectroscopic intensity associated with the excitation to a singly excited state of a particular normal mode. Although coupled-cluster theory can provide vibrational frequencies of high accuracy,22–28 its computational cost prevents its routine use for larger molecules. For this reason, density-functional theory (DFT) has been gaining increasing popularity in recent years and has been used for calculations of Raman spectra of molecules as large as buckminsterfullerene.29 The calculations have often been done in combination with scaling of the frequencies in order to account for anharmonicities and errors inherent in the exchange–correlation functional used.30,31 A typical computational protocol has been to determine the frequencies of the (harmonic) normal modes by DFT using the B3LYP functional coupled with intensity calculations performed at the HF or DFT level,32–34 the choice of level of theory for the intensity calculations depending on the computational tools available for the calculation of the necessary geometric derivative of the pertinent polarization property.