A software implementation of analytic geometric derivatives of electron-repulsion integrals up to second order is presented for the modeling of vibrational spectroscopies at the level of first-principles Kohn–Sham density functional theory (DFT). In line with the general goals of the VeloxChem program, it targets efficient execution in high-performance computing environments with a hybrid MPI/OpenMP parallelization model and is based on the technique of automatic C++ code generation for high versatility. Gradient calculations scale identically with conventional Fock matrix constructions, and also with the prefactor taken into account, the computational cost of the gradient is significantly lower than that of the self-consistent field (SCF) optimization of the reference state. The Hessian calculation shows a scaling of N3.5 with N being the number of contracted Gaussian basis functions. The computational bottleneck in the Hessian calculation is the solving of the coupled-perturbed Kohn–Sham equations that with VeloxChem can be offloaded to GPU-accelerated nodes. The large-scale virtues of the presented software module are demonstrated by the DFT/B3LYP calculation of the IR spectrum of the entire ubiquitin protein with 1,152 atoms in the quantum mechanical (QM) region and TIP3P water in the molecular mechanics (MM) region. The simulated amide I band shows to be in excellent agreement with experiment.