The third-order geometric derivative of the polarizability can be defined from a perturbation tuple (a, b, c, d, e), where perturbations a, b and c correspond to differentiation with respect to geometrical displacements, and d and e to differentiation with respect to a frequency-dependent electric dipole perturbation. Denoting a geometric perturbation as g and the two electric dipole perturbations as f ω and f –ω, where, respectively, each perturbation is associated with a positive or negative frequency ω, eqn (23) takes the form34where the rule choice (k,n) = (2,2) is used because this will give the lowest computational cost. Omitting terms that must be zero straightforwardly or because, in the differentiation carried out, there was lack of dependence on the perturbation operators and, for the sake of brevity, writing contributions that are permutations of identical operators only once, the terms in eqn (34) can be written as35 36 37 38and39where, for example, from eqn (37) is40and where the other differentiated W, Y, and Z terms are of a similar complexity. We consider the length of these expressions, in particular eqn (40), and the corresponding complexity in treating them, as strongly supporting the use of a recursive approach for calculations of the high-order properties required for the GVPT2 treatment, and in a similar manner, automated approaches based on automatic differentiation are needed in order to evaluate the differentiated exchange–correlation energy and kernel .59