  We present a formulation of third-order density-functional perturbation
theory which is manifestly invariant with respect to unitary transfomations
within the occupied-states manifold and is particularly suitable for a
practical implementation of the so called `2n+1' theorem. Our implementation is
demonstrated with the calculation of the third-order anharmonic coupling
coefficients for some high-simmetry phonons in Silicon.
