Consider a polyelectrolyte or a colloidal particle with monomers that can undergo a protonation deprotonation reaction: In a canonical reactive MC simulation there are two types of movements: the bulk movements in which positions of ions are randomly changed with the acceptance probabilities given by the usual Metropolis algorithm; and reaction protonation/deprotonation moves, see Fig. . To construct a MC algorithm for the reaction moves, we first observe that the acid dissociation constant is the inverse of the two body partition function for the formations of a HA molecule. Thus, when a proton is transferred from the bulk to the surface, where it will react with a surface group A -, there are two changes that occur in the system: change of electrostatic energy ∆E and change in the chemical energy k B T ln(K a /c ⊖ ). The probabilities for the old (o) and new (n) configurations during a protonation move are proportional to: