An implementation of a local, linear scaling algorithm for computing electrostatic interactions in molecular dynamics simulations that was recently proposed [J. Rottler and A. C. Maggs, Phys. Rev. Lett. 93, 170201 (2004)] is presented. Coulombic forces are mediated by a propagating electric field that obeys Gauss's law. A detailed description of the numerical implementation is presented, and the accuracy is calibrated by comparing the electrostatic force to that obtained from a standard Ewald summation. A parallelized implementation exhibits excellent scaling behavior on a contemporary low latency compute cluster, and therefore becomes an interesting alternative to existing algorithms for computing electrostatics in large-scale molecular dynamics. The approach opens new possibilities for the study of physical situations that are difficult or impossible to treat with standard reciprocal space methods, in particular nonperiodic boundary conditions and spatially heterogeneous dielectric environments in implicit solvent models.