A detailed presentation of the response theory, which in our approach is fundamental for the analytic calculation of the cubic and quartic force fields and the high-order geometric derivatives of the dipole moment and polarizability that are needed in this work, is too long to show here, and we will therefore restrict ourselves to the most salient features. We refer to the original work48 for a more thorough treatment, and to our recent work47 for a description of the recursive implementation used in the present work.