Physical realism of molecular dynamics (MD) simulations is greatly enhanced by incorporating variable atomic charges which adapt to the local environment dynamically. In the electrostatic plus (ES+) model, atomic charges are determined to equalize electronegativity. However, this model involves costly minimization of the electrostatic energy at each MD step. A preconditioned conjugate-gradient method is developed for this minimization problem by splitting the Coulomb-interaction matrix into short-and long-range components; the computationally less intensive short-range matrix is used as a preconditioner. This preconditioning scheme is found to speed up the convergence significantly. Numerical tests involving up to 26.5 million atoms are performed on a parallel computer, and the preconditioner is shown to improve the parallel efficiency by increasing data locality. The computational cost is further amortized due to the algorithmic similarity to the multiple-time-scale MD. (C) 1997 Elsevier Science B.V.