A general method is presented for the calculation of molecular properties to arbitrary order at theKohn–Sham density functional level of theory. The quasienergy and Lagrangian formalisms arecombined to derive response functions and their residues by straightforward differentiation of thequasienergy derivative Lagrangian using the elements of the density matrix in the atomic orbitalrepresentation as variational parameters. Response functions and response equations are expressedin the atomic orbital basis, allowing recent advances in the field of linear-scaling methodology to beused. Time-dependent and static perturbations are treated on an equal footing, and atomic basissets that depend on the applied frequency-dependent perturbations may be used, e.g.,frequency-dependent London atomic orbitals. The 2n+1 rule may be applied if computationallyfavorable, but alternative formulations using higher-order perturbed density matrices are alsoderived. These may be advantageous in order to minimize the number of response equations thatneeds to be solved, for instance, when one of the perturbations has many components, as is the casefor the first-order geometrical derivative of the hyperpolarizability.
