  This paper presents a re-formulation of the boundary integral method (BIM)
for the Debye-Huckel model of molecular and colloidal electrostatics that
removes the mathematical singularities that have been accepted as an intrinsic
part of the conventional BIM. The essence of the present boundary regularized
integral equation formulation consists of subtracting a known solution from the
conventional BIM in such a way as to cancel out the singularities associated
with the Green's function. This approach better reflects the non-singular
physical behavior of the systems on boundaries with the benefits of (i) the
surface integrals can be evaluated accurately using quadrature without any need
to devise special numerical integration procedures, (ii) being able to use
quadratic or spline function surface elements to represent the surface more
accurately and the variation of the functions within each element is
represented to a consistent level of precision by appropriate interpolation
functions, (iii) being able to calculate electric fields, even at boundaries,
accurately and directly from the potential without having to solve
hypersingular integral equations and this imparts high precision in calculating
the Maxwell stress tensor and consequently, intermolecular or colloidal forces,
(iv) a reliable way to handle geometric configurations in which different parts
of the boundary can be very close together without being affected by numerical
instabilities, therefore potentials, fields and forces between surfaces can be
found accurately at surface separations down to near contact, and (v) having
the simplicity of a formulation that does not require complex algorithms to
handle singularities will result in significant savings in coding effort and in
the reduction of opportunities for coding errors. These advantages are
illustrated using examples drawn from molecular and colloidal electrostatics.
