Analytic High-Order Energy Derivatives for Metal Nanoparticle-Mediated Infrared and Raman Scattering Spectra within the Framework of Quantum Mechanics/Molecular Mechanics Model with Induced Charges and Dipoles

This work is devoted to deriving and implementing analytic second- and third-order energy derivatives with respect to the nuclear coordinates and external electric ﬁeld within the framework of the hybrid quantum mechanics/molecular mechanics method with induced charges and dipoles (QM/DIM). Using these analytic energy derivatives, one can efﬁ-ciently compute the harmonic vibrational frequencies, infrared (IR) and Raman scattering (RS) spectra of the molecule in the proximity of noble metal clusters/nanoparticles. The validity and accuracy of these analytic implementations are demonstrated by the comparison of results obtained by the ﬁnite-difference method and the analytic approaches, and by the full QM and QM/DIM calculations. The complexes formed by pyridine and two sizes of gold clusters (Au 18 and Au 32 ) at varying intersystem distances of 3, 4, and 5 Å are used as the test systems, and Raman spectra of 4,4 (cid:48) -bipyridine in the proximity of Au 2057 and Ag 2057 metal nanoparticles (MNP) are calculated by QM/DIM method and compared with experimental results as well. We ﬁnd that the QM/DIM model can well reproduce the IR spectra obtained from full QM calculations for all the conﬁgurations, while though it properly enhances some of the vibrational modes, it artiﬁcially overestimates RS spectral intensities of several modes for the systems with very short intersystem distance. We show that this could be improved, however, incorporating the hyperpolarizability of the gold metal cluster in the evaluation of RS intensities. the impact of the and


I. INTRODUCTION
In the last several decades, plasmonic metal nanoparticles (MNPs) have generated a great amount of interest due to the dramatic effect of surface localized plasmon resonance. [1][2][3][4] Currently, plasmonic MNPs have been widely used to control and manipulate light at the nanoscale, thus regulating the photophysical and photochemical properties of molecules in their vicinity, such as enhancing molecular optical signals (absorption, 5 Raman, 6,7 and fluorescence [8][9][10] ) and mediating molecular photochemical reactions. [11][12][13][14] Meanwhile, many theoretical and computational methods were developed to unravel the detailed mechanisms of molecular spectral enhancements, 15 plasmon-mediated chemical reactions, 16,17 as well as plasmon-enhanced resonance energy transfer. [18][19][20][21] In many cases, the plasmon enhancement arises mainly from the electromagnetic mechanism, where an external electric field can be substantially magnified by an MNP around its surface. Accordingly, several classical nonatomistic methods 22 had been widely applied.
The hybrid MNPs-molecular systems represent a challenge for theoretical and computational methods because the cona) Electronic mail: yuezhi.mao@stanford.edu b) Electronic mail: yihan.shao@ou.edu c) Electronic mail: liangwz@xmu.edu.cn stituents can not be treated on the same footing by the stateof-the-art methods. The properties of molecules require a quantum mechanics (QM) description thanks to the high accuracy of QM methods. 23,24 However, QM methods are unable to describe the medium to large-sized MNPs due to their steep computational cost. Two kinds of simplified treatments have been usually adopted. One is to simplify the problem by assuming the molecules bonded to very small metal clusters. [25][26][27][28][29][30][31][32] With this treatment, the optical properties of molecule-metal cluster systems can be described by the QM methods. Small MNPs, however, don't support the bulk plasmon. The other is to adopt the mixed quantum/classical approach, which combines the QM approaches with classical mechanics or electrodynamics, such as the Mie theory, 33 discrete dipole approximation, 34 finite difference time domain, 5 and polarizable continuum models. [35][36][37] To more accurately capture the plasmonic effect from MNPs with different sizes, shapes, and compositions, however, atomistic modeling of the MNP-adsorbate system is needed. Specifically, to account for the large polarizability of MNPs, (induced) charges, dipoles, and/or even multipoles would be introduced at each atom site, leading to several (polarizable) force field models for MNPs. These molecular mechanics (MM) methods include the point-dipole interaction model 38,39 (including its combination with either electronegativity equalization 40 or charge-transfer 41 ), charge-dipole interaction model, 42,43 discrete interaction model (DIM) 44,45 and its coordination-dependent variant (cd-DIM), 46 and atomic dipole approximation model. 47 These MM models can be employed in combined quantum mechanics/molecular mechanics (QM/MM) modeling of probe-MNP complexes, where the probe molecule constitutes the QM region so that its vibrational motions (in infrared and Raman spectroscopy) or electronic transitions (in fluorescence) are treated using ab initio QM theories. Mikkelsen et al have investigated the hyperpolarizability changes of organic molecules near gold nanoparticles by a polarizable QM/MM approach recently. 48,49 The DIM model with induced charges and induced dipoles at each metal atom, in particular, has been combined with QM methods by Jensen and others [50][51][52][53][54] to model MNP-mediated one-photon 51,55 and two-photon absorption, 56 Raman, 52,57-59 and fluorescence spectra. 60 In parallel to the QM/MM study of MNP-enhanced electronic transitions, Giovannini, Cappelli, and others have combined density functional theory (DFT) and time-dependent density functional theory (TDDFT) with their FQFµ models for water molecules, which are also based on fluctuating charges and induced dipoles. 61,62 Analytic energy derivatives can provide computational advantages in molecular geometry optimizations, vibrational frequency calculations, and molecular property descriptions. The analytic Raman intensities of the QM method have been derived and implemented decades ago. [63][64][65] Recently, the implementations of analytic Hessian and high-order molecular properties (such as IR and Raman intensities) have been extended to several hybrid (polarizable) QM/MM schemes for solvents or protein systems. [66][67][68] Here, inspired by the encouraging results of the QM/DIM model from Jensen et al, [50][51][52] we present our implementation of this model and its analytic energy derivatives with respect to the nuclear coordinates and perturbed electric field for the study of the adsorbate-MNP systems. Compared to earlier implementations of IR and Raman intensities within a polarizable QM/MM model, [66][67][68] our implementation is applicable to polarizable force fields involving both induced charges and dipoles, such as the DIM model for metal clusters.
The paper is organized as follows. Section II briefly recaps the basic theoretical foundation of the polarizable QM/MM model for the ground-state energy of a molecule-MNP system and presents the analytical expressions for high-order derivatives of the energy with respect to the QM nuclear and field perturbations. It is followed by the computational details in Section III. In Section IV, using the pyridine molecule in the vicinity of three gold clusters as testing systems, the computed IR and Raman spectra are presented and discussed. In addition, the SERS of 4,4 -bipyridine on gold and silver MNPs are simulated by QM/DIM method and compared with experimental spectra. Moreover, the effects of nonlinear response of MNP and charge migration to Raman intensity are analysed. Finally, the concluding remarks and potential future improvements are summarized in the last section.

