A new method for reactive constant pH simulations

Levin, Yan; Bakhshandeh, Amin

doi:10.1063/5.0166840

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.³⁰ 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.³⁸ 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.

Recall that pH of a solution is defined in terms of activity of hydronium ions, pH =

- \log_{10} [a_{+} / c^{⊖}]

⁠, where c^⊖ = 1M is the standard concentration. On the one hand, inside the reservoir, the activity of hydronium ions is

a_{+} = c_{H^{+}} \exp (β μ_{e x}^{r})

⁠, where

c_{H^{+}}

is the concentration of hydronium ions inside the reservoir and μ_ex is their excess chemical potential. On the other hand, inside the system,

a_{+} = ρ_{H^{+}} \exp [β μ_{e x}^{s} + q β φ_{D}]

⁠, where q is the proton charge,

ρ_{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,

μ_{e x}^{r}

is due only to the interaction between the ions, inside the system,

μ_{e x}^{s}

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_{+} = ρ_{H^{+}} \exp [β μ_{e x}^{s}]

⁠. Substituting this into the definition of pH, we find the relation between the canonical and grand canonical pH,

{pH}_{c} = {pH}_{g c} + \frac{β φ_{D}}{\ln (10)} .

(1)

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.

The difficulty in performing simulations of Coulomb systems stems from the long range electrostatic interaction between the particles that cannot be cut off at the cell boundary. Instead, one is forced to periodically replicate the whole simulation cell so that a given ion interacts not only with the ions inside the cell but also with their infinite replicas. For Donnan simulations presented in this Communication, there is an additional difficulty because the simulation cell is no longer charge neutral. The total charge inside the cell, Q_t = ∑_iq_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

ρ_{q} (r) = \sum_{i} q_{i} δ (r - r_{i}) - \frac{Q_{t}}{V} .

(2)

Using the usual approach, the electrostatic potential can be split into long and short range contribution.^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 r inside the simulation cell,

\begin{matrix} φ (r) & = & \sum_{k = 0}^{\infty} \sum_{j = 1}^{N} \frac{4 π q^{j}}{ϵ_{w} V | k |^{2}} \exp [- \frac{| k |^{2}}{4 κ_{e}^{2}} + i k \cdot (r - r^{j})] \\ + \sum_{j = 1}^{N} \sum_{n} q^{j} \frac{erfc (κ_{e} | r - r^{j} - L n |)}{ϵ_{w} | r - r^{j} - L n |} \\ + \frac{1}{V} \sum_{k = 0}^{\infty} {\tilde{φ}}_{b} (k) \exp [i k \cdot r], \end{matrix}

(3)

where n = (n₁, n₂, n₃) are integers,

k = (\frac{2 π}{L} n_{1}, \frac{2 π}{L} n_{2}, \frac{2 π}{L} n_{3})

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

{\tilde{φ}}_{b} (k) = - \frac{4 π Q_{t}}{ϵ_{w} V} \frac{\int_{V} e^{- i k . r} d^{3} r}{k^{2}} .

(4)

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.

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

\begin{matrix} φ (r) & = & \sum_{k \neq 0}^{\infty} \sum_{j = 1}^{N} \frac{4 π q^{j}}{ϵ_{w} V | k |^{2}} \exp [- \frac{| k |^{2}}{4 κ_{e}^{2}} + i k \cdot (r - r^{j})] \\ + \sum_{j = 1}^{N} \sum_{n} q^{j} \frac{erfc (κ_{e} | r - r^{j} - L n |)}{ϵ_{w} | r - r^{j} - L n |} \\ - \frac{π Q_{t}}{ϵ_{w} V κ_{e}^{2}} + \frac{4 π}{3 ϵ_{w} V} r \cdot M + ϕ_{B}, \end{matrix}

(5)

where M = ∑_iq_ir_i is the electric dipole moment inside the cell and ϕ_B is the modified Bethe potential of the cell with a neutralizing background,

ϕ_{B} = - \frac{2 π}{3 ϵ_{w} V} \sum_{i} q_{i} {r_{i}}^{2} + \frac{π Q_{t}}{6 ϵ_{w} L} .

(6)

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 n is the unit normal to the boundary of the macroscopic spherical crystal, resulting in an electrostatic potential $\frac{4 π}{3 ϵ_{w} V} r \cdot M$ in the interior of the crystal; see Fig. 1.

FIG. 1.

View large Download slide

Replicas of the simulation cell. The large sphere is a colloidal particle and the small spheres are ions; the arrows show the instantaneous electric dipole moment M, which, in general, exists inside the cell. Furthermore, the trace of the instantaneous second moment charge distribution tensor of a locally strongly inhomogeneous system of colloidal particles and ions will not be zero. A macroscopic spherical crystal produced by a periodic replication of the simulation cell will then have both surface charge density M · n/V and a dipolar layer, which results in a potential drop across the crystal–reservoir interface, given by the modified Bethe potential, ϕ_B.

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,⁵² 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.

FIG. 2.

View large Download slide

Periodically replicated system separated from the reservoir by a semi-permeable membrane: the electrostatic potential difference between the two is ϕ_B + φ_D.

The electrostatic energy can be calculated using

\begin{matrix} E & = & \frac{1}{2} \int ρ_{q} (r) φ (r) d^{3} r \\ = & \frac{1}{2} \sum_{i} q_{i} \lim_{r \to r_{i}} [φ (r) - \frac{q_{i}}{ϵ_{w} | r - r_{i} |}] \\ - \frac{Q_{t}}{2 V} \lim_{k \to 0} \tilde{φ} (k) . \end{matrix}

(7)

The subtraction inside the square brackets eliminates the self-energy contribution to the total electrostatic energy of all ions and charged groups. The last term is due to the background charge. Noting that the electrostatic energy is invariant with the damping parameter κ_e, the limit in the last term can be calculated using κ_e → ∞, which kills off the real space sum in Eq. . Taking the limit lim_k→0, the last term of Eq. 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

\begin{matrix} E & = & \frac{1}{2} {\sum_{i j}}^{'} \sum_{n} \frac{q_{i} q_{j} erfc (κ_{e} | r_{i} - r_{j} - L n |)}{ϵ_{w} | r_{i} - r_{j} - L n |} \\ + \sum_{k \neq 0} \frac{2 π \exp (- k^{2} / 4 κ_{e})}{ϵ_{w} V k^{2}} (A {(k)}^{2} + B {(k)}^{2}) \\ - \sum_{i} \frac{q_{i}^{2} κ_{e}}{ϵ_{w} \sqrt{π}} - \frac{π Q_{t}^{2}}{2 ϵ_{w} V κ_{e}^{2}} + \frac{2 π}{3 ϵ_{w} V} M^{2}, \end{matrix}

(8)

where

\begin{gathered} A (k) = \sum_{i} q_{i} \cos (k \cdot r_{i}), \\ B (k) = \sum_{i} q_{i} \sin (k \cdot r_{i}), \end{gathered}

(9)

and the prime on the sum indicates that 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,⁵³—which would eliminate k = 0 term from the Ewald summation, isolating the system from the reservoir.

The simulation involves individual insertion/deletion moves, as well as protonation and deprotonation moves. For an ion of charge 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

\min (1, ϕ_{add∕rem})

⁠, where ϕ_add∕rem is

\begin{gathered} ϕ_{add} = \frac{c^{⊖} V 1 0^{- {pX}_{i}}}{N_{i} + 1} e^{- β (Δ E_{ele} + q_{i} [φ_{D} + ϕ_{B}])}, \\ ϕ_{rem} = \frac{N_{i} 1 0^{{pX}_{i}}}{V c^{⊖}} e^{- β (Δ E_{ele} - q_{i} [φ_{D} + ϕ_{B}])}, \end{gathered}

(10)

where ΔE_ele is the change in electrostatic energy upon insertion/deletion of ion of type i and pX_i ≡ −log₁₀ 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.

The protonation/deprotonation moves involve a reaction,

\begin{gathered} HA ⇌ H_{3} O^{+} + A^{-}, \\ K_{a} = \frac{a_{A^{-}} a_{H^{+}}}{a_{HA}}, \end{gathered}

(11)

where A⁻ 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_{p} = \ln (\frac{K_{a}}{c^{⊖}}) - μ_{H^{+}},

(12)

where

μ_{H^{+}} = \ln (c_{H^{+}} / c^{⊖}) + β μ_{e x}

is the total chemical potential of a hydronium ion inside the reservoir. For the deprotonation reaction, the free energy change is reversed so that ΔF_d = −ΔF_p. The acceptance probabilities for protonation and deprotonation moves are then

\min (1, ϕ_{p / d})

⁠, where

\begin{aligned} ϕ_{p} & = e^{- β (Δ E_{ele} + Δ F_{p} + q [φ_{D} + ϕ_{B}])} \\ = 1 0^{{pK}_{a} - {pH}_{a}} e^{- β (Δ E_{ele} + q [φ_{D} + ϕ_{B}])}, \\ ϕ_{d} & = e^{- β (Δ E_{ele} + Δ F_{d} - q [φ_{D} + ϕ_{B}])} \\ = 1 0^{pH - {pK}_{a}} e^{- β (Δ E_{ele} - q [φ_{D} + ϕ_{B}])} . \end{aligned}

(13)

We have defined the intrinsic pK_a of a titration site as pK_a = − log₁₀[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. .

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, $φ_{D}^{new} ≔ φ_{D}^{old} + α Q_{t}$ ⁠, 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⁶ 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²/ϵk_BT = 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).

For simplicity, in the present Communication, we will use only one colloidal particle inside a cubic simulation box of 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,

\begin{gathered} P_{p} = \min [1, \frac{c^{⊖} V 1 0^{{pK}_{a} - pH - pCl}}{N_{Cl} + 1} e^{- β Δ E}], \\ P_{d} = \min [1, \frac{N_{Cl} 1 0^{pH - {pK}_{a} + pCl}}{c^{⊖} V} e^{- β Δ E}] . \end{gathered}

(14)

Since both Bethe and Donnan potentials cancel out in the pair moves, cpH simulation does not permit us to calculate the canonical titration curve. Meanwhile, the grand-canonical curve obtained using the “standard” cpH must agree precisely with the one obtained using our new rGCMCD simulation method. This is exactly what is found, validating the rGCMCD method.

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

FIG. 3.

$Titration curves for an isolated (canonical, solid curve) and semi-grand canonical system (dashed curve) in contact with a salt and acid reservoir, calculated using rGCMCD simulations. In the semi-grand canonical case, the reservoir contains 1 mM of salt; for the canonical case, the concentration in reservoir is fine-tuned so that there is 1 mM of salt inside the system for all pH. This, however, has an imperceptible effect on the canonical titration curve for low salt concentrations. The dots are the results of the “standard” cpH with pair insertions, which agree perfectly with the new rGCMCD method. The colloidal suspension has a volume fraction of 11.3%. The saturated colloidal charge is −212 mCm−2. The pKa of acidic groups is 5.4.$

View large Download slide

Titration curves for an isolated (canonical, solid curve) and semi-grand canonical system (dashed curve) in contact with a salt and acid reservoir, calculated using rGCMCD simulations. In the semi-grand canonical case, the reservoir contains 1 mM of salt; for the canonical case, the concentration in reservoir is fine-tuned so that there is 1 mM of salt inside the system for all pH. This, however, has an imperceptible effect on the canonical titration curve for low salt concentrations. The dots are the results of the “standard” cpH with pair insertions, which agree perfectly with the new rGCMCD method. The colloidal suspension has a volume fraction of 11.3%. The saturated colloidal charge is −212 mCm⁻². The pK_a of acidic groups is 5.4.

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.

FIG. 4.

$Titration curves for an isolated (canonical, solid curve) and semi-grand canonical system (dashed curve) in contact with a salt and acid reservoir, calculated using rGCMCD simulations. In the semi-grand canonical case, the reservoir contains 10 mM of salt; for the canonical case, the concentration in reservoir is fine-tuned so that there is 10 mM of salt inside the system for all pH. The dots are the results of cpH pair insertion simulation. The colloidal suspension has a volume fraction of 11.3%. The saturated colloidal surface charge is −212 mCm−2. The pKa of acidic groups is 5.4.$

View large Download slide

Titration curves for an isolated (canonical, solid curve) and semi-grand canonical system (dashed curve) in contact with a salt and acid reservoir, calculated using rGCMCD simulations. In the semi-grand canonical case, the reservoir contains 10 mM of salt; for the canonical case, the concentration in reservoir is fine-tuned so that there is 10 mM of salt inside the system for all pH. The dots are the results of cpH pair insertion simulation. The colloidal suspension has a volume fraction of 11.3%. The saturated colloidal surface charge is −212 mCm⁻². The pK_a of acidic groups is 5.4.

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,⁵⁷ 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.⁵⁸ 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

The electrostatic potential inside the simulation cell is

\begin{matrix} ϕ (r) & = & \sum_{k = 0}^{\infty} \sum_{j = 1}^{N} \frac{4 π q_{j}}{ϵ_{w} V | k |^{2}} \exp [- \frac{| k |^{2}}{4 κ_{e}^{2}} + i k \cdot (r - r_{j})] \\ + \sum_{j = 1}^{N} \sum_{n} q_{j} \frac{erfc (κ_{e} | r - r_{j} - L n |)}{ϵ_{w} | r - r_{j} - L n |} \\ + \frac{1}{V} \sum_{k = 0}^{\infty} {\tilde{ϕ}}_{b} (k) \exp [i k \cdot r], \end{matrix}

(A1)

where V = L³ is the volume of the cubic simulation cell, Q_t = ∑_iq_i, and

{\tilde{ϕ}}_{b} (k) = - \frac{4 π Q_{t}}{ϵ_{w} V} \frac{\int_{V} e^{- i k . r} d^{3} r}{k^{2}}

(A2)

is the Fourier transform of the background potential. Note that the integral in Eq. vanishes for all k, except for k = 0, so that we can define the singular part of the background potential ϕ_b,s as

{\tilde{ϕ}}_{b} (k) \equiv {\tilde{ϕ}}_{b, s} δ_{k 0}

⁠, where we have defined the Kronecker delta for the zero mode, δ_k0. As k → 0, we can expand the singular part around k = 0,

{\tilde{ϕ}}_{b, s} = - \frac{4 π Q_{t}}{ϵ_{w} V} \frac{\int_{V} (1 - i k \cdot r - \frac{1}{2} {(k \cdot r)}^{2} + \dots) d^{3} r}{k^{2}} .

(A3)

The leading order terms in k are

{\tilde{ϕ}}_{b, s} = - \frac{4 π Q_{t}}{ϵ_{w} k^{2}} + \frac{4 π Q_{t}}{ϵ_{w} k^{2} V} \int_{V} i k . r d^{3} r + \frac{2 π Q_{t}}{ϵ_{w} V k^{2}} \int_{V} {(k . r)}^{2} d^{3} r,

(A4)

with the higher order terms vanishing in the limit k → 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, L_cr, produced by the periodic replication of simulation cell [see Fig. (1)], we introduce a “crystal” delta function as follows:

\begin{matrix} δ_{c r} (k) & = & \frac{1}{{(2 π)}^{3}} \prod_{i = 1}^{3} \int_{- L_{c r} / 2}^{L_{c r} / 2} e^{i k_{i} p_{i}} d p_{i} \\ = & \frac{1}{π^{3}} \prod_{i = 1}^{3} \frac{\sin (k_{i} L_{c r} / 2)}{k_{i}} . \end{matrix}

(A5)

In the limit L_cr → ∞, δ_cr(k) becomes the Dirac delta function. The crystal δ_cr(k) regularizes the divergence and allows us to study the k → 0 limit.⁵⁰ In particular,

\lim_{k \to 0} \frac{k_{i} k_{j}}{k^{2}} = \int δ_{c r} (k) \frac{k_{i} k_{j}}{k^{2}} d^{3} k = \frac{1}{3} δ_{i j},

(A6)

where δ_ij is the Kronecker delta. The last equality follows from the direct use of the representation of δ_cr(k), Eq. , followed by the L_cr → ∞ limit,⁵⁰ or simply from the observation of symmetry. Similarly,

\lim_{k \to 0} \frac{k_{i}}{k^{2}} = \int δ_{c r} (k) \frac{k_{i}}{k^{2}} d^{3} k = 0,

(A7)

since δ_cr(k) is an even function of all k_i. Using Eq. , we see that the second term of Eq. vanishes and the last term becomes

\frac{2 π Q}{3 ϵ_{w} V} \int_{V} (x^{2} + y^{2} + z^{2}) d x d y d z = \frac{π Q L^{2}}{6 ϵ_{w}}

(A8)

so that

{\tilde{ϕ}}_{b, s} = - \frac{4 π Q_{t}}{ϵ_{w} k^{2}} + \frac{π Q_{t} L^{2}}{6 ϵ_{w}} .

(A9)

We now expand the first term of Eq. around k = 0. The singular terms are

\begin{matrix} \frac{4 π}{V ϵ_{w} k^{2}} \sum_{j = 1} q_{j} - \frac{π}{ϵ_{w} V κ_{e}^{2}} \sum_{j = 1} q_{j} + \frac{4 π}{V ϵ_{w}} \sum_{j = 1} q_{j} \frac{i k \cdot (r - r_{j})}{| k |^{2}} \\ - \frac{2 π}{V ϵ_{w}} \sum_{j = 1} q^{j} \frac{{[k \cdot (r - r_{j})]}^{2}}{| k |^{2}}, \end{matrix}

(A10)

which, using Eqs. and , simplify to

\frac{4 π Q_{t}}{ϵ_{w} k^{2} V} - \frac{π Q_{t}}{ϵ_{w} V κ_{e}^{2}} - \frac{2 π}{3 V ϵ_{w}} \sum_{j = 1} q^{j} {(r - r_{j})}^{2} .

(A11)

Meanwhile, the expansion of the last term of Eq. around k = 0, with the help of Eq. , results in

- \frac{4 π Q_{t}}{ϵ_{w} k^{2} V} + \frac{π Q_{t}}{6 ϵ_{w} L} + \frac{2 π Q_{t} r^{2}}{3 V ϵ_{w}} .

(A12)

Combining everything, we see that the diverging 1/k² terms cancel out and we obtain the electrostatic potential inside the simulation cell,

\begin{matrix} φ (r) & = & \sum_{k \neq 0}^{\infty} \sum_{j = 1}^{N} \frac{4 π q_{j}}{ϵ_{w} V | k |^{2}} \exp [- \frac{| k |^{2}}{4 κ_{e}^{2}} + i k \cdot (r - r_{j})] \\ + \sum_{j = 1}^{N} \sum_{n} q_{j} \frac{erfc (κ_{e} | r - r_{j} - L n |)}{ϵ_{w} | r - r_{j} - L n |} \\ - \frac{π Q_{t}}{ϵ_{w} V κ_{e}^{2}} + \frac{4 π}{3 ϵ_{w} V} r \cdot M - \frac{2 π}{3 ϵ_{w} V} \sum_{i} q_{i} r_{i}^{2} + \frac{π Q_{t}}{6 ϵ_{w} L} . \end{matrix}

(A13)

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

The electrostatic potential, Eq. , is invariant with respect to κ_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

\frac{1}{V} \int_{V} φ (r) d^{3} r = ϕ_{B} .

(B1)

Following Bethe, we can also calculate this as the mean potential over the whole spherical crystal of radius V_cr = 4πR³/3,

ϕ_{B} = \frac{1}{ϵ_{w} V_{c r}} \int_{V_{c r}} \int_{V_{c r}} \frac{ρ_{q} (r^{'})}{| r^{'} - r |} d^{3} r^{'} d^{3} r .

(B2)

Interchanging the orders of integration, the integral over r 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,

\int_{V_{c r}} \frac{1}{| r^{'} - r |} d^{3} r = \frac{2 π}{3} (3 R^{2} - r^{' 2}) .

(B3)

The average potential inside the simulation cell is then

ϕ_{B} = \frac{2 π}{3 V} \int_{V} (3 R^{2} - r^{' 2}) ρ_{q} (r^{'}) d^{3} r^{'} = - \frac{2 π}{3 V} \int_{V} r^{' 2} ρ_{q} (r^{'}) d^{3} r^{'},

(B4)

where in the last equality, we have used the total ions + sites + background charge neutrality inside the cell. Finally, substituting

ρ_{q} (r) = \sum_{i} q_{i} δ (r - r_{i}) - \frac{Q_{t}}{V}

into the above-mentioned equation, we obtain

ϕ_{B} = - \frac{2 π}{3 ϵ_{w} V} \sum_{i} q_{i} r_{i}^{2} + \frac{π Q_{t}}{6 ϵ_{w} L} .

(B5)

REFERENCES

1.

Y.

Avni

,

R.

Podgornik

, and

D.

Andelman

, “

Critical behavior of charge-regulated macro-ions

,”

J. Chem. Phys.

153

,

024901

(

2020

).

https://doi.org/10.1063/5.0011623

Google Scholar

Crossref

PubMed

2.

D.

Prusty

,

R. J.

Nap

,

I.

Szleifer

, and

M.

Olvera de la Cruz

, “

Charge regulation mechanism in end-tethered weak polyampholytes

,”

Soft Matter

16

,

8832

–

8847

(

2020

).

https://doi.org/10.1039/d0sm01323d

Google Scholar

Crossref

PubMed

3.

R.

Podgornik

, “

General theory of charge regulation and surface differential capacitance

,”

J. Chem. Phys.

149

,

104701

(

2018

).

https://doi.org/10.1063/1.5045237

Google Scholar

Crossref

PubMed

4.

T.

Nishio

, “

Monte Carlo studies on potentiometric titration of (carboxymethyl)cellulose

,”

Biophys. Chem.

57

,

261

–

267

(

1996

).

https://doi.org/10.1016/0301-4622(95)00057-3

Google Scholar

Crossref

PubMed

5.

Y.

Avni

,

D.

Andelman

, and

R.

Podgornik

, “

Charge regulation with fixed and mobile charged macromolecules

,”

Curr. Opin. Electrochem.

13

,

70

–

77

(

2019

).

https://doi.org/10.1016/j.coelec.2018.10.014

Google Scholar

Crossref

6.

R.

Lunkad

,

A.

Murmiliuk

,

Z.

Tošner

,

M.

Štěpánek

, and

P.

Košovan

, “

Role of pKA in charge regulation and conformation of various peptide sequences

,”

Polymers

13

,

214

(

2021

).

https://doi.org/10.3390/polym13020214

Google Scholar

Crossref

PubMed

7.

F. L. B.

da Silva

and

B.

Jönsson

, “

Polyelectrolyte–protein complexation driven by charge regulation

,”

Soft Matter

5

,

2862

–

2868

(

2009

).

https://doi.org/10.1039/b902039j

Google Scholar

Crossref

8.

T.

Curk

and

E.

Luijten

, “

Charge regulation effects in nanoparticle self-assembly

,”

Phys. Rev. Lett.

126

,

138003

(

2021

).

https://doi.org/10.1103/physrevlett.126.138003

Google Scholar

Crossref

PubMed

9.

M.

Lund

,

T.

Åkesson

, and

B.

Jönsson

, “

Enhanced protein adsorption due to charge regulation

,”

Langmuir

21

,

8385

–

8388

(

2005

).

https://doi.org/10.1021/la050607z

Google Scholar

Crossref

PubMed

10.

R. J.

Nap

,

A. L.

Božič

,

I.

Szleifer

, and

R.

Podgornik

, “

The role of solution conditions in the bacteriophage PP7 capsid charge regulation

,”

Biophys. J.

107

,

1970

–

1979

(

2014

).

https://doi.org/10.1016/j.bpj.2014.08.032

Google Scholar

Crossref

PubMed

11.

M.

Lund

, “

Electrostatic chameleons in biological systems

,”

J. Am. Chem. Soc.

132

,

17337

–

17339

(

2010

).

https://doi.org/10.1021/ja106480a

Google Scholar

Crossref

PubMed

12.

A.

Majee

,

M.

Bier

,

R.

Blossey

, and

R.

Podgornik

, “

Charge regulation radically modifies electrostatics in membrane stacks

,”

Phys. Rev. E

100

,

050601

(

2019

).

https://doi.org/10.1103/physreve.100.050601

Google Scholar

Crossref

PubMed

13.

R.

Netz

, “

Charge regulation of weak polyelectrolytes at low- and high-dielectric-constant substrates

,”

J. Phys.: Condens. Matter

15

,

S239

(

2002

).

https://doi.org/10.1088/0953-8984/15/1/331

Google Scholar

Crossref

14.

J. N.

Israelachvili

, “

Units, symbols, useful quantities and relations

,” in

Intermolecular and Surface Forces

,

3rd ed.

(

Academic Press

,

San Diego, CA

,

2011

).

Google Scholar

15.

S. H.

Behrens

,

D. I.

Christl

,

R.

Emmerzael

,

P.

Schurtenberger

, and

M.

Borkovec

, “

Charging and aggregation properties of carboxyl latex particles: Experiments versus DLVO theory

,”

Langmuir

16

,

2566

–

2575

(

2000

).

https://doi.org/10.1021/la991154z

Google Scholar

Crossref

16.

I.

Palaia

,

A.

Goyal

,

E.

Del Gado

,

L.

Šamaj

, and

E.

Trizac

, “

Like-charge attraction at the nanoscale: Ground-state correlations and water destructuring

,”

J. Phys. Chem. B

126

,

3143

–

3149

(

2022

).

https://doi.org/10.1021/acs.jpcb.2c00028

Google Scholar

Crossref

PubMed

17.

M.

Lund

and

B.

Jönsson

, “

Charge regulation in biomolecular solution

,”

Q. Rev. Biophys.

46

,

265

–

281

(

2013

).

https://doi.org/10.1017/s003358351300005x

Google Scholar

Crossref

PubMed

18.

M.

Lund

and

B.

Jönsson

, “

On the charge regulation of proteins

,”

Biochemistry

44

,

5722

–

5727

(

2005

).

https://doi.org/10.1021/bi047630o

Google Scholar

Crossref

PubMed

19.

J.

Landsgesell

,

L.

Nová

,

O.

Rud

,

F.

Uhlík

,

D.

Sean

,

P.

Hebbeker

,

C.

Holm

, and

P.

Košovan

, “

Simulations of ionization equilibria in weak polyelectrolyte solutions and gels

,”

Soft Matter

15

,

1155

–

1185

(

2019

).

https://doi.org/10.1039/c8sm02085j

Google Scholar

Crossref

PubMed

20.

N. A.

Pérez-Chávez

,

A. G.

Albesa

, and

G. S.

Longo

, “

Thermodynamic theory of multiresponsive microgel swelling

,”

Macromolecules

54

,

2936

–

2947

(

2021

).

https://doi.org/10.1021/acs.macromol.0c02885

Google Scholar

Crossref

21.

Y.

Burak

and

R. R.

Netz

, “

Charge regulation of interacting weak polyelectrolytes

,”

J. Phys. Chem. B

108

,

4840

–

4849

(

2004

).

https://doi.org/10.1021/jp036367i

Google Scholar

Crossref

22.

T.

Hiemstra

,

W. v.

Riemsdijk

, and

M.

Bruggenwert

, “

Proton adsorption mechanism at the gibbsite and aluminium oxide solid/solution interface

,”

Neth. J. Agric. Sci.

35

,

281

–

293

(

1987

).

https://doi.org/10.18174/njas.v35i3.16725

Google Scholar

Crossref

23.

T.

Hiemstra

and

W.

Van Riemsdijk

, “

Physical chemical interpretation of primary charging behaviour of metal (hydr)oxides

,”

Colloids Surf.

59

,

7

–

25

(

1991

).

https://doi.org/10.1016/0166-6622(91)80233-e

Google Scholar

Crossref

24.

A.

Bakhshandeh

,

A. P.

Dos Santos

, and

Y.

Levin

, “

Interaction between charge-regulated metal nanoparticles in an electrolyte solution

,”

J. Phys. Chem. B

124

,

11762

–

11770

(

2020

).

https://doi.org/10.1021/acs.jpcb.0c09446

Google Scholar

Crossref

PubMed

25.

A.

Bakhshandeh

,

D.

Frydel

, and

Y.

Levin

, “

Theory of charge regulation of colloidal particles in electrolyte solutions

,”

Langmuir

38

,

13963

(

2022

).

https://doi.org/10.1021/acs.langmuir.2c02313

Google Scholar

Crossref

PubMed

26.

Y.

Levin

, “

Electrostatic correlations: From plasma to biology

,”

Rep. Prog. Phys.

65

,

1577

(

2002

).

https://doi.org/10.1088/0034-4885/65/11/201

Google Scholar

Crossref

27.

Y.

Levin

and

M. E.

Fisher

, “

Criticality in the hard-sphere ionic fluid

,”

Physica A

225

,

164

–

220

(

1996

).

https://doi.org/10.1016/0378-4371(95)00336-3

Google Scholar

Crossref

28.

C.

Labbez

and

B.

Jönsson

,

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

(

Springer

,

2006

), pp.

66

–

72

.

Google Scholar

29.

C.

Labbez

,

B.

Jönsson

,

M.

Skarba

, and

M.

Borkovec

, “

Ion−ion correlation and charge reversal at titrating solid interfaces

,”

Langmuir

25

,

7209

–

7213

(

2009

).

https://doi.org/10.1021/la900853e

Google Scholar

Crossref

PubMed

30.

Y.

Levin

and

A.

Bakhshandeh

, “

Comment on “Simulations of ionization equilibria in weak polyelectrolyte solutions and gels” by J. Landsgesell, L. Nová, O. Rud, F. Uhlík, D. Sean, P. Hebbeker, C. Holm and P. Košovan, Soft Matter, 2019, 15, 1155–1185

,”

Soft Matter

19

,

3519

–

3521

(

2023

).

https://doi.org/10.1039/d2sm01393b

Google Scholar

Crossref

PubMed

31.

W.

Smith

and

B.

Triska

, “

The reaction ensemble method for the computer simulation of chemical and phase equilibria. I. Theory and basic examples

,”

J. Chem. Phys.

100

,

3019

–

3027

(

1994

).

https://doi.org/10.1063/1.466443

Google Scholar

Crossref

32.

J. K.

Johnson

,

A. Z.

Panagiotopoulos

, and

K. E.

Gubbins

, “

Reactive canonical Monte Carlo: A new simulation technique for reacting or associating fluids

,”

Mol. Phys.

81

,

717

–

733

(

1994

).

https://doi.org/10.1080/00268979400100481

Google Scholar

Crossref

33.

C.-Y.

Leung

,

L. C.

Palmer

,

S.

Kewalramani

,

B.

Qiao

,

S. I.

Stupp

,

M.

Olvera de la Cruz

, and

M. J.

Bedzyk

, “

Crystalline polymorphism induced by charge regulation in ionic membranes

,”

Proc. Natl. Acad. Sci. U. S. A.

110

,

16309

–

16314

(

2013

).

https://doi.org/10.1073/pnas.1316150110

Google Scholar

Crossref

PubMed

34.

G. M.

Ong

,

A.

Gallegos

, and

J.

Wu

, “

Modeling surface charge regulation of colloidal particles in aqueous solutions

,”

Langmuir

36

,

11918

–

11928

(

2020

).

https://doi.org/10.1021/acs.langmuir.0c02000

Google Scholar

Crossref

PubMed

35.

H.

Wennerström

,

B.

Jönsson

, and

P.

Linse

, “

The cell model for polyelectrolyte systems. Exact statistical mechanical relations, Monte Carlo simulations, and the Poisson–Boltzmann approximation

,”

J. Chem. Phys.

76

,

4665

–

4670

(

1982

).

https://doi.org/10.1063/1.443547

Google Scholar

Crossref

36.

J.

Landsgesell

,

P.

Hebbeker

,

O.

Rud

,

R.

Lunkad

,

P.

Košovan

, and

C.

Holm

, “

Grand-reaction method for simulations of ionization equilibria coupled to ion partitioning

,”

Macromolecules

53

,

3007

–

3020

(

2020

).

https://doi.org/10.1021/acs.macromol.0c00260

Google Scholar

Crossref

37.

T.

Curk

,

J.

Yuan

, and

E.

Luijten

, “

Accelerated simulation method for charge regulation effects

,”

J. Chem. Phys.

156

,

044122

(

2022

).

https://doi.org/10.1063/5.0066432

Google Scholar

Crossref

PubMed

38.

T.

Bollhorst

,

K.

Rezwan

, and

M.

Maas

, “

Colloidal capsules: Nano- and microcapsules with colloidal particle shells

,”

Chem. Soc. Rev.

46

,

2091

–

2126

(

2017

).

https://doi.org/10.1039/c6cs00632a

Google Scholar

Crossref

PubMed

39.

F. G.

Donnan

, “

Theorie der membrangleichgewichte und membranpotentiale bei vorhandensein von nicht dialysierenden elektrolyten. Ein beitrag zur physikalisch-chemischen physiologie

,”

Z. Elektrochem. Angew. Phys. Chem.

17

,

572

–

581

(

1911

).

https://doi.org/10.1002/bbpc.19110171405

Google Scholar

40.

S.

Barr

and

A.

Panagiotopoulos

, “

Grand-canonical Monte Carlo method for Donnan equilibria

,”

Phys. Rev. E

86

,

016703

(

2012

).

https://doi.org/10.1103/physreve.86.016703

Google Scholar

Crossref

41.

R. v.

Roij

, “

Defying gravity with entropy and electrostatics: Sedimentation of charged colloids

,”

J. Phys.: Condens. Matter

15

,

S3569

(

2003

).

https://doi.org/10.1088/0953-8984/15/48/016

Google Scholar

Crossref

42.

A. P.

Philipse

, “

Remarks on the Donnan condenser in the sedimentation–diffusion equilibrium of charged colloids

,”

J. Phys.: Condens. Matter

16

,

S4051

(

2004

).

https://doi.org/10.1088/0953-8984/16/38/020

Google Scholar

Crossref

43.

P.

Warren

, “

Electrifying effects in colloids

,”

Nature

429

,

822

(

2004

).

https://doi.org/10.1038/429822a

Google Scholar

Crossref

PubMed

44.

T.

Eckert

,

M.

Schmidt

, and

D.

de las Heras

, “

Sedimentation of colloidal plate-sphere mixtures and inference of particle characteristics from stacking sequences

,”

Phys. Rev. Res.

4

,

013189

(

2022

).

https://doi.org/10.1103/physrevresearch.4.013189

Google Scholar

Crossref

45.

M.

Sullivan

,

K.

Zhao

,

C.

Harrison

,

R. H.

Austin

,

M.

Megens

,

A.

Hollingsworth

,

W. B.

Russel

,

Z.

Cheng

,

T.

Mason

, and

P.

Chaikin

, “

Control of colloids with gravity, temperature gradients, and electric fields

,”

J. Phys.: Condens. Matter

15

,

S11

(

2002

).

https://doi.org/10.1088/0953-8984/15/1/302

Google Scholar

Crossref

46.

E. R.

Smith

, “

Electrostatic energy in ionic crystals

,”

Proc. R. Soc. London, Ser. A

375

,

475

–

505

(

1981

).

https://doi.org/10.1098/rspa.1981.0064

Google Scholar

Crossref

47.

D.

Frenkel

and

B.

Smit

,

Understanding Molecular Simulation: From Algorithms to Applications

(

Elsevier

,

2001

), Vol.

1

.

Google Scholar

48.

M. P.

Allen

and

D. J.

Tildesley

,

Computer Simulation of Liquids

(

Oxford University Press

,

2017

).

Google Scholar

Crossref

49.

Z.

Hu

, “

Infinite boundary terms of Ewald sums and pairwise interactions for electrostatics in bulk and at interfaces

,”

J. Chem. Theory Comput.

10

,

5254

–

5264

(

2014

).

https://doi.org/10.1021/ct500704m

Google Scholar

Crossref

PubMed

50.

A. P.

dos Santos

,

M.

Girotto

, and

Y.

Levin

, “

Simulations of Coulomb systems with slab geometry using an efficient 3d Ewald summation method

,”

J. Chem. Phys.

144

,

144103

(

2016

).

https://doi.org/10.1063/1.4945560

Google Scholar

Crossref

PubMed

51.

S.

Yi

,

C.

Pan

, and

Z.

Hu

, “

Note: A pairwise form of the Ewald sum for non-neutral systems

,”

J. Chem. Phys.

147

,

126101

(

2017

).

https://doi.org/10.1063/1.4998320

Google Scholar

Crossref

PubMed

52.

R.

Euwema

and

G.

Surratt

, “

The absolute positions of calculated energy bands

,”

J. Phys. Chem. Solids

36

,

67

–

71

(

1975

).

https://doi.org/10.1016/0022-3697(75)90115-8

Google Scholar

Crossref

53.

S.

De Leeuw

and

J.

Perram

, “

Computer simulation of ionic systems. Influence of boundary conditions

,”

Physica A

107

,

179

–

189

(

1981

).

https://doi.org/10.1016/0378-4371(81)90031-5

Google Scholar

Crossref

54.

J. S.

Høye

and

E.

Lomba

, “

Mean spherical approximation (MSA) for a simple model of electrolytes. I. Theoretical foundations and thermodynamics

,”

J. Chem. Phys.

88

,

5790

–

5797

(

1988

).

https://doi.org/10.1063/1.454753

Google Scholar

Crossref

55.

N. F.

Carnahan

and

K. E.

Starling

, “

Equation of state for nonattracting rigid spheres

,”

J. Chem. Phys.

51

,

635

–

636

(

1969

).

https://doi.org/10.1063/1.1672048

Google Scholar

Crossref

56.

A.

Bakhshandeh

,

D.

Frydel

, and

Y.

Levin

, “

Reactive Monte Carlo simulations for charge regulation of colloidal particles

,”

J. Chem. Phys.

156

,

014108

(

2022

).

https://doi.org/10.1063/5.0077956

Google Scholar

Crossref

PubMed

57.

E. A.

Guggenheim

, “

On the conception of electrical potential difference between two phases. II

,”

J. Phys. Chem.

34

,

1540

–

1543

(

1930

).

https://doi.org/10.1021/j150313a014

Google Scholar

Crossref

58.

S.

Craxford

,

O.

Gatty

, and

T.

Teorell

, “

XCV. A note on surface ph

,”

London, Edinburgh Dublin Philos. Mag. J. Sci.

25

,

1061

–

1066

(

1938

).

https://doi.org/10.1080/14786443808562090

Google Scholar

Crossref

2023

Author(s)

A new method for reactive constant pH simulations

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

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

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

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

A new method for reactive constant pH simulations

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

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

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

REFERENCES

Citing articles via

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only