Differently from simulations of simple fluids, implementation of computer simulations for Coulomb systems requires use of advanced techniques to circumvent the problems associated with the long-range nature of the Coulomb force , which prevents application of the minimal-image periodic boundary condition. Instead, an infinite number of replicas of the system must be considered, alongside with the long-range Coulomb interaction between the ions in these replicas. Ewald summation method can then be used to efficiently perform the sum over the replicas. Nevertheless, a large computational cost is required to simulate dense charged systems, when compared to short-range interacting systems . The overall energy of the replicated system depends on the choice of the boundary condition for the resulting crystal of the periodic replicas (tin-foil, vacuum, slab), and care must be taken to choose the boundary condition appropriate to the physical situation at hand. Moreover, it is sometimes convenient (if not necessary) to study systems that lack charge neutrality, in which case a proper regularization of the Coulomb interaction must be perform to avoid infinite energies . The presence of acid/base groups poses additional challenges . When performing MC simulations, the possibility of a reaction between a proton and an adsorption site must be taken into account, in addition to the standard MC moves. This defines a reactive Monte Carlo (rMC) simulation, in which equilibrium is established considering MC moves that account for the protonation/deprotonation reactions, in a way that satisfies the detailed balance condition . Furthermore, in a canonical simulation, "potential of hydrogen" defined as pH = -log 10 (a H ) -where the proton activity is a H = c H c e -βµex , c H is the average concentration of free hydroniums inside the simulation at equilibrium, c = 1 M is the standard concentration, µ ex is the excess electrochemical potential of a hydronium ion, and β = 1/k B T is the inverse thermal energy -is not known a priori and must be calculated separately after the equilibrium has been achieved. In principle, a H can be obtained via the traditional Widom insertion method . In practice, however, the calculation becomes inaccurate due to the typically small number of hydronium ions inside the simulation cell for intermediate and large pH values. Moreover, insertion of a hydronium ion results in a charge non-neutral periodic system, requiring special techniques to avoid divergence of the excess chemical potential. A recently proposed alternative simulation method, designed to circumvent these difficulties, relies on the surface Widom insertion algorithm with a special technique for computing the change in the Coulomb energy when a virtual proton is inserted into the system .