We present a simulation method that allows us to calculate the titration curves for systems undergoing protonation/deprotonation reactions—such as charged colloidal suspensions with acidic/basic surface groups, polyelectrolytes, polyampholytes, and proteins. The new approach allows us to simultaneously obtain titration curves both for systems in contact with salt and acid reservoir (semi-grand canonical ensemble) and for isolated suspensions (canonical ensemble). To treat the electrostatic interactions, we present a new method based on Ewald summation—which accounts for the existence of both Bethe and Donnan potentials within the simulation cell. We show that the Donnan potential dramatically affects the pH of a suspension. Counterintuitively, we find that in suspensions with a large volume fraction of nanoparticles and low ionic strength, the number of deprotonated groups can be 100% larger in an isolated system, compared to a system connected to a reservoir by a semi-permeable membrane—both systems being at *exactly* the same pH.

In many physics, chemistry, and biology applications, one needs to understand systems with active surface groups.^{1–16} Depending on the pH and salt concentration inside the solution, a surface group can be either protonated or deprotonated.^{3,17–21} The state in which a group finds itself depends on its intrinsic acid dissociation constant, the local concentration of hydronium ions, the electrostatic potential in the vicinity of the group, the steric repulsion between ions, etc. For some simple surfaces—such as metal oxide or colloidal particles with regular arrangement of surface groups, it is possible to develop a theoretical approach, which can fairly accurately predict the resulting surface charge and its dependence on salt and pH inside the solution.^{22–25} Unfortunately, no purely theoretical approach is possible for more complex systems, such as polyelectrolytes or proteins. Furthermore, even in the case of simple surfaces in solutions containing multivalent ions, one finds that theoretical approaches begin to fail.^{26,27} For such complex systems, one is forced to rely on computer simulations.^{19,28–30} The constant pH (cpH) simulation method is the most popular tool to obtain titration curves of such systems.^{31,32} There is, however, a fundamental difficulty in using such an approach indiscriminately. In cpH simulations, proteins, colloidal particles, or polyelectrolytes are confined inside a simulation box, while ions can be freely exchanged with the reservoir of acid and salt.^{33–35} The cpH simulation method is, therefore, intrinsically semi-grand canonical.^{30} When performing a cpH simulation, a proton is implicitly brought into the system from an external reservoir at a fixed pH. To keep the charge neutrality inside the simulation cell, one of the cations or protons inside the bulk of the cell is arbitrarily deleted. However, such an arbitrary deletion does not respect the detailed balance and leads to incorrect results, unless the system contains a lot of salt and is very dilute in polyelectrolyte.^{28,30} Fortunately, it is easy to correct the standard cpH algorithm by combining a protonation move with a simultaneous grand canonical insertion of an anion and a deprotonation move with a simultaneous grand canonical deletion of an anion.^{19,28–30,36,37} This restores the detailed balance of the cpH algorithm, making it internally consistent. This corrected version of the cpH algorithm will be used in the present Communication.

In a real physical system, to confine polyelectrolytes or colloidal particles within some region of space requires a semi-permeable membrane, which would separate the system from the reservoir containing acid and salt. Since the counterions are at larger concentration inside the system than they are in the reservoir, they will tend to diffuse across the membrane.^{38} The net ionic flux will end when a sufficiently large electrostatic potential difference is established between the reservoir and the system, forcing the reversal of the flow. This is known as the Donnan potential.^{39,40} In fact, one does not need to have a membrane to establish a finite Donnan potential. For example, a colloidal column in a gravitational field^{41–43} will spontaneously establish a finite Donnan potential along the vertical direction. This happens because, in general, colloidal particles have a finite buoyant mass,^{44,45} which forces them to be inhomogeneously distributed inside the suspension. Meanwhile, the ions are not affected by the gravitational field and are free to diffuse throughout all space, including the top portion of suspension, where there are no colloidal particles. This results in Donnan potential along the vertical direction of colloidal suspension. In this case, the gravitational field plays the role similar to that of a membrane, separating the region in which colloidal particles can be found. The important point that is often lost when performing constant pH simulations is that the systems (simulation cell) are always at a different electrostatic potential from the reservoir. Meanwhile, the electrochemical potentials and, consequently, the pH inside the system and the reservoir are the same.