II.A. Total Energy of Hybrid QM/MM System
Here the DIM method is briefly described with a compact notation and the reader is referred to the original publications. 44,45 The energy of the MM region is given by the electrostatic interaction, where M ind represents the collection of induced charges and dipoles (q ind , µ ind ) on MM atoms; T the interaction matrices of charge-charge, charge-dipole, and dipole-dipole (T (0) , T (1) , T (2) ) as where r i j is the distance between the i-th and j-th MM atoms, and σ is the width of Gaussian functions used to damp the electrostatic interaction between neighboring atoms. 44 The diagonal block of T is formed by the atomic capacitance (c i ) and polarizability (α i ) for charge and dipole self-interactions at i-th MM atom, respectively. F tot is the total external generalized field exerted on the MM sites, including both potentials and fields (−V, E). Using the variational condition, the stationary solution is given by solving a linear equation, where the Lagrangian multiplier λ is introduced to constrain the net charge for the MM region as Q. 1 refers to a row vector (1, 0, 0, 0, 1, 0, 0, 0, · · · ) with a value of 1 for only the induced charge on each MM atom. When Q is equal to 0 (as for all test cases in this work), the solution to Eq. (3) could be formally written as, where T −1 represents the upper-left 4N MM × 4N MM block of the inverse of the matrix in Eq. (3). Therefore, the MM energy is reformulated as The last expression for the MM energy will be used in the derivation of the high-order derivatives of the energy. When the MM region is adjacent to the QM region (as described on the DFT level of theory), the total generalized field F tot would include the contributions from the QM region in addition to the external source, where the matrix V contains integrals φ µ T φ ν . Namely, for each pair of basis function (φ µ and φ ν ) and the j-th MM atom, it contains the interaction between the basis function pair and a Gaussian MM charge (with width σ ), and its first derivatives.
In Eq. (6), P is the density matrix with elements P µν , Z I labels the I-th QM nuclear charge and T I is the corresponding generalized external field on MM atoms, and R f ext represents the generalized field on each MM atom at position (R jx , R jy , R jz ): With the definition of the total generalized field F tot in Eq. (6), the total electronic energy of the system can then be expressed as (10) where the first term is the well-known DFT core Hamiltonian h 0 and the two-electron integral supermatrix II in atomic basis adopting approximate coefficient C HFX and operator f ( r 12 ) for hybrid and range-separated functionals. Further, E xc = d r f xc ( r) is the DFT exchange-correction energy, and µ QM is the dipole moment of the QM region, which leads to the interaction between QM and external fields within the dipole approximation. Here M is the dipole moment matrix and R I the I-th QM nuclear coordinate. The corresponding Fock matrix to be used in the SCF cycles could be obtained by differentiating the total energy given by Eq. (10) with respect to the density matrix (P): where the elements of F xc can be computed as (14) with variables, ξ , being electron density or its gradient. 69 With this Fock matrix, we could define a total core Hamiltonian as, and then the expression (16) will be used for the calculation of the total energy. The last term is the van der Waals (vdW) interaction between QM and MM regions. Here we adopt the distance-dependent variance of the classical Lennard-Jones (LJ) 12-6 potential according to Jensen et al. 52 The expression is presented in Section S1 of the supporting information (SI).

