Discrete sums of the form $\sum_{k=1}^N q_k \cdot \exp\left( -\frac{t s_k}{2 \cdot \sigma^2} \right)$ where $\sigma>0$ and $q_1, \dots, q_N$ are real numbers and $s_1, \dots, s_N$ and $t$ are vectors in $R^d$, are frequently encountered in numerical computations across a variety of fields. We describe an algorithm for the evaluation of such sums under periodic boundary conditions, provide a rigorous error analysis, and discuss its implications on the computational cost and choice of parameters. While the algorithm itself was introduced before (and is closely related to a class of algorithms for the evaluation of non-uniform discrete Fourier Transforms), the error analysis and its consequences appear to be novel. We illustrate our results via numerical experiments.