*c*

^{⊖}= 1M is the standard concentration. On the one hand, inside the reservoir, the activity of hydronium ions is $a+=cH+\u2061exp(\beta \mu exr)$, where $cH+$ is the concentration of hydronium ions inside the reservoir and

*μ*

_{ex}is their excess chemical potential. On the other hand, inside the system, $a+=\rho H+\u2061exp[\beta \mu exs+q\beta \phi D]$, where

*q*is the proton charge, $\rho H+$ is the concentration of hydronium ions inside the system, and

*φ*

_{D}is the Donnan potential difference between the reservoir and the system. Note that the concentration of hydronium ions inside the system is different from their concentration inside the reservoir. Furthermore, while inside the reservoir, $\mu exr$ is due only to the interaction between the ions, inside the system, $\mu exs$ also includes contributions coming from electrostatic and steric interactions of hydronium ions with the colloidal particles. Nevertheless, equivalence between electrochemical potentials inside the system and the reservoir requires that pH of both be the same. (For the additional discussion of the definition of pH in inhomogeneous systems see the Conclusions of this Communication.) We see the important role that the Donnan potential plays for determining pH inside a system that is in contact with a reservoir. If one uses hydrogen and calomel (reference) electrodes to measure pH inside a system connected to a reservoir by a semi-permeable membrane—with calomel placed in the reservoir—the Donnan potential will be an integral part of the measurement. Suppose now that the system is equilibrated and then disconnected from the reservoir—so that the number of ions of each type remains fixed and no longer fluctuates. Clearly based on the ensemble equivalence of statistical mechanics, the number of protonated groups in such a

*canonical*system will remain exactly the same as in the original semi-grand canonical system. However, the pH of such a canonical system will no longer be the same as that of the semi-grand canonical system connected to the reservoir. The reason for this is that inside a canonical system, the potential drop between system and reservoir will no longer exist and the activity of hydronium ions will be $a+=\rho H+\u2061exp[\beta \mu exs]$. Substituting this into the definition of pH, we find the relation between the canonical and grand canonical pH,

Equation (1) opens a way for us to use the grand-canonical Monte Carlo (GCMC) simulation to simultaneously calculate titration curves for both canonical and grand-canonical systems. However, to do this, the usual GCMC algorithm, which relies on cation–anion pair insertion into the simulation cell, in order to preserve the charge neutrality, must be modified to allow for the charge fluctuation inside the cell, in order to calculate the Donnan potential. Note that the Donnan potential is not accessible in the normal GCMC since it cancels out for pair moves involving oppositely charged particles. In this Communication, we will present a reactive GCMC-Donnan (rGCMCD) simulation method, which will allow us to calculate the number of protonated groups and the Donnan potential in a single simulation. Combining this with Eq. (1), we are able to obtain the number of protonated group for a given pH in both canonical and grand-canonical systems. We should stress that at this time, there is no alternative method that permits us to calculate titration curves for isolated (canonical) systems. In principle, one could use canonical reactive MC simulations and then attempt to use the Widom insertion method to obtain the excess chemical potential of hydronium ions inside the simulation cell. The problem, however, is that for moderate and large pH, there might not be any hydronium ions inside the simulation cell, thus preventing us from obtaining an accurate canonical pH.

*Q*

_{t}=

*∑*

_{i}

*q*

_{i}, where the sum is over all the ions and the charged groups, is not necessarily zero at a given instant of simulation. Clearly, such an imbalance of charge in an infinitely replicated system leads to the divergence of the electrostatic energy. To avoid this, we introduce a uniform neutralizing background of charge density

*ρ*

_{b}= −

*Q*

_{t}/

*V*, where

*V*is the volume of the simulation cell. The charge density inside the system can then be written as

^{46–48}The long range can be efficiently summed in the Fourier space, while the short range is summed in the real space,

^{49–51}resulting in the electrostatic potential at position

**inside the simulation cell,**

*r***= (**

*n**n*

_{1},

*n*

_{2},

*n*

_{3}) are integers, $k=(2\pi Ln1,2\pi Ln2,2\pi Ln3)$ are the reciprocal lattice vectors, and

*κ*

