A general numerical method is presented to compute the electric potential for a macromolecule of arbitrary shape in a solvent with nonzero ionic strength. The model is based on a continuum description of the dielectric and screening properties of the system, which consists of a bounded internal region with discrete charges and an infinite external region. The potential obeys the Poisson equation in the internal region and the linearized Poisson-Boltzmann equation in the external region, coupled through appropriate boundary conditions. It is shown how this three-dimensional problem can be presented as a pair of coupled integral equations for the potential and the normal component of the electric field at the dielectric interface. These equations can be solved by a straightforward application of boundary element techniques. The solution involves the decomposition of a matrix that depends only on the geometry of the surface and not on the positions of the charges. With this approach the number of unknowns is reduced by an order of magnitude with respect to the usual finite difference methods. Special attention is given to the numerical inaccuracies resulting from charges which are located close to the interface; an adapted formulation is given for that case. The method is tested both for a spherical geometry, for which an exact solution is available, and for a realistic problem, for which a finite difference solution and experimental verification is available. The latter concerns the shift in acid strength (pK-values) of histidines in the copper-containing protein azurin on oxidation of the copper, for various values of the ionic strength. A general method is given to triangulate a macromolecular surface. The possibility is discussed to use the method presented here for a correct treatment of long-range electrostatic interactions in simulations of solvated macromolecules, which form an essential part of correct potentials of mean force. © 1991.