AbstractWe report here an efficient implementation of the finite difference Poisson–Boltzmann solvent model based on the Modified Incomplete Cholsky Conjugate Gradient algorithm, which gives rather impressive performance for both static and dynamic systems. This is achieved by implementing the algorithm with Eisenstat's two optimizations, utilizing the electrostatic update in simulations, and applying prudent approximations, including: relaxing the convergence criterion, not updating Poisson–Boltzmann‐related forces every step, and using electrostatic focusing. It is also possible to markedly accelerate the supporting routines that are used to set up the calculations and to obtain energies and forces. The resulting finite difference Poisson–Boltzmann method delivers efficiency comparable to the distance‐dependent dielectric model for a system tested, HIV Protease, making it a strong candidate for solution‐phase molecular dynamics simulations. Further, the finite difference method includes all intrasolute electrostatic interactions, whereas the distance dependent dielectric calculations use a 15‐Å cutoff. The speed of our numerical finite difference method is comparable to that of the pair‐wise Generalized Born approximation to the Poisson–Boltzmann method. © 2002 Wiley Periodicals, Inc. J Comput Chem 23: 1244–1253, 2002