where, N H is the number of free hydronium ions inside the simulation cell at a given moment and N A is the Avogadro number. Note that for canonical (closed) systems, the acceptance probabilities do not explicitly depend on pH. The pH of the system, after the equilibrium is established, has to be calculated using a separate procedure that require determination of the chemical potential of hydronium ions inside the simulation cell. The conventional approach of employing Widom's particle insertion method encounters practical challenges for moderate and high pH, when cell's interior has few or maybe even no free hydroniums at all, resulting in a very inaccurate reading of the chemical potential of hydronium ions. To overcome this difficulty, we have recently introduced a surface Widom insertion algorithm, which allows us to easily calculate the pH inside the system after equilibrium has been established. One difficulty when working with Coulomb systems is that we can not cutoff the long range electrostatic interactions between particles and use simple periodic boundary conditions. Instead the system has to be infinitely replicated. To efficiently sum the replicas, we use Ewald summation method, which effectively creates a spherical macroscopic crystal comprised of periodically replicated microscopic simulation cells. Typically a simulation cell will have a net electric dipole moment M M M = i q i r r r i , and a finite second moment tensor of the charge density. Electrostatically, the field produced by a dielectric sphere of uniform polarizability M M M is equivalent to that of a sphere with non-uniform surface charge density M M M • n n n/V , where n n n represents the unit normal to the surface of the macroscopic sphere. The existence of an effective surface charge induces an electric field within the crystal's interior.