The QM/DIM method and its analytic energy derivatives are implemented in a locally modified version of the Q-CHEM package. The parameters used in this QM/DIM method are listed in Section S1 of SI. Pyridine (Py) molecule and two gold clusters including 18 and 32 gold atoms are chosen as the test system. Based on these components, three complex configurations shown in Fig. are formed: surficial Py-Au 18 -S, vertical Py-Au 18 -V, and Py-Au 32 with the nearest N-Au distances being 3, 4, and 5 Å. The geometries are obtained from full QM restrained optimization with a harmonic potential, so that each of the structures has only one imaginary frequency along with the labeled N-Au direction; their coordinates are provided in the SI. Restrained optimizations are also carried out by QM/DIM method and the obtained coordinates show no apparent difference from full QM results (with largest RMSD value as 0.007 Bohr in Table of SI). The IR and Raman spectra are calculated by full QM and hybrid QM/DIM methods. In the latter, pyridine is treated by QM and the gold atoms form the DIM region. The partial Hessian is used in the QM/DIM method, whose potential errors are discussed in Section S3 of SI. To validate the implementation of the analytic derivatives of the QM/DIM method, analytical and finitedifferent approaches are used to compute the frequencies and Raman spectral intensities of the Py-Au 18 -V with the labeled N-Au distance as 3 Å and the values are collected in Tab. S7 of the SI. Furthermore, to demonstrate the capability of current analytic derivatives approaches within QM/DIM, we calculate the normal surface enhanced Raman scattering spectra (SERS) of 4,4 -bipyridine (4,4 -BPy), which is set close to the surfaces of Au 2057 and Ag 2057 icosahedral MNPs, respectively (their geometric arrangements are shown in Fig. of SI). The simulated SERS spectra are compared with experimental ones. All the calculations are carried out with PBE0 functional and 6-31+G(d) for N, C, and H atoms and LanL2DZ ECP and basis set for Au atoms using the locally modified Q-CHEM 5.2 package. B3LYP functional is also used for 4,4 -BPy system for the comparison with previous theoretical works. The convergence thresholds of SCF and CPKS are both 10 -8 au while a 10 -14 au threshold is used for two-electron integrals. SG-1 DFT grid 79 is utilized for DFT numerical integration. The vibrational spectra are plotted using a Lorentz function with a width of 10 cm -1 .