The application of the particle–particle—particle–mesh (P3M) method for computing electrostatic interactions in molecular simulations relies on the use of a charge-shaping function to split the potential into two contributions, evaluated in real and reciprocal space, respectively. Although the charge-shaping function is traditionally taken to be a Gaussian, many other choices are possible. In the present study, we investigate the accuracy of the P3M method employing, as charge-shaping functions, polynomials truncated to a finite spacial range (TP functions). We first discuss and test analytical estimates of the P3M root-mean-square force error for both types of shaping functions. These estimates are then used to find the optimal values of the free parameters defining the two types of charge-shaping function (width of the Gaussian or coefficients of the TP function). Finally, we compare the accuracy properties of these optimized functions, using both analytical estimates and numerical results for a model ionic system. It is concluded that the use of specific TP functions instead of the traditional Gaussian function leads to improvements in terms of computational speed, simplicity of use, and accuracy of results.