Localized-orbital computation of linear and nonlinear susceptibilities

We present a method to compute high-order derivatives of the total energy which can be used in the framework of density functional theory. We provide a proof of the $2n+1$ theorem for a general class of energy functionals in which the orbitals are not constrained to be orthonormal. Furthermore, by combining this result with a recently introduced Wannier-like representation of the electronic orbitals, we find expressions for the static linear and nonlinear susceptibilities which are much simpler than those obtained by standard perturbative expansions. We test numerically the validity of our approach with a 1D model Hamiltonian.