We present an implementation of analytical second-order geometric derivatives for density functional methods using Gaussian basis sets. Key features include a stable and efficient numerical quadrature, the direct iterative solution of the coupled perturbed Kohn-Sham equations, integral prescreening based on rigorous estimates, and exploitation of point group symmetry for all finite groups. Benchmark results indicate a moderate cubic growth of CPU and storage requirements with system size; low symmetry molecules with up to 100 heavy atoms can be treated on personal computers. The performance of gradient corrected functionals in predicting IR spectra of larger molecules is exemplified for transition metal carbonyl complexes. (C) 2002 Published by Elsevier Science B.V.