_{e}is an arbitrary damping parameter. The last term is the electrostatic potential due to the uniform background, of which the Fourier transform is

Note that the electrostatic potential, Eq. (A1), is invariant to the specific value of the damping parameter *κ*_{e}. It is usual to choose *κ*_{e} sufficiently large so that simple periodic boundary conditions can be used for the short range part of the electrostatic potential.

**→ 0 in the Fourier sum of Eq. (A1) is singular and must be performed with great care. The limiting procedure is shown in Appendix A, where we obtain**

*k***=**

*M**∑*

_{i}

*q*

_{i}

*r*_{i}is the electric dipole moment inside the cell and

*ϕ*

_{B}is the modified Bethe potential of the cell with a neutralizing background,

We recognize the ** M** dependent term in Eq. (A13) as arising from the uniform polarization of the macroscopic crystal composed of spherically replicated simulation cells. From continuum electrostatics, such uniform polarization is equivalent to the surface charge density

**·**

*M***/**

*n**V*, where

**is the unit normal to the boundary of the macroscopic spherical crystal, resulting in an electrostatic potential $4\pi 3\u03f5wVr\u22c5M$ in the interior of the crystal; see Fig. 1.**

*n*The nature of *ϕ*_{B} is more subtle. In general, a simulation cell containing colloidal particles or polyelectrolytes will have a non-zero trace of the second moment charge density tensor. Spherical replication of such cells will lead to a crystal with a dipole surface layer,^{52} across which there will be a potential drop precisely given by Eq. (6). In Appendix B, we also provide an alternative real space derivation of Eq. (6).

Since ions can diffuse across the membrane, while colloidal particles are confined within the system, a Donnan potential exists between the reservoir and the system. While the Bethe potential can be thought of as the average of the local electrostatic potential inside the periodically replicated crystal, the Donnan potential is inherently an interfacial effect that arises from ionic diffusion between the system and the reservoir. The total electrostatic potential drop between the system and the reservoir is shown in Fig. 2.

*κ*

_{e}, the limit in the last term can be calculated using

*κ*

_{e}→ ∞, which kills off the real space sum in Eq. (A13). Taking the limit lim

_{k→0}, the last term of Eq. (7) then evaluates to

*Q*

_{t}

*ϕ*

_{B}/2 and the total electrostatic energy of a system with ions + sites of total charge

*Q*

_{t}and a neutralizing background can be written as

*i*=

*j*is excluded from the summation when

**n**= 0. Note that the Bethe potential canceled out and does not contribute to the total electrostatic energy. Furthermore, since the system is implicitly connected to a reservoir, it is not possible to use the tinfoil boundary condition,

^{53}—which would eliminate

**= 0 term from the Ewald summation, isolating the system from the reservoir.**

*k**q*

_{i}to enter the simulation cell, it must cross the macroscopic boundary of the crystal. This will result in electrostatic energy change

*q*

_{i}(

*φ*

_{D}+

*ϕ*

_{B}). The grand canonical acceptance probabilities for insertion/deletion moves can then be written as $min1,\varphi add\u2215rem$, where

*ϕ*

_{add∕rem}is

*E*

_{ele}is the change in electrostatic energy upon insertion/deletion of ion of type

*i*and pX

_{i}≡ −log

_{10}

*X*

_{i}/

*c*

^{⊖}, in which

*X*

_{i}is the activity of ion

*i*inside the reservoir. To relate pX

_{i}to the concentration of salt, one can perform a separate GCMC simulation for the reservoir. In practice, however, we find that for monovalent ions, the mean spherical approximation (MSA) combined with the Carnahan–Starling (CS) expression for the excluded volume

^{54,55}provides an excellent approximation for the relation between pX

_{i}, pH, and the concentrations of salt and acid in the reservoir.

^{−}is the deprotonated site. We recognize the acid dissociation constant to be the inverse of the internal partition function of HA molecule. The free energy change due to the removal of hydronium from the reservoir and its reaction with an

*isolated*A

^{−}group is then

*F*

_{d}= −Δ

*F*

_{p}. The acceptance probabilities for protonation and deprotonation moves are then $min1,\varphi p/d$, where

_{a}of a titration site as pK

_{a}= − log

_{10}[K

_{a}/

*c*

^{⊖}]. During the protonation/deprotonation moves, a titration site is randomly chosen and its state is changed with an acceptance probability given by Eq. (13).

The simulation starts with an initial guess for the value of *φ*_{D}. After each MC step, the Donnan potential is automatically updated so as to drive the system toward a charge neutral state, $\phi Dnew\u2254\phi Dold+\alpha Qt$, where *α* > 0 is an arbitrary small parameter that controls convergence. The simulation stops when ⟨*Q*_{t}⟩ ≈ 0, to the desired accuracy. In practice, the system very quickly converges to a state in which *Q*_{t} has small oscillations about zero. To make sure that the system is fully equilibrated, we monitor the energy. When ⟨*E*⟩ becomes constant, we know that the system has reached equilibrium. Usually, we use 5 × 10^{6} particle moves for equilibration. After this, we collect samples, which include colloidal charge, Donnan potential, and ion distribution. The samples are collected with intervals of ∼3000 particle moves, to make sure that they are fully uncorrelated. We use 20 000 such samples to calculate the averages.

As an example of the algorithm developed in the present Communication, we will use it to study a colloidal suspension of finite volume fraction. The colloidal particles have radius 60 Å and contain *Z* = 600 active surface groups of intrinsic pK_{a} = 5.4 uniformly distributed over the surface. All ions have radius 2 Å, and the solvent is modeled as a dielectric continuum of Bjerrum length, *λ*_{B} = *q*^{2}/*ϵk*_{B}*T* = 7.2 Å. The reservoir contains acid of an arbitrary pH and 1:1 salt of concentration 1 mM. The electrostatic energy is calculated using Eq. (A11), with the damping parameter set to *κ*_{e} = 5/*L*, where *L* is the size of the cubic simulation box. Such a value of *κ*_{e} is sufficiently large that simple periodic boundary conditions can be used to evaluate the short range (sum over erfc) contribution to the electrostatic energy in Eq. (A11).

*L*= 200 Å. There is, however, no conceptual difficulty in using an arbitrary number of colloidal particles inside a simulation box, except for the CPU time constraint. Finally, to test the new rGCMCD algorithm, we perform a “standard” cpH in which a protonation move is combined with a grand canonical insertion of anion and a deprotonation move with a grand canonical deletion of an anion.

^{28,30,56}The acceptance probabilities for these pair moves are, respectively,

In Fig. 3, we present the titration curves obtained using our rGCMCD simulations for both canonical (isolated) and semi-grand canonical suspensions.

As can be seen from the figure, there is a very significant difference between titration curves of an isolated colloidal suspension and those of suspension connected to a salt and acid reservoir. For example, at pH = 7.5, 28% of surface groups of semi-grand canonical system are deprotonated, while for an isolated canonical system, this jumps to 57%. When the concentration of salt increases, the difference between canonical and semi-grand canonical systems diminishes; see Fig. 4. This justifies the use of cpH methods to study biologically relevant systems that contain physiological concentrations of salt of *c*_{s} ≈ 150 mM, which are usually isolated from any external reservoir.

We have presented a new simulation method that enables us to simultaneously obtain titration curves for both isolated systems (canonical ensemble) and systems connected to an external reservoir by a semi-permeable membrane. Counterintuitively, we find that at low ionic strength, the number of deprotonated groups can be 100% larger in an isolated system, compared to a system connected to a reservoir by a semi-permeable membrane—both systems at exactly the same pH! The difference can be even greater for more concentrated (large colloidal volume fraction) suspensions in deionized solutions. As the concentration of salt increases, the difference between the ensembles becomes less important, justifying the use of cpH methods for studying biologically relevant systems that are usually closed and contain large concentrations of salt. It is important to stress that at the moment there is no alternative method available to study canonical systems at moderate and high pH. The usual canonical reactive simulations fail under such conditions since the simulation cell contains none or very few hydronium ions, preventing one from accurately calculating the pH inside the system.

In our definition of pH for heterogeneous systems, we have followed the Gibbs–Guggenheim principle,^{57} which prohibits the separation of the electrochemical potential into chemical and electrostatic potential contributions for the definition of measurable quantities. Meanwhile, in the electrochemistry literature, the Gibbs–Guggenheim principle is often disregarded and pH is defined in terms of only the chemical part of the total electrochemical potential. From the practical point of view, pH is determined by measuring the electromotive force (EMF) between a glass or hydrogen electrode and a saturated calomel (reference) electrode. In our definition of pH, the reference electrode is always kept inside the reservoir, while the hydrogen electrode is moved between the reservoir and the system. With this setup, we will obtain the same EMF independent of the position of the hydrogen electrode, indicating the same pH inside the system and in the reservoir.^{58} Meanwhile, if both electrodes are moved together—changing the reference value of the electrostatic potential—we will obtain different pH values inside the system and in the reservoir. With such a definition, the system’s pH will be different from that of the reservoir but will be the same as that of an isolated canonical system. Whichever definition of pH is chosen, one needs to know the Donnan potential in order to calculate the titration isotherms of isolated (canonical) systems using semi-grand canonical simulations. The “standard” cpH simulation methods do not allow us to calculate titration isotherms for isolated systems. One must either use canonical reactive MC—with Widom particle insertion, which is very inaccurate even for intermediate to high values—or the approach presented in this Communication, which is accurate for any pH.

This work was partially supported by the CNPq, the CAPES, and the National Institute of Science and Technology Complex Fluids (INCT-FCx).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Yan Levin**: Conceptualization (lead); Supervision (lead); Writing – original draft (equal); Writing – review & editing (equal). **Amin Bakhshandeh**: Data curation (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

### APPENDIX A: DETAILS FOR THE LIMITING PROCEDURE FOR SINGULAR TERM OF ELECTROSTATIC POTENTIAL

*V*=

*L*

^{3}is the volume of the cubic simulation cell,

*Q*

_{t}=

*∑*

_{i}q

_{i}, and

**, except for**

*k***= 0, so that we can define the singular part of the background potential**

*k**ϕ*

_{b,s}as $\varphi \u0303b(k)\u2261\varphi \u0303b,s\delta k0$, where we have defined the Kronecker delta for the zero mode,

*δ*

_{k0}. As

**→ 0, we can expand the singular part around**

*k***= 0,**

*k**k*are

**→ 0. We note that the limit**

*k***→ 0 corresponds to a large distance behavior of the electrostatic potential. To clearly see how the divergence scales with the size of the macroscopic crystal,**

*k**L*

_{cr}, produced by the periodic replication of simulation cell [see Fig. (1)], we introduce a “crystal” delta function as follows:

*L*

_{cr}→ ∞,

*δ*

_{cr}(

**) becomes the Dirac delta function. The crystal**

*k**δ*

_{cr}(

**) regularizes the divergence and allows us to study the**

*k***→ 0 limit.**

*k*^{50}In particular,

*δ*

_{ij}is the Kronecker delta. The last equality follows from the direct use of the representation of

*δ*

_{cr}(

**), Eq. (A5), followed by the**

*k**L*

_{cr}→ ∞ limit,

^{50}or simply from the observation of symmetry. Similarly,

*δ*

_{cr}(

**) is an even function of all**

*k**k*

_{i}. Using Eq. (A7), we see that the second term of Eq. (A4) vanishes and the last term becomes

**= 0. The singular terms are**

*k***= 0, with the help of Eq. (A9), results in**

*k**k*

^{2}terms cancel out and we obtain the electrostatic potential inside the simulation cell,

### APPENDIX B: REAL SPACE DERIVATION OF THE BETHE POTENTIAL INSIDE A CELL WITH A NEUTRALIZING BACKGROUND

*κ*

_{e}. In the limit

*κ*

_{e}→ ∞, the sum over erfc vanishes, and the electrostatic potential inside the cell can be written in terms of only the Fourier components. It is clear then that

*V*

_{cr}= 4

*πR*

^{3}/3,

**is the potential at position**

*r***′ inside a uniformly charged sphere of unit charge density, which can be easily obtained from the solution of Poisson equation,**

*r*## REFERENCES

*Intermolecular and Surface Forces*

*A New Monte Carlo Method for the Titration of Molecules and Minerals*

*Understanding Molecular Simulation: From Algorithms to Applications*

*Computer Simulation of Liquids*