Improvements in boundary element methods are described for solution of reaction field equations that incorporate important dielectric effects of solvation, including influences of volume polarization, into electronic structure calculations on solute properties. Most current implementations assume constant boundary elements on the cavity surface separating solvent from solute, often employing an empirical parameter to enhance slow convergence associated with the treatment of singularities. In this work we describe a scheme for the linear interpolation of boundary elements and the analytic treatment of singularities that improves convergence without the need for any empirical parameter. Another advance is described for isodensity surface triangulation that succeeds even with molecular surfaces having prominent pockets, which cause the failure of previous simpler methods. Numerical examples are presented to demonstrate the efficacy of these new procedures in practice. (C) 2003 American Institute of Physics.