  Density-functional perturbation theory (DFPT) is nowadays the method of
choice for the accurate computation of linear and non-linear response
properties of materials from first principles. A notable advantage of DFPT over
alternative approaches is the possibility of treating incommensurate lattice
distortions with an arbitrary wavevector, ${\bf q}$, at essentially the same
computational cost as the lattice-periodic case. Here we show that ${\bf q}$
can be formally treated as a perturbation parameter, and used in conjunction
with established results of perturbation theory (e.g. the "2n+1" theorem) to
perform a long-wave expansion of an arbitrary response function in powers of
the wavevector components. This provides a powerful, general framework to
accessing a wide range of spatial dispersion effects that were formerly
difficult to calculate by means of first-principles electronic-structure
methods. In particular, the physical response to the spatial gradient of any
external field can now be calculated at negligible cost, by using the response
functions to $\mathit{uniform}$ perturbations (electric, magnetic or strain
fields) as the only input. We demonstrate our method by calculating the
flexoelectric and dynamical quadrupole tensors of selected crystalline
insulators and model systems.