II.B. Energy Gradient and Hessian Matrix
The first derivatives of the energy in Eq. (16) with respective to the QM coordinates, x, can be expressed as where superscripts with [] refer to the explicit derivatives excluding the contributions from orbital rotations. Here with V . The energy weighted density matrix is The contributions from MM part could be found by differentiating the last two terms in Eq. (10) as where derived from Eq. (6).
The Hessian can then be obtained by directly differentiating the QM gradient in Eq. (17) with respect to the QM coordinates, (22) where with V [xy] . The MM contributions to the Hessian matrix could be summarized as Here P y includes the derivatives of orbital rotations. The first term represents the MM charges and dipoles induced by the seond-order QM potentials and fields with respect to QM nuclear coordinate changes, while the second term is the interaction between first-order QM potentials and fields and MM induced charges and dipoles from another first-order QM potentials and fields perturbed by QM atom displacements. After the construction of the Hessian matrix in the mass-weighted coordinates, and projecting the rotational and translational degrees of freedom out of this harmonic force constant matrix, one could obtain the normal modes and vibrational frequencies by diagonalizing the projected matrix. 70 Note that here the Hessian matrix only includes the QM-QM block, and hence, it is actually partial Hessian within the QM region.

II.C. Infrared and Raman Spectral Intensities
IR and Raman spectra provide information about the molecular vibrations. IR spectra can be obtained from light absorption whereas Raman spectra reflect light scattering process. The IR and Raman spectral intensities are proportional to the squares of nuclear derivatives of molecular electric dipole and polarizability, respectively. 65,70,71 When the perturbation is the external field f m , the derivatives of total energy in Eq. (10) would give the total dipole moment where the first two terms produce the QM dipole moment as in Eq. (12), and R m is the m-th column of matrix R in Eq. (9). For m = 0 (the x component), for instance, with j as the MM atom index. The first-order changes of the total dipole in Eq. (25) with respect to nuclear displacements could be derived as Here the last term refers to the monopolar/dipolar response in the MM region. Namely, as one displaces the QM atoms, it leads to a variation in the corresponding potential and field (F QM,x ) and thus causes a change to the induced charges and dipoles on MM atoms.
The electronic polarizability is obtained by differentiating the total dipole moment given by Eq. (25) with respect to the external field f n , and its nuclear derivative is derived as where F QM,nx · T −1 · R m represents the contribution of the MM part, which indicates that the two perturbations both act on the QM region (i.e., F QM,nx ) in the current framework. It is incorporated into the dipole moment matrix and its derivative matrix in the last line. The derivatives of density matrix are derived in Appendix A. As in the case of the first-order derivatives of density matrix P x or P n , only the first-order derivatives of orbital rotations are needed in the second-order derivatives of density matrix P nx , which is known as the 2n + 1 rule. 72,73 The derivatives of orbital rotations are computed by the coupledperturbed Hartree-Fock equation or z-vector equation, which is presented in Appendix B.

