Continuum solvent models, particularly those based on the Poisson-Boltzmann equation (PBE), are widely used in the studies of biomolecular structures and functions. Existing PBE developments have been mainly focused on how to obtain more accurate and/or more efficient numerical potentials and energies. However to adopt the PBE models for molecular dynamics simulations, a difficulty is how to interpret dielectric boundary forces accurately and efficiently for robust dynamics simulations. This study documents the implementation and analysis of a range of standard fitting schemes, including both one-sided and two-sided methods with both first-order and second-order Taylor expansions, to calculate molecular surface electric fields to facilitate the numerical calculation of dielectric boundary forces. These efforts prompted us to develop an efficient approximated one-dimensional method, which is to fit the surface field one dimension at a time, for biomolecular applications without much compromise in accuracy. We also developed a surface-to-atom force partition scheme given a level set representation of analytical molecular surfaces to facilitate their applications to molecular simulations. Testing of these fitting methods in the dielectric boundary force calculations shows that the second-order methods, including the one-dimensional method, consistently perform among the best in the molecular test cases. Finally, the timing analysis shows the approximated one-dimensional method is far more efficient than standard second-order methods in the PBE force calculations. © 2017 Wiley Periodicals, Inc.