Primitive model electrolytes. A comparison of the HNC approximation for the activity coefficient with Monte Carlo data

Accuracy of the mean activity coefficient expression (Hansen-Vieillefosse-Belloni equation), valid within the hypernetted chain (HNC) approximation, was tested in a wide concentration range against new Monte Carlo (MC) data for +1:-1 and +2:-2 primitive model electrolytes. The expression has an advantage that the excess chemical potential can be obtained directly, without invoking the time consuming Gibbs-Duhem calculation. We found the HNC results for the mean activity coefficient to be in good agreement with the machine calculations performed for the same model. In addition, the thermodynamic consistency of the HNC approximation was tested. The mean activity coefficients, calculated via the Gibbs-Duhem equation, seem to follow the MC data slightly better than the Hansen-Vieillefosse-Belloni expression. For completeness of the calculation, the HNC excess internal energies and osmotic coefficients are also presented. These results are compared with the calculations based on other theories commonly used to describe electrolyte solutions, such as the mean spherical approximation, Pitzer's extension of the Debye-H\"uckel theory, and the Debye-H\"uckel limiting law.


Introduction
Knowledge of the physico-chemical properties of electrolyte solutions contributes toward better understanding of the processes in biology, chemical sciences, and various technologies. It is of no surprise that electrolyte solutions are, after a hundred years, still extensively studied in search of theories providing quantitative description of these systems. The advances in modeling, together with the review of experimental data for electrolyte solutions, are summarized in books and thematic papers (see, for example, references [1][2][3][4][5][6]). Although it has been established that only the models that explicitly include solvent are capable of providing a proper microscopical description of electrolyte solutions [7,8], one can describe thermodynamic properties of these systems on the McMillan-Mayer level of approximation by a proper choice of ion size parameters [9,10]. In this paper we will restrict ourselves to the latter case.
Seminal efforts to understand the properties of electrolyte solutions were made by Debye and Hückel (DH) [11]. The DH theory is based on a linearized version of the Poisson-Boltzmann equation which yields important insights into dilute electrolyte solutions. Unfortunately, this approach can only be used to quantitatively describe very dilute solutions of size symmetric +1:−1 electrolytes [12]. The deficiencies of the Debye and Hückel theory were analyzed in many contributions (see, e.g., references [1,3,[13][14][15][16]). Various continuations were suggested to extend the range of applicability of the theory, among these the one proposed by Pitzer [17] is very useful for engineering purposes. This approach is simple, fully analytical, and capable of semi-quantitatively describing the osmotic coefficient behavior of +1:−1 electrolyte up to the concentration close to 1.0 mol dm −3 [17].
Theoretical extension, correcting for the deficiency of the Debye and Hückel theory, is provided by the modified Poisson-Boltzmann approach [14,[18][19][20][21]. For example, the individual activity coefficients of pure electrolytes were calculated by Molero et al. (cf. reference [22]) and for ternary systems containing neutral hard spheres recently by Outhwaite et al. [23]. The theory yields excellent agreement with computer simulations.
Another approach to calculate the properties of solutions involve a class of integral equation theories based on the Ornstein-Zernike (OZ) integral equation [24,25] and approximate closures. The so-called mean spherical approximation (MSA) was solved analytically for the primitive model electrolytes, with the thermodynamic quantities written in a closed form [26][27][28][29][30][31][32][33]. Through the energy route, one obtains the osmotic and activity coefficients that are in good agreement with Monte Carlo computer simulations for +1:−1 electrolytes [12]. Due to its analytical nature, the MSA theory provides useful insights and is accordingly extensively used (see, e.g., references [34][35][36][37][38][39][40]). Another approximate but very robust closure of the Ornstein-Zernike theory was proposed by Kovalenko and Hirata (KH) [41,42]. By changing the radii of the ions [32,[43][44][45] and/or including the association between unlike ions [46,47] these methods provide good fits to experimental data. An integral equation approximation that is widely used in describing the thermodynamics of symmetric, as well as highly asymmetric electrolytes, is the hypernetted-chain (HNC) closure [12,24,48,49].
An advantage of the OZ hypernetted-chain approximation is that it yields accurate structural description of electrolyte solutions in terms of the pair-distribution functions, g ij (r). Once this information is known, the standard statistical-mechanical equations, connecting the pair-distribution functions with thermodynamic properties, can be applied [24]. The reduced excess internal energy of the system is readily obtained via the expression where β = 1/k B T (k B being the Boltzmann constant and T the absolute temperature), ρ i (ρ j ) are the number densities of the species i (j), ρ is the total number density of the system, u ij (r) the interaction pair potential (see section 2), and dr = 4πr 2 dr. The osmotic coefficient, Φ HNC , can be conveniently calculated using the virial route [24]. For the primitive electrolyte, the pair potential contains a discontinuity at the distance of closest approach of two ions (see section 2). The problem can be surmounted by the introduction of the so-called background correlation (referred to as the cavity distribution) function, y ij (r) = exp[βu ij (r)]g ij (r), and the formal division of the u ij (r) to the Coulombic, u C ij (r), and the hard sphere part. In this way, the equation gets the form applicable to the primitive model electrolyte [24] Here a ij denotes the distance of closest approach of two spherical particles i and j, and u C ij (r) ∝ r −1 for all r, is the usual Coulomb potential.
The mean activity coefficient, γ ± , can be obtained by integration of the Gibbs-Duhem equation [2] ln γ ±,GD = Φ HNC − 1 + c c=0 for example DHLL+B2, at very low electrolyte concentrations [50]. The mean activity coefficient calculation via the Gibbs-Duhem equation becomes particularly cumbersome and time consuming in ternary mixtures containing additional solutes [51]. The expressions written above are general and can be in principle, apart from the HNC theory, applied to any other approximate closure.
An alternative way of calculating the activity coefficients, but valid only within the HNC approximation, has been proposed by Verlet and Levesque [52]. In its current form the formula was written by Hansen and Vieillefosse [53], and was successfully applied by Belloni [54] for asymmetric electrolytes. The expression, here referred to as the Hansen-Vieillefosse-Belloni (HVB) equation, reads Here γ i is the individual activity coefficient of the species i (+ or −), h(r) denotes the total, and c(r) the direct correlation function, while c (s) (0) denotes the Fourier transform of the short-range part of the direct correlation function at k = 0. Further, ρ j is the number density of species j.
To our best knowledge, the accuracy of the activity coefficient expression given by HVB equation (1.4) has not been thoroughly tested for the primitive model electrolytes. In a few cases, where the Monte Carlo data for the same system were available [55], the expression was found to be in good agreement with the "exact" machine calculations. Therefore, the purpose of this work is to systematically examine the validity of HVB equation (1.4) for the size symmetric and asymmetric +1:−1, and size symmetric +2:−2 electrolytes. Notice that the formula given by equation (1.4) is only valid within the HNC approximation and therefore is not generally applicable. In other words, different expressions for the mean activity coefficients can be derived consistently with the closure conditions used [56].
Even so, the HVB equation is very useful because it avoids long and cumbersome evaluation of the excess chemical potential of the solute via the Gibbs-Duhem route (see, e.g., reference [51]). Notice that the mean activity coefficient of bulk electrolyte is often an input information for studying heterogeneous systems, electrical double-layer, or Donnan equilibrium, i.e. wherever the bulk electrolyte is in equilibrium with charged surfaces. The results obtained by means of HNC approximation, coupled with the HVB equation, are in the present study compared with the new Monte Carlo simulations, and with the results of some other electrolyte theories used in describing electrolyte solutions. In addition, the internal consistency of the HNC theory is tested with respect to the two different routes of the mean activity coefficient calculation.
The paper is organized as follows: following a brief introduction, the model and methods are presented, followed by the results and discussion section. Numerical results are presented in form of figures and tables. Conclusions are given at the end. Appendices summarize all the relevant equations used in these calculations: Debye-Hückel theory (appendix A), Pitzer's approach (appendix B), and the mean spherical approximation study (appendix C).

Model and methods
The model most frequently used in describing the electrolyte solutions is called the primitive model. In this description, an ion is represented as a hard sphere carrying positive or negative charge in its center, while the solution as a whole is treated as a dielectric continuum with a relative permittivity ε r (McMillan-Mayer level of approximation) of pure solvent at pressure (p) and temperature (T ) of observation [24]. Despite its simplicity, the model is capable of explaining many experimentally determinable properties of real electrolytes [9,57], as well as their mixtures [10].
The interaction pair potential between two ions of valence z i (z j ), separated by a distance r, contains the hard sphere and Coulomb part, u C ij (r)

33003-3
where a ij = 1 2 (a i + a j ), a i (a j ) is the diameter of species i (j), and λ B is the Bjerrum length Here β = 1/k B T , where k B is the Boltzmann constant, T is the absolute temperature, ε 0 is the permittivity in vacuum and e 0 is the elementary charge. For aqueous solutions at 25 • C studied here, the Bjerrum length assumes the value λ B = 7.14 Å. Three different models of aqueous electrolyte solutions at 25 • C were considered: These values are often used as representatives for simple ions in aqueous solutions.

Monte Carlo simulation
The Monte Carlo calculations were performed in the canonical and grand canonical (GCMC) ensemble, using the standard Metropolis sampling algorithm [58]. Periodic boundary conditions with the minimum image (MI) convention, and Ewald summation (ES) method were used in the canonical ensemble simulations to minimize the finite sample effects, while GCMC was performed only in combination with the Ewald summation. The calculations were carried out with equal number of cations and anions, the total number of particles N being between 200 and 1000. After an equilibration run of (1 − 10) · 10 7 configurations, each production run consisted of (1 − 10) · 10 7 attempted configurations. From four to ten independent runs were performed for simulation in canonical ensemble for each concentration studied. The results given in tables and figures are the average values over these runs, after the equilibration runs were discarded.
While the excess internal energy and osmotic coefficient were obtained as simple canonical ensemble averages [58], the Widom's test particle insertion method [59, 60] was used in the canonical ensemble simulation to calculate the excess chemical potential of the electrolyte in solution In the expression above, · · · denotes the canonical ensemble average, while U pair is the interaction energy between the system and a non-perturbing test pair of oppositely charged ions, inserted at random locations in the system. During the canonical ensemble simulation, for every N (number of particles in the simulation box) configuration, a hundred of (Widom's) insertions were attempted, to obtain the canonical ensemble average requested by equation (2.2). The results for the mean chemical potential of the electrolyte solution as obtained by Widom's method may, at higher concentrations, depend on the number of particles in the system [61]. For this reason, in many cases, especially for more concentrated electrolytes, the canonical ensemble simulations were supplemented by the grand canonical ensemble Monte Carlo method. In the latter case, the excess chemical potential is obtained directly, without invoking the insertion method [62,63].

Hypernetted chain approximation
The hypernetted chain (HNC) approximation is based on the Ornstein-Zernike (OZ) equation. For multi-component mixtures, the OZ equation reads [24] where h(r) and c(r) are the total and direct correlation functions, respectively, and the integral is of the convolution type. A general closure relation between h(r) and c(r) for the OZ equation is [48] ln [h ij (r) is in literature known as the "bridge graph" and cannot be written as a closed form function of the distribution functions h(r) and c(r). B ij (r) are all the graphs in the representation of the background correlation function, which are neither series nor parallel graphs. At least two field points are bridged, i.e. connected by the Mayer bond. The HNC approximation assumes that these graphs mutually cancel and sets B ij (r) to zero. Due to the long-range nature of the Coulomb interaction, the set of equations (2.3) and (2.4) has to be re-normalized before it can be solved numerically [48,64]. The re-normalized form of the integral equation was solved by direct iteration using the fast Fourier transform routine on a linear grid with 2 18 division points separated by the distance of ∆r=0.005 Å. It should be noted that the results strongly depend on the number of points and separation interval ∆r. The decisive criteria in our case was the smallest zeroth (electroneutrality condition) and second moment condition that should both in theory be equal to zero [65,66].
Equations (1.1) and (1.2) were used to calculate the excess internal energy and the osmotic coefficients, respectively, while the equation (1.4) was used to calculate the activity coefficient. The mean activity coefficient, γ ± , for charge symmetric electrolytes studied here, is finally calculated from the expression γ 2 ± = γ + γ − . In order to check the thermodynamic consistency of the HNC approximation, the mean activity coefficients were (for a limited number of cases) calculated with the help of the Gibbs-Duhem equation (1.3), using the HNC osmotic coefficient (equation (1.2)) results.

Results and discussion
Thermodynamic properties, e.g. excess internal energies, osmotic coefficients, and activity coefficients, are for three examples: The excess internal energy, βE ex /N , chemical potential, ln γ ± , and osmotic coefficient, Φ, at different concentrations of size symmetric +1:−1 electrolyte are presented in tables 1 and 2. While all the GCMC results were obtained using the Ewald summation technique (table 1), the comparison between the canonical Monte Carlo results obtained using the periodic boundary conditions with minimum image convention and results obtained using Ewald summation are given in table 2. For simulations in the canonical ensemble, the Widom method was used to calculate the mean activity coefficients. As seen from table 2, for +1:−1 electrolyte in the concentration range studied here (from 0.0001 mol dm −3 to 1.5 mol dm −3 ), both ways of accounting for the finite sample effects (minimum image, Ewald summation), yield the results, which agree within 5% or (most often) better. This appears to be within the sum of numerical uncertainties of separate calculations. Figure 1 shows the results presented in tables 1 and 2 (ln γ ± and 1 − Φ; i.e. deviations from ideality), as well as the results obtained by the frequently used analytical theories (DH approximation, Pitzer's approach, MSA theory -see equations in the appendices). ln γ ± (panel (a)) and 1 − Φ (panel (b)) are shown for different methods as a function of square root of the concentration. It is clear from this figure that the DH theory quantitatively describes the osmotic and activity coefficients of the model electrolyte only up to approximately c = 0.05 mol dm −3 . Pitzer's continuation extends the range of the DH approach up to 0.75 mol dm −3 . By contrast, the MSA and HNC approximations are equally good in the whole concentration range, only at concentrations higher than 1 M, the HNC performs slightly better.
The main purpose of this work is to test the validity of the HVB equation (1.4). First we checked the thermodynamic consistency of the HNC approach by comparing the mean activity coefficients obtained by the Gibbs-Duhem equation with those obtained directly (equation (1.4)). The former 33003-5 results in table 2 are denoted as γ ±,GD . One can see that the two sets of results agree very well. It is our impression, however, that the mean activity coefficients calculated via the Gibbs-Duhem equation are slightly closer to the MC data than those calculated by the HVB equation. The latter results for γ ± seem to be systematically too high, but the differences are small, most likely within the experimental uncertainties.
The results for size asymmetric +1:−1 electrolyte are presented in tables 3 and 4 and figure 2. Notice that in the DH theory calculations the parameter a was taken to be a = 1 2 (a + + a − ) since the DH theory cannot take the size asymmetry of the ions into account. The conclusions are the same as for the case (a): the size asymmetry considered here does not dramatically effect the performance of these theories in the applied concentration range.  The excess internal energy, βE ex /N , excess chemical potential, ln γ ± , and osmotic coefficient, Φ, for different concentrations of size symmetric +2:−2 electrolyte are presented in tables 5 and 6 and figure 3 (only ln γ ± and Φ). Interestingly, different treatments of the boundary conditions do not effect the results for the excess chemical potential obtained by Widom's method. This is not entirely true for the excess internal energy and osmotic coefficient calculations at concentrations 33003-6 above 1 mol dm −3 . The discrepancies may become quite large at concentrations ≈ 1.5 mol dm −3 , and the Ewald summation method should be used to obtain correct results. For this reason the results of computer simulations, using the minimum image convention at 1.5 mol dm −3 , are not included in table 6. At such a high concentration of the +2:−2 electrolyte, the minimum image method yields unphysical pair distribution functions, and consequently should not be used for accurate simulations. As seen from figure 3, the differences among different theories become, as expected, more pronounced in the case of +2:−2 electrolyte than for +1:−1 electrolyte examined before. The results presented in this paper suggest, in close agreement with many previous calculations, that the HNC approximation accurately describes the thermodynamic properties of the primitive model electrolytes in a wide concentration range.

Conclusions
The expression for the mean activity coefficient valid within the HNC approximation, Hansen-Vieillefosse-Belloni equation, was for primitive model +1:−1 and +2:−2 electrolytes tested against the newly obtained Monte Carlo simulation data. For the sake of completeness, the thermodynamic properties of +1:−1 and +2:−2 electrolytes, calculated by means of some other theories, often used in describing the electrolyte solutions, are also presented.
Special attention is paid to the numerical accuracy of simulations and the HNC calculations. Although the HVB equation (expression (1.4)) has been used before, to our best knowledge, a systematic test of its accuracy has not been performed so far. For all the examples studied here, the HVB equation, compared to the canonical or grand canonical ensemble Monte Carlo simulations, yields accurate results for the mean activity coefficients and is, especially at higher concentration and +2:−2 electrolytes, superior to other approximations examined here. At the same time, it appears to be simple in usage and hence represents an excellent tool in describing the excess chemical potential of bulk electrolyte solutions. Notice, that this information is needed [55,70] whenever the membrane equilibria involving electrolyte solutions are studied.

A. Debye-Hückel theory
In the framework of the Debye-Hückel theory, the expression for the reduced excess internal energy reads βE ex where z + and z − are valencies of cations and anions, respectively, κ −1 is the so-called Debye screening length defined as κ 2 = 4πλ B ρ i z 2 i (λ B being the Bjerrum length and ρ i the number density of the ionic species i), and a is the distance of closest approach of two ions (assumed to be same for all pairs).
The osmotic coefficient is given by where σ(κa) is a function defined as [17,67] σ(x) = 3 The mean activity coefficient is calculated via The expression for the osmotic coefficient reduces in this case to [68]

B. Pitzer's approach
The osmotic coefficient is obtained via the virial route, which yields [17] (B1) In this expression ρ = ρ + + ρ − is the total number density of ionic species (for significance of other quantities see appendix A). An important improvement above the Debye-Hückel theory is the last term in equation (B1). The term provides the correction to the Debye-Hückel theory at intermediate concentrations and it causes the osmotic coefficient to increase at higher concentrations, exactly as observed experimentally. The equation for the mean activity coefficient is [17] ln γ ±,P = − λ B |z + z − | 6 where Φ P is given by the equation (B1).

C. Mean spherical approximation
Within the MSA closure, one calculates the excess internal energy as [30] βE ex MSA = − where α 2 = 4πλ B , ∆ = 1 − π 6 i ρ i a 3 i , and (C4) The electrostatic contribution to the osmotic coefficient is given by [30] Φ where ρ = i ρ i is the total density of the system. The electrostatic contribution to the activity coefficient of the species i has the form [30,33] ln γ el i = where b i = α 2 z i − π 2∆ a 2 i P n 2Γ(1 + Γa i ) , The hard sphere contribution to the osmotic and activity coefficients follow from the equation of state of a mixture of hard spheres. Here we use the Mansoori-Carnahan-Starling-Leland equation of state [69]. In this case, the osmotic coefficient reads [69] Φ hs = (1 + η + η 2 ) − 3η(y 1 + y 2 η) − η 3 y 3 ∆ 3 (C9) whereas the hard sphere contribution to the activity coefficient of the species i is [29] ln γ hs

33003-12
Notations in equations (C9) and (C10) denote the following , and as usual The osmotic and activity coefficients of the primitive model electrolyte solution are then obtained by summing up the electrostatic and hard sphere contributions ln γ i,MSA = ln γ el i + ln γ hs i .