  We propose a harmonic surface mapping algorithm (HSMA) for electrostatic
pairwise sums of an infinite number of image charges. The images are induced by
point sources within a box due to a specific boundary condition which can be
non-periodic. The HSMA first introduces an auxiliary surface such that the
contribution of images outside the surface can be approximated by the
least-squares method using spherical harmonics as basis functions. The
so-called harmonic surface mapping is the procedure to transform the
approximate solution into a surface charge and a surface dipole over the
auxiliary surface, which becomes point images by using numerical integration.
The mapping procedure is independent of the number of the sources and is
considered to have a low complexity. The electrostatic interactions are then
among those charges within the surface and at the integration points, which are
all the form of Coulomb potential and can be accelerated straightforwardly by
the fast multipole method to achieve linear scaling. Numerical calculations of
the Madelung constant of a crystalline lattice, electrostatic energy of ions in
a metallic cavity, and the time performance for large-scale systems show that
the HSMA is accurate and fast, and thus is attractive for many applications.
