During the evolution of the system we extract periodically two quantities: energy Eq.(\[Hamiltonian\]) and magnetization Eq.(\[magnetization\]). The kernel responsible for this job is [`calculateCUDA`]{}. It first partitions the cells into CUDA blocks. In each block we have easy access to barrier synchronization and shared memory among its threads. Each block within its cells adds the local energies and accumulates in a partial vector $(n_1, n_2,\ldots, n_q)$ the number of spins in each state. This is performed in shared memory using atomic increments to avoid race conditions. After that, those blocks’ results are added up in parallel using a butterfly-like algorithm [@Kirk10] by kernels [`sumupECUDA`]{} and [`sumupMCUDA`]{}, but none of the known optimizations [@harris07cuda] are applied, since it implies obfuscating the code for a marginal global speedup. Previous kernels end up with up to approximately a thousand partial energies and vectors of spin counters, that are finally added in the CPU.

It has to be noticed that device memory consumption in this part is linear not only in $N$, but also in $q$.

Defensive Programming Techniques and Code Availability
------------------------------------------------------

Writing scientific code that is maintainable, robust and repeatable is of utmost importance for the fields of science where computer simulation 