The CASPT2 second-order energy derivatives are computed by the finite difference of the firstorder derivatives calculated using the analytical nuclear gradient program recently developed by the authors. We note in passing that the use of the analytical gradient program is essential in realizing the CASPT2 vibrational analysis because it reduces the computational cost by two-to three-orders of magnitude for the molecules in this study. An embarrassingly parallel algorithm has been included that allows the user to select the number of MPI processes required for each gradient calculation and distribute the 6N displacements among different tasks. The Hessian is symmetrized, mass weighted, and the translational and rotational degrees of freedom are set to zero using the usual projection (see Ref. 30 for a concise summary). The normal modes are obtained by diagonalizing the projected mass-weighted Hessian yielding the normal modes (eigenvectors), and the corresponding harmonic vibrational frequencies can be calculated from the eigenvalues. The infrared intensity is also computed from the transition dipole moments. Raman intensities are not yet implemented.