A procedure is presented that allows us to simulate from first principles the normalized spectra of nuclear inelastic scattering ͑NIS͒ of synchrotron radiation by molecular crystals containing a Mössbauer isotope. Neglecting intermolecular vibrations the NIS spectrum is derived from the normal modes of the free molecule, that are calculated with the density-functional method B3LYP. At low temperatures the inelastic part of the calculated NIS spectrum is a superposition of peaks that correspond to the individual vibrational modes of the molecule. The area of each peak is proportional to that part of the mean-square displacement of the Mössbauer isotope that is due to the corresponding vibrational mode. Angular-dependent NIS spectra have been recorded for a guanidinium nitroprusside single crystal and temperature-dependent NIS spectra for the spin-crossover system ͓Fe(tpa)(NCS) 2 ͔ ͓tpaϭtris͑2-pyridylmethyl͒amine͔. Qualitative agreement is achieved between measured and simulated spectra for different crystal orientations of guanidinium nitroprusside. A remarkable increase of the iron-ligand bond stretching upon spin crossover has unambiguously been identified by comparing the measured NIS spectra of ͓Fe(tpa)(NCS) 2 ͔ with the theoretical simulations.