For the IR and (regular) Raman spectroscopies, programs that allow for the analytic calculation of the required first-order geometric derivatives of the dipole moment and polarizability, respectively, have been available for some time.49–51 The calculation of anharmonic corrections to the intensities in these spectroscopies requires both the development of the necessary vibrational perturbation theory38 to obtain expressions for these corrections and the possibility of calculating second- and third-order geometric polarization property derivatives, as well as the cubic and quartic force constants, that enter into these expressions. Programs that would allow for the analytic calculation of some of these properties are available, but such calculations have mainly been restricted to the HF level of theory, and for some of the properties (and more so if a DFT description is desired), the researcher has had to resort to numerical differentiation. Analytic calculation offers several advantages over numerical methods such as higher attainable accuracy and ease of computation,51 as numerical derivatives are sensitive to the finite perturbation/geometry displacements employed, and this can have significant effects on the results if not managed carefully.52–54 For these reasons, analytic methods are preferred over numerical ones.