Among the most popular approaches used for simulating plasmonic systems, the discrete dipole approximation suffers from poorly scaling volume discretization and limited near-field accuracy. We demonstrate that transformation to a surface integral formulation improves scalability and convergence and provides a flexible geometric approximation allowing, e.g., to investigate the influence of fabrication accuracy. The occurring integrals can be solved quasi-analytically, permitting even rapidly changing fields to be determined arbitrarily close to a scatterer. This insight into the extreme near-field behavior is useful for modeling closely packed particle ensembles and to study "hot spots" in plasmonic nanostructures used for plasmon-enhanced Raman scattering. (C) 2009 Optical Society of America