III. COMPUTATIONAL DETAILS
The QM/DIM method and its analytic energy derivatives are implemented in a locally modified version of the Q-CHEM package. 74 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. 1 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 S3 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. S5 of SI). The simulated SERS spectra are compared with experimental ones. 75,76 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. 74 B3LYP functional is also used for 4,4 -BPy system for the comparison with previous theoretical works. 77,78 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 .

IV.A. Infrared Spectra
The infrared (IR) spectra of pyridine molecule and three different pyridine-gold complexes (Py-Au 18 -S, Py-Au 18 -V, Py-Au 32 ) computed by QM and QM/DIM methods are displayed in Fig. 2. The top row shows the IR spectra of the gas-phase pyridine minimum-energy structure, while the other three rows show the spectra calculated using the optimized configurations (with the distances restrained at 3, 4, and 5Å). In these three rows, the green curves correspond to the spectra of the pyridine molecule at the geometry extracted from each of these complexes. Overall, the QM/DIM method (blue curves) reproduces similar trends in the harmonic frequency shifts and comparable IR intensities relative to the full QM profiles (red curves) for the tested systems.
Overall, the harmonic vibrational frequencies of the pyridine molecule are shifted by up to 20 cm −1 from the gas phase to the complexes. Within each of the three complexes, however, the frequencies of the adsorbed pyridine molecule have marginal shifts towards higher frequency (smaller than 7 cm −1 ) upon the binding to Au clusters. This is not surprising because the heavy Au atoms vibrate very slowly and have a negligible effect on the nearby pyridine molecule. For a representative structure, Py-Au 18 -V (3 Å), the numerical values are provided in Table S13 of the SI, and information for several key vibrational modes are collected in Table I. When the pyridine molecule approaches the Au clusters, the IR peaks of three vibrational modes (ring breathing and two C-H wags) within 1000∼1300 cm −1 become increasingly more intense as the distance reduces from 5 to 3 Å. The second C-H wag at 1238 cm −1 is predicted by full QM calculations to have the largest enhancement of 24 times (from 2.7 to 48.7 km/mol) for the S configuration, 33 times for V configuration, and 24 times for Py-Au 32 , as shown in the panels (j), (k), and (l) of Fig. 2, respectively. In contrast, some peaks are significantly weakened after the binding to Au clusters. This is especially obvious with the two out-of-plane C-H bends at 744 and 780 cm −1 , which have comparable intensities in the gas phase. While only the IR peak at 744 cm −1 is reduced in Py-Au 18 -S (3Å) and Py-Au 18 -V (3Å) complexes, both peaks (744 cm −1 and 780 cm −1 ) lose some strength in the Py-Au 32 (3Å) structure.
These changes in the IR intensities are reproduced reasonably well by QM/DIM calculations, as illustrated in Fig. 2. As shown in Table I    When the distance is 5 Å, although QM/DIM slightly underestimates the intensity of mode 12, it nearly reproduces the full QM spectra as Figs. 3(d), 3(e), and 3(f) show. As R decreases, however, QM/DIM gives less satisfying Raman profiles. Taking Py-Au 18 -V system with R = 4 Å as an example (see Fig. 3(h)), we observe that QM/DIM significantly underestimates the intensity of mode 12 but overestimates this of mode 11 compared with the full QM results. When R is further shortened to 3 Å, the discrepancy between full QM and QM/DIM results in all frequency ranges increases. As shown in Table I, QM/DIM significantly overestimates the intensities of all the modes except mode 12. In other words, the QM/DIM approach could provide artificial enhancements for some vibrational modes, which can be attributed to the neglect of chemical enhancement and nonlinear responses of the MNP in the current QM/DIM model. Moreover, we calculate the normal Raman scattering spectra of 4,4 -bipyridine (4,4 -BPy) in the gas phase, and in the proximity of Au 2057 and Ag 2057 with respect to the DFT XC functional PBE0 and B3LYP, respectively. Fig. 4 and Fig. S6 in SI show the calculated results, where the gas-phase Raman intensities are scaled by a factor of 20. It is noted that the spectra calculated with PBE0 and B3LYP are almost identical, indicating that the impact of functional on the spectra is small.
To compare with the experimental off-resonance SERS spectra with the incident light of 535 nm, here the Raman scattering intensity (in unit of cm 2 /sr) for the vibrational mode k is calculated as based on short-time and Placzek approximations. 80 Here ω 0 is the frequency of the incident light and S k (in unit of Å 4 /AMU) refers to the Raman scattering factor, where isotropic and anisotropic polarizability derivatives,ᾱ and γ , are composed by the components of nuclear derivatives of molecular polarizability α as, respectively. The constants in Eq. (29) are provided in Table  S2. Fig. 4 demonstrates that, for the adsorbed configurations, QM/DIM predicts a 20 times enhancement for most of the Raman peaks due to electromagnetic enhancement. Both the theoretical calculation and experimental measurement 75,76 yield a consistent result trend on the MNP-induced changes on the relative spectral intensities. For example, the gold nanoparticle makes the strongest Raman peak of 4,4 -BPy shift to 1678 cm −1 from 1336 cm −1 , the relative intensity at around 1024 cm −1 decrease, and the two Raman peaks in the 600-800 cm −1 range disappear. However, the QM/DIM predictions hardly alter the molecular vibrational frequencies while the experimental frequency shifts could be as large as 18 cm −1 . 75,76 It is unsurprising because after the 4,4 -BPy is adsorbed on the MNP, the QM/DIM optimization yields a similar geometry as the gas-phase one, thus producing only small changes in the vibrational frequencies. This is in line with the Raman spectra of the Py cases. As shown in Table I, the frequencies calculated by full QM and QM/DIM methods have similar values. The large frequency changes after molecule adsorbed on MNP come from geometry differences between Py gas-phase minimum and full QM optimized Py-Au 18 complex.

