What is the best method for calculating partial charges and protonation/deprotonation states of a small molecule such as phenol or caffeine, depending on the pH? For instance, it is well known that above pH=10 or so phenol is deprotonated, but there is data about its protonation around pH=2. For more complex molecules such as caffeine, similar things happen. So I am more interested in calculating partial charges of protonated and deprotonated forms than on pKas (sometimes, experimental data is available).
I head about Gaussian and the like, but I also wonder about software or methods that are easy to use.