  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.
