We present the theory and the implementation of analytical free energy second derivatives with respect to nuclear displacements for a molecular solute described within the framework of the polarizable continuum model. The formulation applies to a cavity with an accurately modeled molecular shape and it permits a complete consideration of all aspects of the solvation model. In particular, the implementation uses the recently proposed method known as the integral equation formalism (IEF), and it can be applied to Hartree-Fock and to density functional calculations. The analysis of both formal and technical features is reported as well as some numerical applications to solvatochromic shifts in IR vibrational frequencies and to transition state searches for reactions in solutions.