IV.C. Contribution of Nonlinear Response of MNP to Raman Intensity
To shed light on the issue of the QM/DIM approach in the Raman calculations, we calculate the nuclear derivatives of the polarizability through the finite-difference (FD) method based on classical turning points (CTP). Two vibrational modes are chosen: the ring breathing modes around 1020 and 1050 cm −1 (modes 11 and 12). As shown in Fig. 5, the pyridine nitrogen atom shows little participation in the vibrational pattern of mode 12 while it gets involved in that of mode 11 (the same for other artificially enhanced modes). As Eq. (27) shows, the molecular electric polarizability α mn is related to the density matrix derivatives with respect to the external field and molecular dipole moment matrix. We can divide the value of α mn into four terms with respect to the atomic basis sets contributed by its two components: the Py molecule and Au cluster. As shown in Table II, for mode 11, the Au cluster contributes 8.877 au to the final polarizability change, which largely cancels the other contributions; while the value (-0.519 au) is negligible for mode 12. This distinction could be understood from the vibrational motions of the pyridine. The nitrogen atom plays a major role in the vibration of the ring breathing mode 11 in Fig. 5(a) and changes the local fields exerted on the Au clusters, which requires the contribution from high-order nonlinear response terms, i.e., the hyperpolarizability-related terms of the Au atoms, to be incorporated. On the other hand, as the vibration of ring stretch mode 12 in Fig. 5(b) is symmetric along the chosen z axis, there is no change in the polarizability of the Au cluster. Furthermore, there is no explicit contribution from the metal atoms to the Raman intensities in the current QM/DIM method, which can be revealed by comparing Eqs. (26) and (28).
As stated for Eq. (28), in the current QM/DIM scheme, the contribution of the MM region only comes from the induced charges and dipoles by the second-order potentials and fields of the QM counterpart, where the induced dipole moments are proportional to the first power of the field intensity. MNP may possess a strong nonlinear response, and the induced dipole moments should include the nonlinear terms that are proportional to the second and higher powers of the field intensity. Therefore, it is reasonable to account for the nonlinear response of MNPs.
The incorporation of the hyperpolarizability has been proposed by Andrews et al in their analysis of molecular Ra- where R AB is the inter-distance between two fragments, A and B;ê k is a unit number at k-th direction; µ A,x is the nuclear derivatives of the ground state dipole moment of A; and β β β B represents the electric hyperpolarizability tensor of B. It corresponds to the event that the incident and scattered photons of the Raman scattering are on the same fragment A, meanwhile one virtual photon is exchanged between the two fragments, i.e., excitation energy transfer (EET) occurs between two locally excited states of the two molecules. We note that to arrive at the hyperpolarizability-dipole interaction expression in Eq. (32), two approximations are made for the coupling between the two locally excited states: the exchange integrals are ignored and a multipole expansion of the Coulomb interaction operator is adopted. (The third one is that charge-transfer excitations are not into consideration.) This expression, however, is incomplete because the distance could also be differentiable rather than taking the first-order term in dipole moment (shown in Eq. (37)). On other hand, if the hyperpolarizability β β β MM is introduced in the QM/DIM at the beginning, the additional term of the energy would be whose second-order derivative with respect to external field is Then its nuclear derivative is derived as Here we had assumed f ext = 0 as in the case of no external field. The first line at the right-hand-side in the above expression closely resembles Eq. (32): Provided that the electric field is generated by a dipole as f QM Then the interaction of the first term and the hyperpolarizability yeilds Eq. (32), while the second term is from the derivatives of the distance rather than molecular dipole. The rest contributions to the nuclear derivatives of polarizability α x i j in Eq. (36) includes the expression β MM mnl f QM l , referring to the first-order change in metal's polarizability induced by QM field. As in the QM/DIM method, each metal atom is induced by QM potential and field by atomic capacitance and polarizability and the induced charges and dipoles interact with each other, it wold be difficult to distinguish this mutual polarization from the contributions in Eq. (36) with the use of molecular hyperpolarizability parameter for metal cluster. Therefore, we use Eq. (35) to account the effect of hyperpolarizability of metal in the following. This correction is sufficient for our purpose because the gold cluster has large hyperpolarizability and it has little participation in the vibrational modes of the Raman spectral region we are interested in. The hyperpolarizability of the Au 18 cluster calculated by the finite-difference method with a field strength of 0.001 au is given in Table III. Consistent values are obtained with a smaller field strength (0.0002 au). As demonstrated in Table III, the calculated hyperpolarizability of the Au 18 cluster has large components along the z axis, comparable with the values reported in the literature. 83,84 The corrected QM/DIM Raman intensities of selected modes of Py-Au 18 -V (3 Å) complex is shown in Table IV. There are five variations (a-e) of the hyperpolarizabilityrelated correction for the QM/DIM method since some arbitrariness arises when only the molecular hyperpolarizability parameter is available rather than the atomic ones. The (a), (b), and (c) columns are calculated by taking the distance R AB in Eq. (32) as the distance between the two centers of mass (COM), the COM of pyridine and the location of the nearest gold atom, and the shortest distance of the two fragments, respectively. In other words, the distance R AB reduces gradually when moving from (a) to (c). The (d) and (e) columns collect the Raman intensities corrected using Eq. (35) without and with QM induced fields, respectively. Even though the QM/DIM Raman intensity of mode 11 is still larger than the full QM one, it is unsurprisingly reduced with the correction of the hyperpolarizability dependent term. On the other hand, the intensity of mode 12 is also slightly increased by the (a), (b), and (c) corrections.  For approach (c), even though more correct intensity is reproduced for mode 11, several modes are wrongly enhanced such as 5, 16, and 19. It reflects that the distance is too small within this calculation. When the distance is larger, for instance, in the case of approach (b), the green curve moderately approaches the QM one in Fig. 6. On the other hand, the (e) profile (blue line) improves the intensities of almost all the vibrational modes except that the intensity of mode 16 is overestimated and that of mode 12 is further reduced relative to the uncorrected QM/DIM result (from 236 to 168 Å 4 /AMU, while the full QM reference is 294 Å 4 /AMU). It is acceptable though that these corrections could not produce overall satisfying intensities for all the modes, because only the molecular hyperpolarizability of the gold cluster is employed in these calculations rather than the distributed atomic hyperpolarizabilities while the fields and their derivatives act on individual gold atoms.

