We develop an efficient simulation method to study suspensions of charged spherical colloids using the primitive model. In this model, the colloids and the co- and counterions are represented by charged hard spheres, whereas the solvent is treated as a dielectric continuum. In order to speed up the simulations, we restrict the positions of the particles to a cubic lattice, which allows precalculation of the Coulombic interactions at the beginning of the simulation. Moreover, we use multiparticle cluster moves that make the Monte Carlo sampling more efficient. The simulations are performed in the semigrand canonical ensemble, where the chemical potential of the salt is fixed. Employing our method, we study a system consisting of colloids carrying a charge of 80 elementary charges and monovalent co- and counterions. At the colloid densities of our interest, we show that lattice effects are negligible for sufficiently fine lattices. We determine the fluid-solid melting line in a packing fraction eta-inverse screening length kappa plane and compare it with the melting line of charged colloids predicted by the Yukawa potential of the Derjaguin-Landau-Verwey-Overbeek theory. We find qualitative agreement with the Yukawa results, and we do not find any effects of many-body interactions. We discuss the difficulties involved in the mapping between the primitive model and the Yukawa model at high colloid packing fractions (eta>0.2).