Quantum mechanics/molecular mechanics (QM/MM) calculations are applied for anharmonic vibrational analyses of biomolecules and solvated molecules. The QM/MM method is implemented into a molecular dynamics (MD) program, GENESIS, by interfacing with external electronic structure programs. Following the geometry optimization and the harmonic normal-mode analysis based on a partial Hessian, the anharmonic potential energy surface (PES) is generated from QM/MM energies and gradients calculated at grid points. The PES is used for vibrational self-consistent field (VSCF) and post-VSCF calculations to compute the vibrational spectrum. The method is first applied to a phosphate ion in solution. With both the ion and neighboring water molecules taken as a QM region, IR spectra of representative hydration structures are calculated by the second-order vibrational quasi-degenerate perturbation theory (VQDPT2) at the level of B3LYP/cc-pVTZ and TIP3P force field. A weight-average of IR spectra over the structures reproduces the experimental spectrum with a mean absolute deviation of 16 cm-1. Then, the method is applied to an enzyme, P450 nitric oxide reductase (P450nor), with the NO molecule bound to a ferric (FeIII) heme. Starting from snapshot structures obtained from MD simulations of P450nor in solution, QM/MM calculations have been carried out at the level of B3LYP-D3/def2-SVP(D). The spin state of FeIII(NO) is likely a closed-shell singlet state based on a ratio of N-O and Fe-NO stretching frequencies (νN-O and νFe-NO) calculated for closed- and open-shell singlet states. The calculated νN-O and νFe-NO overestimate the experimental ones by 120 and 75 cm-1, respectively. The electronic structure and solvation of FeIII(NO) affect the structure around the heme of P450nor leading to an increase in νN-O and νFe-NO.