I want to run classical molecular dynamics simulations of a periodically replicated surface (rutile TiO2$\ce{TiO2}$ with grooves). In order to do so, I first need to solve for the partial charges residing on each atom of the surface. What electronic structure methods are available for determining the partial charges on a surface?
I am vaguely familiar with Mulliken population analysis, but my understanding is that this method is not the most accurate. What other methods are commonly used, and how accurate are each of these methods?
Finally, what are some good (and free) packages for performing such calculations? I've used SIESTA in the past, is that well-suited for determining partial charges of a surface?