Continuum modeling of electrostatic interactions based upon numerical solutions of the Poisson-Boltzmann equation has been widely adopted in biomolecular applications. To extend their applications to molecular dynamics and energy minimization, accurate and robust methodologies to compute solvation energies and forces must be developed. Indeed, energy conservation is still not possible in molecular dynamics simulations with any numerical Poisson-Boltzmann methods at realistic spatial discretization resolutions of 1/4 to 1/2 Angstrom. It was observed a while ago that the use of surface polarization charges in the computation of solvation energies may greatly enhance the energy convergence in the finite-difference solution of the Poisson equation. In this talk, I will discuss our recent efforts to generalize the strategy to the finite-difference solutions of the full Poisson-Boltzmann equation and to the computation of solvation forces. Our preliminary studies show that a combination of the charge-based strategy, the harmonic average treatment of the dielectric interface, and the charge singularity removal offers the convergence quality in energies and forces on a par with higher-order finite-difference numerical solutions.