Numerical treatment of singular charges is a grand challenge in solving Poisson-Boltzmann (PB) equation for analyzing electrostatic interactions between the solute biomolecules and the surrounding solvent with ions. For diffuse interface models in which solute and solvent are separated by a smooth boundary, no effective algorithm for singular charges has been developed, because the fundamental solution with a space dependent dielectric function is intractable. In this research work, regularization formulations are introduced to capture the singularity analytically, which are the first of their kind for diffuse interface Poisson's equation and PB models. The success lies in a dual decomposition -- besides decomposing the potential into Coulomb and reaction field components, the dielectric function is also split into a constant base plus space changing part. Using the constant dielectric base, the Coulomb potential is represented analytically via Green's functions. After removing the singularity, the reaction field potential satisfies a regularized PB equation with a smooth source. Some diffuse interface models including a Gaussian convolution surface (GCS) are also introduced. The GCS efficiently generates a diffuse interface for three-dimensional realistic biomolecules. The performance of the proposed regularization is examined by considering both analytical and GCS diffuse interfaces, and compared with the trilinear method. Moreover, the proposed GCS-regularization algorithm is validated by calculating electrostatic free energies for a set of proteins and by estimating salt affinities for seven protein complexes. The results are consistent with experimental data and estimates of sharp interface PB models.