Electrostatic forces play a crucial role in biomolecular interactions such as protein-ligand association, protein induced membrane deformation or DNA curvature and deformation in protein-DNA complexes. Although the electrostatic forces derived from the potential solution of the Poisson-Boltzmann equation (PBE) have been widely used in studying these interactions with Brownian dynamics simulations, the modeling and computation of these forces in the continuum framework are not well established. In this talk I will propose the models of the electrostatic body force and surface force compatible with the PBE, and discuss the stable computation methods of these forces using a new stable regularization scheme of the PBE. Numerical experiments will be presented to verify these models and their usefulness of in studying the rigid body motion of biomolecules.