We present an implementation of variational perturbation theory in the framework of density functional theory. We use an ab initio pseudopotential scheme with a plane wave basis set and expand the energy functional up to second order in the perturbation. The approach is fairly general and does not rely on the representativeness of the perturbation through a Hamiltonian operator and does not require the use of canonical orbitals. Instead, a functional formulation is used to characterize the perturbation. Several types of applications are presented which illustrate the variety of linear response phenomena that can be treated with our method (vibrational modes, Raman scattering, and nuclear magnetic resonance chemical shift computations). In combination with advanced gradient correction formulas, an accurate description of second order effects in periodic and isolated systems can be achieved.