where φ B is the Bethe potential of the periodically replicated system, φ B = -2π 3εwV i q i r r r i 2 . The surface Widom method involves introduction of a virtual proton into the system, which then reacts with one of the surface sites. If the randomly chosen site is already protonated, this results in ∆E = ∞, otherwise the energy of the virtual protonate state is calculated neutralizing the charge of the selected site, -q → 0. Since the virtual proton is periodically replicated when it is introduced into the system, we must also include the corresponding neutralizing background. The resulting change in the electrostatic energy, ∆E, can be calculated using Eq. ( ). The intrinsic pK a of the reaction is defined as pK a = -log 10 (K a /c ), where c = 1 M is the standard concentration. The subscript 0 on the brackets in Eqn. (6) indicates that the evolution of the system in between the virtual proton insertions is performed using the Hamiltonian of the unperturbed (without virtual proton) system.