A general answer to the question is just impossible. This answer is limited to modeling the electrostatic term.
We have recently published a detailed comparison of different methods to compute atomic charges, using two sets of 100+ penta-alanine conformers. (link) We tested both the robustness of the charges with respect to conformational changes and the accuracy with which they reproduce the dipole moment (as a simple way of checking the accuracy of the ESP). For the parametrization of a force field, you want a method to assign charges that is successful in both tests.
The Hirshfeld-I method gives the best trade-off between the two tests. However, we also found that for ionic systems (e.g. oxides) Hirshfeld-I overestimates the ESP surrounding isolated molecules, so make sure you double check this for the system you are interested in. (link) If your molecule is rigid, ESP-fitted charges are a better choice because the erratic conformational dependence is no longer an issue.
Anyway, whatever method you use, without including atomic dipoles, one can never parameterize electrostatic interactions up to chemical accuracy, i.e. with errors below 1 kcal mol-1. For rigid molecules, it is possible to make a multipole expansion of each atom-in-molecule density -- with any partitioning scheme such as Hirshfeld-I, Voronio, Mulliken or QTAIM -- to reach any desired level of accuracy for the electrostatic term. (link)