We studied the effect of using different heme group charge parametrization methods and schemes (Merz-Kollman, CHELPG, and single- and multiconformational RESP) on the quality of the results produced by the constant-(pH,E) MD method, applied to the redox titration of Desulfovibrio vulgaris Hildenborough cytochrome c(3). These new and more accurate charge sets enabled us to overcome the previously reported dependence of the method's performance on the dielectric constant, ε, assigned to the protein region. In particular, we found a systematic, clear shift of the E(mod) toward more negative values than those previously reported, in agreement with an electrostatics based reasoning. The simulations showed strong coupling between protonating/redox sites. We were also able to capture significant direct and, especially, indirect interactions between hemes, such as those mediated by histidine 67. Our results highlight the importance of having a good quantum description of the system before deriving atomic partial charges for classic force fields.