  Accurate calculation of electrostatic potential and gradient on the molecular
surface is highly desirable for the continuum and hybrid modeling of large
scale deformation of biomolecules in solvent. In this article a new numerical
method is proposed to calculate these quantities on the dielectric interface
from the numerical solutions of the Poisson-Boltzmann equation. Our method
reconstructs a potential field locally in the least square sense on the
polynomial basis enriched with Green's functions, the latter characterize the
Coulomb potential induced by charges near the position of reconstruction. This
enrichment resembles the decomposition of electrostatic potential into singular
Coulomb component and the regular reaction field in the Generalized Born
methods. Numerical experiments demonstrate that the enrichment recovery
produces drastically more accurate and stable potential gradients on molecular
surfaces compared to classical recovery techniques.
