on GPU based Statistical Mechanics simulations of lattice spin systems. In particular, we study the metastability problem in the ferromagnetic $q$-state Potts model [@Wu1982] in two dimensions when $q > 4$. While this phenomenon is clearly observed in finite size systems, its persistence in the thermodynamics limit is still an unsolved problem and subject of debate [@Bi1981; @MeMo2000; @PeIbLo2008; @BaBeDu2008; @LoFeGrCa2009]. In an earlier work, Binder proposed a numerical criterion to determine whether metastability remains in the thermodynamic limit or not, based on the scaling properties of the average energy in the vicinity of the transition temperature [@Bi1981]. However, the narrow range of temperature values of the metastable region requires high precision calculations for the criterion to work. Hence, to reduce finite size bias and statistical errors down to an appropriated level, large enough system sizes are needed. The computation capabilities required to carry out such calculations in a reasonable time were unavailable until recently.

We developed an optimized algorithm to perform Monte Carlo numerical simulations of the $q$-state Potts model on GPU cards. This algorithm allowed us to simulate systems up to $N=32768\times32768 \sim 1.073 \times 10^9$ spins with a lower bound time of 0.147ns per spin flip using 