Calculations of anharmonic contributions for the purpose of correcting vibrational frequencies have regularly been carried out,17,23,39–41 requiring at least third-order geometric derivatives of the molecular energy, from here on referred to as the cubic force constants. We will refer to the corresponding fourth-order derivatives as the quartic force constants. Calculations of the cubic and quartic force constants have previously almost without exception been done using numerical differentiation.23,41–43 The only exception is the analytic calculation of cubic and quartic force constants at the HF level reported by Handy and coworkers.44,45 Recently, we presented an analytic implementation of cubic and quartic force constants at the DFT level46 by the use of a newly developed recursive code47 for the calculation of molecular properties by response theory.48