In this work, we present a first step towards the efficient implementation of polarizable molecular mechanics force fields with GPU acceleration. The computational bottleneck of such applications is found in the treatment of electrostatics, where higher-order multipoles and a self-consistent treatment of polarization effects are needed. We have coded these sections, for the case of a non-periodic simulation, with the CUDA programming model. Results show a speedup factor of 21 for a single precision GPU implementation, when comparing to the serial CPU version. A discussion of the optimization and parameterization steps is included. Comparison between different graphic cards and a shared memory parallel CPU implementation is also given. The current work demonstrates the potential usefulness of GPU programming in accelerating this field of applications.