The screened Coulomb (Yukawa or Debye-Hückel) potential, = exp(-κr )/r , where r is the separation distance and κ is the Debye-Hückel screening parameter, gives a good description of the electrostatic interactions in a variety of biologically and physically important charged systems. It is well known that the direct calculation of the energy and forces due to a collection of N charged particles involves the pairwise summation of all charged particle interactions and exhibits an O(N 2 ) computational complexity which severely restricts maximum problem size. This has prompted the development of fast summation algorithms that allow the electrostatic energy and forces to be obtained in only O(N log N ) operations. To date, however, practically all such implementations have been limited exclusively to pure Coulombic potentials (κ = 0), and the central contribution of the present method is to extend this capability to the entire range of the inverse Debye length, κ ≥ 0. The basic formulation and computational implementation of the spherical modified Bessel function-based multipole expansions appropriate for the screened Coulomb kernel are first presented. Next, a simple model system consisting of a single source charged particle is studied to show that the maximum electrostatic energy error incurred by an M-order multipole expansion for the Yukawa potential is bounded above by the error of the equivalent multipole expansion for the Coulombic potential. Finally, timing and accuracy studies are presented for a variety of charged systems including polyelectrolyte chains, random distributions of charges inside a cube, and face-centered-cubic lattice charge configurations containing up to 103,823 charges.