  We develop a parallel rejection algorithm to tackle the problem of low
acceptance in Monte Carlo methods, and apply it to the simulation of the
hopping conduction in Coulomb glasses using Graphics Processing Units, for
which we also parallelize the update of local energies. In two dimensions, our
parallel code achieves speedups of up to two orders of magnitude in computing
time over an equivalent serial code. We find numerical evidence of a scaling
relation for the relaxation of the conductivity at different temperatures.