IV.D. Population Analysis
Raman signals can be enhanced by the chemical enhancement and the electromagnetic enhancement. In the current QM/DIM scheme, we don't account for the intermolecular charge transfer (CT) effect, which may be partially responsible for the limited accuracy of the QM/DIM scheme for Raman spectral calculations, especially when the intermolecular distance is short. To check this effect, we calculated the amount of charges on pyridine molecule within the complexes as shown in Table V.
Among the four used population methods, only the electrostatic potential (ESP) and fragment-based Hirshfeld (FBH) charge schemes consistently predict the expected trend in the fragment charge population, namely an enhanced CT with a shorter distance between two fragments. These values clearly show that the net transferred charge between Py and Au cluster is almost zero when the nearest distance between these two fragments is 5 Å, while it can be as large as 0.16 e − when the distance shortens to 3 Å in Py-Au 18 -V according to the FBH calculation. It suggests that the ground state charge migration between the Py molecule and Au cluster might also play an important role in the limited accuracy of our current QM/DIM Raman spectral simulations.

V. CONCLUDING REMARKS
In this work, we derived and implemented the QM/DIM method and its high-order energy derivatives for the calculation of noble MNP-mediated molecular vibrational spectra, including IR and Raman scattering. Taking the complexes composed of a pyridine molecule with MNPs consisting of 18 or 32 gold atoms as testing examples, we assessed the accuracy of the current QM/DIM scheme in the calculations of molecular vibrational frequencies and IR/Raman intensities. The QM/DIM method demonstrates the ability to reproduce the molecular vibrational frequencies and IR spectral intensities obtained using the full QM approach. The Raman profiles from QM/DIM calculations behave well for large intermolecular distances such as 5 Å for the tested complexes, while it shows artificial enhancements at shorter molecule-MNP distances for some of the vibrational modes. Even though, with large MNPs such as Au 2057 and Ag 2057 , the simulated normal SERS of 4,4 -BPy by QM/DIM method are comparable with experimental spectra.
The deviation of QM/DIM results from full QM results for small metal clusters, however, can be partially resolved. By incorporating the term that incorporates the first-order nonlinear polarizability of Au clusters and the dipole (or electric field) derivatives of the probe molecule, we can capture the non-negligible three-virtual-photon process of Raman scattering within this hybrid method. With this correction, the Raman spectra for the short-distance configurations of the Py-Au clusters calculated using QM/DIM are improved. We thus suggest that the atomic (distributed) hyperpolarizability of metal atoms should be parameterized within the hybrid QM/DIM method to enable it to describe Raman scat-  tering spectra more accurately. Besides, as investigated by Schatz and the coworker, 89 the other high-order properties such as dipole-quadrupole and dipole-magnetic-dipole polarizabilities may also play a rule in the Raman intensity if the electric field derivatives are comparable with the field strength.
Furthermore, at short intermolecular distances, the overlap of wave functions opens the possibility of a charge migration between Py and the Au cluster. In those cases, the effect of intermolecular CT on the Raman signal should not be neglected, which may otherwise limit the applicability of the current QM/DIM scheme for Raman spectral simulations. It is therefore desirable to develop a scheme to account for the CT effect in future QM/MM methods. The use of only QM-QM block Hessian matrix in the current work is also a potential source of small errors in the simulated vibrational spectra. In addition, the vdW interaction between the molecule and MNP could affect the optimized molecular structure and then its vibrational frequencies predicted by using the QM/DIM method.
Finally, we note that the plasmon resonance effect of MNPs has not been taken into account in this work. We expect to include the resonance effect by adopting the frequencydependent parameters for metal atoms in future publications, which requires solving complex-valued response equations.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available within the article and its supplementary material.

Appendix A: Derivatives of Density Matrix
The occupied MO coefficients could be expanded to the second-order in orbital rotations Θ vo as Ref. 90 suggested.
Its first-order derivative is without second-order terms while the second-order derivative could be further written as where we set Θ vo = 0 as the MOs are optimized. Here When the perturbation (x) in Eq. A2 is the external field in the n direction, it is simplified to Similarly, if the second perturbation (y) in Eq. (A3) is the external field in the n-th direction, it reduces to As the derivatives of density matrix are dependent on the derivatives of occupied MOs and vo ) are determined using the procedure outlined in the next section. where F x = F [x] + F (2) · P x . After moving the P x dependent terms to one side, we arrive at where F (2) = II + Ω Ω Ω + VT −1 V. Note that Ω Ω Ω is the exchangecorrelation portion of the response kernel as defined in Eq.
for QM nuclear perturbations, where the last line at the righthand-side represents the contribution from MM part. And one uses to solve for Θ [n] vo with external field perturbations. The products of property vectors and response properties are interchangeable, for instance, Θ  by interchanging the variables. It corresponds to the wellknown 2n + 1 rule.