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. 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. 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, discrete dipole approximation, finite difference time domain, and polarizable continuum models. 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 (including its combination with either electronegativity equalization or charge-transfer ), charge-dipole interaction model, discrete interaction model (DIM) and its coordination-dependent variant (cd-DIM), and atomic dipole approximation model. 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. 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 to model MNP-mediated one-photon and two-photon absorption, Raman, and fluorescence spectra. 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. 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. 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. Here, inspired by the encouraging results of the QM/DIM model from Jensen et al, 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, our implementation is applicable to polarizable force fields involving both induced charges and dipoles, such as the DIM model for metal clusters.