Creating accurate, analytic atom−atom potentials for small organic molecules from first principles can be a time-consuming and computationally intensive task, particularly if we also require them to include explicit polarization terms, which are essential in many systems. We describe how the CamCASP suite of programs can be used to generate such potentials using some of the most accurate electronic structure methods currently applicable. We derive the long-range terms from monomer properties and determine the short-range anisotropy parameters by a novel and robust method based on the iterated stockholder atom approach. Using these techniques, we develop distributed multipole models for the electrostatic, polarization, and dispersion interactions in the pyridine dimer and develop a series of many-body potentials for the pyridine system. Even the simplest of these potentials exhibits root mean square errors of only about 0.6 kJ mol −1 for the low-energy pyridine dimers, significantly surpassing the best empirical potentials. Our best model is shown to support eight stable minima, four of which have not been reported before in the literature. Further, the functional form can be made systematically more elaborate so as to improve the accuracy without a significant increase in the human-time spent in their generation. We investigate the effects of anisotropy, rank of multipoles, and choice of polarizability and dispersion models.