We introduce a novel procedure to parametrize biomolecular force fields. We perform finite-temperature quantum mechanics/molecular mechanics (QM/MM) molecular dynamics simulations, with the fragment or moiety that has to be parametrized being included in the QM region. By applying a force-matching algorithm, we derive a force field designed in order to reproduce the steric, electrostatic, and dynamic properties of the QM subsystem. The force field determined in this manner has an accuracy that is comparable to the one of the reference QM/MM calculation, but at a greatly reduced computational cost. This allows calculating quantities that would be prohibitive within a QM/MM approach, such as thermodynamic averages involving slow motions of a protein. The method is tested on three different systems in aqueous solution:  dihydrogenphosphate, glycyl−alanine dipeptide, and a nitrosyl−dicarbonyl complex of technetium(I). Molecular dynamics simulations with the optimized force field show overall excellent performance in reproducing properties such as structures and dipole moments of the solutes as well as their solvation pattern.