Scientific paper Thermodynamics of Asymmetric Primitive Model Electrolytes via the Hypernetted Chain Approximation Tomaž Mohoric, Miha Luksic* and Barbara Hribar-Lee Faculty of Chemistry and Chemical Technology, University of Ljubljana, Aškerčeva c. 5, SI-1000 Ljubljana, Slovenia * Corresponding author: E-mail: miha.luksic@fkkt.uni-lj.si Received: 12-12-2011 Dedicated to Prof. Dr. Gorazd Vesnaver on the occasion of his 70h birthday Abstract The accuracy of the activity coefficient expression (Hansen-Vieillefosse-Belloni (HVB) equation), valid within the hypernetted-chain (HNC) approximation, was tested in a wide concentration range against newly obtained grand canonical Monte Carlo data for the size and charge asymmetric primitive model electrolytes. In some cases, uncharged hard sphere component was also present. The HVB expression enables a direct calculation of the excess chemical potential, without invoking the time consuming calculation via the Gibbs-Duhem relation. We found the Ornstein-Zernike (OZ)/HNC results for the mean activity coefficient, as well as for the reduced excess internal energy and osmotic coefficient, to be in good agreement with the machine calculations performed for the same model. The accuracy of the results was found to be dependent on the packing fraction of the solutions. The mean spherical approximation calculations were also used to describe the thermodynamics of these systems and compared with the OZ/HNC and simulation results. Keywords: Primitive model electrolyte, hypernetted-chain approximation, mean spherical approximation, Hansen-Vi-eillefosse-Belloni equation, mean activity coefficient, osmotic coefficient, Monte Carlo simulation 1. Introduction Electrolyte solutions play an important role in a number of different areas of chemistry, biology, and their related fields. Knowledge of ionic equilibria is of great significance in a variety of processes, ranging from elementary cases to technologically demanding operations. To be able to describe chemical equilibria of electrolyte solutions under different conditions (concentration, temperature etc.), it is necessary to know the mean activity coefficients under these conditions. It is therefore not surprising that the topic has been extensively studied; the advances in modelling, together with the review of experimental data for electrolyte solutions, are summarized in books and thematic papers (see, for example, references 1-7). Since it has been established that by a proper choice of the ion size parameters, the thermodynamic properties of electrolyte solutions can be adequately described even by simple models based on the McMillan-Mayer level of description,8,9 i.e. the solvent is treated implicitly, we will here restrict ourselves to this case. The most famous attempt to describe the thermodynamics of electrolyte solutions resulted in the Debye-Huckel theory (DH).10 The theory is based on a linearized Poisson-Boltzmann equation and 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.11 The deficiencies of the Debye-Huckel theory were analyzed in many contributions and various improvements were suggested to extend the range of applicability (for review see references 1,3,12). Theoretical extension, correcting for the deficiency of the DH theory, is provided by the modified Poisson-Boltzmann approach.1317 For example, the individual activity coefficients of pure electrolytes were calculated by Molero et al. (cf. reference 18) and for ternary systems containing uncharged hard spheres recently by Outhwaite et al.19 This theory yields excellent agreement with computer simulations. Another approach to calculate the properties of solutions involves a class of integral equation theories based on the Ornstein-Zernike (OZ) integral equation20 and approximate closures. The so-called mean spherical approximation (MSA) was, for primitive model electrolytes, solved analytically, with the thermodynamic quantities written in a closed form.21-28 Through the energy route, one obtains the osmotic and mean activity coefficients that are in good agreement with Monte Carlo computer simulations for +1:-1 electrolytes.11 A more recent approximate closure of the Ornstein-Zernike theory was proposed by Kovalenko and Hirata (KH).29-30 Another approximation that is widely used for describing the thermodynamics of symmetric, as well as of highly asymmetric electrolytes, is the so-called hypernet-ted-chain (HNC) closure.11,20,31-32 Solving the OZ equation within the HNC approximation provides the pair-distribution functions33 which contain all the information about the fluid's structure. Once this information is known, the standard statistical-mechanical equations connecting the pair-distribution functions with thermodyna-mic properties (i.e. reduced excess internal energy and osmotic coefficient) can be applied.20,33 The mean activity coefficient can be, in general, calculated from the concentration dependence of the osmotic coefficient, by integrating the Gibbs-Duhem equation.34 Although the numerical procedure is straightforward, it may be time consuming.2 An alternative way of calculating the activity coefficients, although valid only within the HNC approximation, was proposed by Verlet and Levesque.35 The equation was written in the current form by Hansen and Vieillefosse36 and successfully applied to asymmetric electrolytes by Belloni.37 It is therefore referred as the Hansen-Vieillefos-se-Belloni (HVB) equation:12 lay;,HVB = /ciJ)Wdr+ i ) ¡{hijCr^hijir)- c,y(r)]}dr (1) Here y is the individual activity coefficient of species i, Pj is the number density of j-th species, h and c denote the total and the direct correlation functions, respectively, while c(s) stands for the short-range part of c. The distance between particles i and j is given by r, and the volume element dr equals 4m2 dr due to spherical symmetry. In this work, we were primarily concerned with two issues: First, what is the validity of the OZ/HNC theory for asymmetric primitive model electrolytes. We considered the cases of different asymmetric electrolytes and compared the values of the reduced excess internal energies and osmotic coefficients predicted by the OZ/HNC theory with the data obtained from Monte Carlo computer simulations. The latter data represent the exact values for a given model. Second, we wanted to check the validity of the HVB equation for the above mentioned cases of electrolytes. Both issues were recently considered in our labo- ratory for the cases of size symmetric and asymmetric +1:-1, and size symmetric +2:-2 electrolytes.12 In the last part of this work, we extended our study of electrolyte solutions also to the case of a mixture of the electrolyte with a neutral (uncharged) component (the so-called solvent primitive model for electrolyte solutions). It has been established earlier19,38 that the HNC approximation adequately describes such systems only for low volume fractions of added neutral component. We considered the case of a restricted primitive model electrolyte (i.e. a size symmetric +1:-1 electrolyte) in a mixture with uncharged hard spheres of the same size as ions, and examined the behaviour of the OZ/HNC theory for different concentrations of the electrolyte and the uncharged component. The paper is organized as follows: proceeding the brief introduction, the model and methods are presented, followed by the results and discussion section. Numerical results are presented in a form of tables and figures. Conclusions are given at the end. Appendix summarizes the relevant equations used for the MSA approximation (OZ/MSA theory). 2. Model and Methods 2.1. Primitive Model Electrolyte The so-called primitive model for electrolyte solutions considered in this work is widely used for it is very simple but still able to capture many experimental observations of single electrolyte solutions8,39 and of their mixtures.9,39 The ions in this model are represented as charged hard spheres while the medium (solvent) is assumed to be a dielectric continuum with relative permittivity of a true solvent at a given temperature and pressure.20 The pair interaction potential u.. is composed of a hard-sphere part and an electrostatic (Coulomb) part, and for ions i and j with nominal charges zi and z, and diameters ai and a, i J 1 J respectively, separated by a distance r it reads: pun = r < ai tJ 'pufj = ZiZj-f, r > a-ij k r (2) Here ajj = 1/2(ai + a) and XB is the Bjerrum length defined as: An = Pel 47T£0£r (3) where e0 is the fundamental unit of charge, e0 is the permittivity of the vacuum and er is the relative permittivity of a medium. P = 1/kBT, where kB is the Boltzmann constant and T the absolute temperature. The Bjerrum length at the temperature of 25 °C has the value of 7.14 A for aqueous solutions (er = 78.5). For the ionic diameters of size asymmetric +2:-1, +3:-1 and +4:-1 electrolytes we took the semi-empirical values as determined by J. P. Simonin et al.8 The values correspond to ionic sizes at infinite dilution. Since we compared only the theory and the computer simulations, there was no need to take experimentally fitted ionic sizes (they are actually changing with the concentration of the solution). To complete our study, we also examined the case of a size symmetric +2:-1 electrolyte and the case of a highly size and charge asymmetric +10:-1, and +10:-2 electrolyte. Modelled electrolytes and their ionic diameters are summarized in Table 1. In each case we surveyed the concentration range from 10-4 mol dm-3 up to 1.5 mol dm-3 (except when the upper concentration was limited by the occupation of space). Table 1. Modelled electrolytes, corresponding ionic diameters, and approximate concentration ranges studied. Electrolyte [A] a-[A] +2:-1 4.25 4.25 +2:-1 (MgCy 6.71 3.62 +3:-1 (LaCl3) 7.75 3.62 +4:-1 (ThCl4) 12.1 3.62 + 10:-1 20.0 4.25 + 10:-2 20.0 4.25 Conc. range [mol dm-3] ,-4 10 10-4 -"1-4 1.5 1.5 10-4- 1.5 10-4 - 0.5 5 • 10-3 - 5 • 10-2 2.5 • 10-2- 0.1 for electrolytes, as well as for the electrolytes in a mixture with a neutral component, is well established and extensively described in several previous papers, and is therefore not repeated here.41-43 The details of the simulations are as follows: The total number of particles in the simulations varied between 200 and 1200, depending on the electrolyte. For each concentration studied, the simulation consisted of an equilibration run plus four independent production runs. Each run consisted of 107 attempted configurations, where a new configuration was formed by either changing the position of a randomly chosen particle (canonical step), or insertion or annihilation of a neutral combination of particles (grand canonical step). The ratio between canonical and grand canonical steps in the simulation depended on the system studied. Averages of the ther-modynamic properties obtained from the production runs are collected in Tables 2-9. The errors given in parentheses were calculated as the maximum difference between the run's value in the series and the average of that series. 2. 3. Hypernetted-chain (HNC) Approximation Modern fluid theories are based on the Ornstein-Zernike equation that, for an isotropic multicomponent system, reads:20,33 ~ cij(rlj) + J ctk(rik)hkj(rkj)dTk (5) In the cases of mixtures of an electrolyte with an uncharged component we kept the primitive model description for the +1:-1 electrolyte. The uncharged component consisted of uncharged hard spheres. All particles (charged and uncharged) were of the same size (a+ = a = aN = 4.25 A). The pair interaction potential between ions remained the same as given above (equation 2). The expression for the pair interaction potential between the uncharged particle and any other particle is: ico, 0, r < ay r > a, j (4) We set the concentration of the electrolyte constant and were changing the concentration of the uncharged component. Concentration range studied was limited with the convergence of the computer simulations. 2. 2. Monte Carlo Computer Simulation We performed the Monte Carlo computer simulations in the grand canonical ensemble (GCMC), using the standard Metropolis sampling algorithm40 and the periodic boundary conditions.40 To take into account the long range nature of electrostatic forces we applied the Ewald summation technique.40 The methodology of the method As in equation 1, h denotes the total correlation function and c the direct correlation function, pk is the number density of a component k, and the sum goes over all the components in the mixture. The correlation functions h and c are both unknown. The general closure relation that represents the connection between the two reads:33 (6) where /3 = 1/kBT, uij is the pair potential energy between particles i and j separated by a distance r, and Bij is a set of integrals known as the "bridge graphs". Since Bj cannot be written as a closed form function of h and c, we are forced to use some approximation in order to solve the system of equations 5 and 6. The hypernetted chain approximation sets the Bit term equal to zero.3132 ij Due to the long range nature of electrostatic forces equations 5 and 6 need to be re-normalized, i.e. split into the long- and short-range part, before solved iteratively.13 The re-normalized form of the integral equation was in our case solved by the direct iteration using the fast Fourier transform routine on a linear grid with 218 division points separated by the distance of Ar = 0.005 A. Results depend strongly on these two parameters. The criteria of the smallest zeroth (electroneutrality condition) and second moment (in theory both moments should be equal to zero44-45) were used to determine the parameters. The numerical solution provides the h (and c) correlation functions from which the pair distribution functions g, can be obtained as g = h + 1. Thermodynamic properties (like reduced excess internal energy and osmotic coefficient) can be expressed in terms of these pair-distribution functions. The expression for the reduced excess internal energy per particle (N) of the system in terms of pair-distribution functions reads:33 (7) where pi and pj are the number densities of components i and j, p is the total number density of the system, and u. is the pair interaction potential. The osmotic coefficient, corresponding to the pair interaction potential given by equations (2) and (4), can be calculated via the virial route:20 (8) where a.. is the distance of closest approach of two particles i and j, and uC. is the Coulomb part of the pair interaction potential (see equation 2). For a general Av+ Bv electrolyte, the individual activity coefficients are calculated via equation 1 while the mean activity coefficient of the electrolyte solution is defined as yl++V- = Y++Yv-, where v+ and v are stechiometric factors in the formula of the electrolyte. 3. Results and Discussion Thermodynamic properties - reduced excess internal energies per particle, reduced excess chemical potentials (P/jex = lnY±), and the osmotic coefficients - are for various systems studied here presented in Tables 2-9 and Figures 1-5 (only lnY± and 1-0), as a function of the electrolyte concentration. Note that the molar concentration of species Cj is related to its number density in the following way: pi = ciNA, where NA is the Avogadro number. The results are grouped into four subsections (a-d) with respect to the size and charge asymmetry of the electrolyte, and regarding the presence of the uncharged (neutral) component. The last subsection (e) discusses the overall trends in OZ/HNC performance. (a) Size symmetric but charge asymmetric electrolytes (z+ = +2, z_ = -1, a+ = a_ = 4.25 A). The reduced excess internal energy per particle, jEEXIN, the excess chemical potential, ln/±, and the osmotic coefficient, 0 at different concentrations of a size symmetric +2:-1 electrolyte, c, are presented in Table 2 and in Figure 1 (ln/± and 1-0). In Figure 1, ln/± (panel (a)) and 1-0 (panel (b)) are shown for different theoretical methods as a function of the square root of the concentration: symbols represent the GCMC data, the continuous line corresponds to the OZIHNC results while the dashed line shows the results of the OZIMSA. The error bars in GCMC values roughly correspond to the size of the symbols and are therefore not shown explicitly. Figure 1. lny± (panel (a)) and 1-0 (panel (b)) as a function of c112 for a +2:-1 primitive model electrolyte, a+ = a = 4.25 Â. Symbols denote the grand canonical Monte Carlo data, and theoretical values are given by lines (continuous: OZIHNC, dashed: OZIMSA). Results apply to Ag = 7.14 Â. It is clear from Table 2, that the MSA and HNC approximations are equally good in describing the reduced excess internal energy, as well as the osmotic coefficient (see also Figure 1) of this electrolyte in the whole concentration range. Only at concentrations higher than 1 mol dm-3, the OZ/HNC performs slightly better. This confirms the previous observations that OZ/HNC describes the behaviour of simple electrolytes extremely well.11,20,31-32 e e o as CO vo vo t- in VO CO t- (N in CO 00 t- (N (N 00 t- VO (N o t- CO CO O 00 c^ 00 00 00 00 00 o o o O O o o o o o o o oooo «COrH Onrn IN in vo as in CO 00 I I 00 "^t VO OOO I I I Ocn»N coco-^in lllll as 00 (N t- (N o CO VO o 00 VO VO (N (N VO in as 00 t- VO (N O t- in CO as C^ 00 00 00 00 00 00 o O O O O O o o o o o O VOVOCO ohm (N CO as CO in in in (N 00 VO in O t- 00 OS (N CO CO o (N (N (N (N VO t- in o 00 CO 1—1 in T-H IN CO in vo t- 00 00 00 t- VO O o o o o O O O o o O o o llllllllllll CN o CN o «Ï ^ o (N ¿a CO o i 00 1 (N OS O O o CO in (N O 00 i—i t- 1—1 (N vo in OS CO 00 T-H IN co VO t- O O O O o o O CO 00 o 00 ^tOOVM OOin^n COCO-^in lllll lllll COOVO(NOO 0 00 CO p^ ci it IN o s_1 o\ o\ O i—i it 00 o\ 00 00 o\ o\ OS 00 as in o in o as as p^ o IN ■^t cJ p^ C3 C3 C3 The results of the HVB equation for the excess chemical potential of the electrolyte, i.e. ln/±, are also in excellent agreement with the computer simulations, and superior to the MSA results in the whole concentration range studied here (cf. Table 2 and panel (a) of Figure 1). (b) Size and charge asymmetric electrolytes (z+ = +2 (a+ = 6.71 À), +3 (a+ = 7.75 À), and +4 (a+ = 12.10 À), z =-1 (a = 3.62 À)). The ionic radii for the model electrolytes studied in this subsection were taken from reference 8 to capture the behaviour of MgCl2 (+2:-1), LaCl3 (+3:-1), and ThCl4 (+4:-1) in aqueous solutions. The results for the reduced excess internal energy per particle, fiEEX/N, the excess chemical potential, ln/±, and the osmotic coefficient, 0, at different concentrations of a size and charge asymmetric electrolytes (with the charge and size parameters given above) are presented in Tables 3-5 and Figure 2. Similar as in Figure 1, ln/±, (panel (a)) and 1-0 (panel (b)) are shown for different theoretical methods as a function of square root of the concentration: symbols represent the GCMC data, the continuous lines correspond to the OZ/HNC results while the dashed lines show the results of the OZ/MSA. Figure 2. Same as in Figure 1, but for size and charge asymmetric electrolytes. Symbols denote the grand canonical Monte Carlo data (+2:-1 (•), +3:-1 (O), and +4:-1 (x) electrolyte). O co 00 t- 00 >n o o t- (N co o co as \o co co 00 t- so co 00 r- r- O 00 as ^ as 00 00 00 00 c^ c^ o o o o o cá cá cá cá cá cá < ce 5 o o z a 5 e 5 ca o >n >n SO as i—i (N o . . . . . . O 00 (N CO en oo I I t- 00 co r- O, Tt en oo I I o o >n 00 (N co >n >n i—1 (N co >n i—1 >n o (N i—i CO O 1—1 o O >n i—i 00 (N c^ SO c^ . o o cá cá cá cá cá cá Ó O cá O llllllllllll co co o co cá cá cá (N 00 co O 00 co >n O co >n SO ^ co co SO 00 c^ SO SO 00 cá cá cá cá cá Ó Ó ó cá >n a\ >n 00 (N o >n o . . . . . . I lililí lllll OOOOOCOCOCOCOOO 00000000000 >n SO r- o (N (N O . . . . . . so o n >n 00 co co SO 00 o (N (N r- co o CO co o >n c^ SO so Ó cá cá cá cá cá cá Ó ó cá cá cá CN CN 0 o o >n o co 00 a\ co a\ . . . . . c^ . iiiii i i co co co c^ co c^ o 00 00 00 00 c^ ^ (N (N 00 >n (N (N o co 00 (N __1 00 o c^ c^ c^ c^ c^ 00 00 00 00 c^ c^ cá cá o cá o cá cá o cá O cá c^ c^ 5? c^ c^ c^ r- CD l/o l/o 1/0 1—1 o . . . c^ . l/o . 00 . CM Tt CO co I I c? I/o oí ó? ■ít I/o CD o ■ít I/o (N ■a I/o I/o so co o co o a\ so co 1—1 co l/o so so l/o cá cá o cá o 1 cá 1 cá Ó ó o cá 1 cá lili lili t- so o 00 I I ^co^oocooco I/O^O^^I/ooi^ i/o-^tococoi/ooo (NcO^i/O^O^OCOO dddddddo lililí c^ l/o a\ co CD 00 c^ O . . . . co . c^ . lllll 000000000 1/0 ^CO cocoi/OÍN^OCOI/O^OOCOI/OÍNCO <><><><>r>co^qco<>(Ncoi/oi>co 0>nr- 0nr-oincoco cmcovoo cicicici II iNiNO-^tO ^hcocovooo 7 7 7 7 7 o ta IN CO t- ta o\ O t- 00 00 00 o\ in 00 in IN o\ CO .—1 o IN SO o 00 00 00 00 00 00 ^ o o o ci ci ci ci ci ci o\ in CO .-H o\ CO SO SO SO o\ 00 IN 00 ta 00 IN t- ta. IN so ■—1 in .—1 r- i—| CM CO VO r- —. —. ^ VO i—i ci ci ci ci ci ci i i 1 ci i—H I III OCOOOfN Hrn^rnr-COOVOinO ^HCMCOVOOO cicicicici III (N (N CO "tf t-t- 00 VO 00 lllll s *—1 IN in ta. O t- o\ r- 00 in IN o\ CO o ^ 00 00 00 00 ci ci ci ci ci ci c^ VO ta CD IN in O VO 00 00 a\ CM ci ci ci i—H oocoo t- ^HCMCO^ cicici9 III1 ta. s C? ^ 00 O-^t^q CMCMCM I I I O^OCOO 0000 miflrH r-in-^t ^HCMCM I I I COfNCO ocoo cmcmcm I I I OOO 000^0;^ CM I I (N I rococo ^ooco ^HCMCM I I I CN up 7 o CN CN o o 'a 3 o "o e S e e co in a\ H in rn CO in co OOO q ^ >q in ^o ^o ^o II cK lar challenge since a neutral combination of ions to be inserted or deleted becomes 11 (6) in this particular case, and the probability for such an insertion is, depending on the polyelectrolyte concentration, extremely small (less than 1%). However, the convergence of the presented results was tested by repeating the same simulation with an increased number of attempted insertions. The results in both cases differed no more than different blocks of the same run. At very low concentrations (< 0.005 mol dm-3 for +10:-1, and < 0.01 mol dm-3 for +10:-2 electrolyte) the OZ/HNC did not provide a convergent solution for these systems, while at high concentrations (> 0.075 mol dm-3 for +10:-1, and > 0.0115 mol dm-3 for +10:-2 electrolyte) the ergodic condition was not satisfied for the simulation method used. The results are collected in Tables 6 and 7, and presented in graphical form in Figure 3. Again, symbols are used to present GCMC data: filled symbols correspond to the +10:-1 electrolyte, and empty Figure 4. ln^± (panel (a)) and 1-0 (panel (b)) in the case of a mixture of a +1:-1 electrolyte with uncharged hard sphere component as a function of the concentration of uncharged component, cN (aN = 4.23 A) The concentration of a +1:-1 electrolyte present in the systems (a+ = a = 4.23 A, XB = 7.14 A) was approximately 1.0 mol dm-3 in all cases. Symbols are grand canonical Monte Carlo data, and theoretical values are given by lines (continuous: OZ/HNC, dashed: OZ/MSA). < ITS O 5 O z a 5 6 CO. O O 5 ¿a CN CN CN 230% of the 0 value while for simple electrolytes only up to 50%. (d) Electrolyte in a mixture with uncharged component, i.e. the solvent primitive model (z+ = -z_ = +1, a+ = a = 4.25 A, zN = 0, aN = 4.25 A) In the last case studied, we tested the performance of the HNC approximation for the system where a neutral (uncharged) component was added to a size symmetric +1:_1 electrolyte of the approximate concentration cEL = 1.0 mol dm3. (Note that due to the use of the GCMC method, the concentration of the electrolyte fluctuated around the value of 1.00 mol dm3 for ± 0.07 mol dm3.) The same thermodynamic properties of the solution as given in previous subsections were calculated at different concentrations of the added neutral component, cN, modelled as hard spheres with diameter aN = 4.25 A. The results are presented in Tables 8 and 9, and in Figure 4. Table 9. The reduced excess chemical potential of the uncharged component, \nyN, as obtained from GCMC, OZ/HNC, and OZ/MSA theory for a model +1:_1 electrolyte in a mixture with uncharged component (a+ = a = aN = 4.25 A). ^ = 7.14 A. The concentration of the electrolyte varied within the range of 1.00 ± 0.07 mol dm-3. GCMC HNC MSA cN [mol dm 3] ln?N ln7N ln?N 0.1000(2) 0.445(2) 0.4447 0.4488 0.2472(3) 0.492(1) 0.4924 0.4837 0.501(1) 0.540(2) 0.5417 0.5451 0.752(1) 0.603(2) 0.6052 0.6073 1.004(3) 0.666(2) 0.6698 0.6713 1.255(2) 0.732(1) 0.7367 0.7367 1.506(2) 0.800(1) 0.8063 0.8037 2.516(4) 1.092(2) 1.106 1.091 5.11(1) 1.989(3) 2.069 1.983 7.76(5) 3.239(6) 3.480 3.204 10.53(7) 5.018(6) 5.675 4.974 It is clear from Tables 8 and 9 and Figure 4 that as the concentration of the neutral component in the system increases, the discrepancy of the HNC approximation results with respect to the GCMC data increases. This is consistent with previous observations19,38,43 that the HNC approximation adequately describes such systems only for low volume fractions of the added neutral component. The MSA approximation is superior to the HNC in the regions of high volume fraction of the neutral component and performs very well in the whole concentration region studied here. The contribution of the neglected "bridge graphs" in the HNC approximation becomes important under these conditions. It is worth mentioning that the calculations with the added neutral component were repeated also for a different electrolyte concentrations (approximately 0.01 mol dm-3) but since the trends in the thermodynamic functions hardly differ from the ones presented here we did not include them in this paper. In Figure 5, we present the dependence of the excess chemical potential of the uncharged component, lnY^, on the concentration of the uncharged species, cN. It is clear that the results of the OZ/MSA theory are in an excellent agreement with GCMC data in the whole concentration range studied, while OZ/HNC describes the ln^ rather poorly at high concentrations of the neutral component. However, at low concentrations of the uncharged component (cN < 5 mol dm 3) OZ/HNC predicts the ln^ very close to the simulation data. (e) The analysis of the performance of the HVB equation in the primitive model asymmetric electrolyte systems. From the results given above, it is clear that in different systems the range of the validity of the HNC approximation, as well as the range of applicability of the HVB equation, is different for each electrolyte studied, as well o---T-1-1-1-1- 0 1 4 6 3 8 10 12 i ^ [mol (l:ii ] Figure 5. ln^ as a function of cN (aN = 4.25 Â). The concentration of a +1:-1 electrolyte present in the systems (a+ = a = 4.25 Â, Ag = 7.14 Â) was approximately 1.0 mol dm 3 in all cases. Symbols are grand canonical Monte Carlo data, and theoretical values are given by lines (continuous: OZ/HNC, dashed: OZ/MSA). as it depends on the concentration of the neutral component added to the system. In attempt to determine the parameter that is decisive in the applicability of the method, the relative error in the osmotic coefficient and mean activity coefficient was calculated as a function of the packing fraction of the solution, ii = L ^a-- Here the sum runs over all the species in the solution, as before, is the diameter of the species i, and pi is its number density. The error was calculated with respect to the GCMC simulation results and is shown in Figure 6. Figure 6. Relative error in the OZIHNC results (with respect to the GCMC data) for the mean activity coefficient, Y± (panel (a)), and osmotic coefficient, 0 (panel (b)), as a function of the packing fraction of the solution, n Different symbols correspond to different cases studied: +1:-1 model electrolyte in a mixture with uncharged component (+), size symmetric +2:-1 (x), MgCl2 (□), LaCl3 ThCl4 (■), +10:-1 (•), and +10:-2 (O) model electrolyte. The relative error in the OZIHNC results for the mean activity coefficient (obtained through the HVB equation (equation 1)) and for the osmotic coefficient (equation 8) depends, regardless of the nature of the electrolyte or the presence of the neutral component, uniquely on the packing fraction of the solution. This trend in the error, displayed in Figure 6, suggests that the packing fraction must have a major impact on the performance of the HNC ap- proximation. These figures could also be used to estimate the relative error of OZ/HNC calculations under given conditions. Somehow large difference between OZ/HNC and GCMC value for the osmotic coefficient for the +10:-1 and +10:-2 electrolyte solutions is a consequence of a poorly determined osmotic coefficient in the GCMC, and not the inaccuracy of OZ/HNC. This issue can be solved by calculating the exact contact values of the pair distribution functions using the Widom technique.48 The electrostatic contribution to the individual activity coefficient of the species i reads:25'28 ci ■ ' 2AZ-; 1 + raj UZ, p =1 y ^ai " D.Aj 1 + Tat 4P Pi 11 2 r> 1 + ra; (10) (11) (12) The electrostatic contribution to the osmotic coefficient is given by:25 I ff~ /K ' ^ 0EL 3np 8p IA ) 3 np 8 p where p = Xp is the total density of the system. (13) 2rij, (16) The hard sphere contribution to the osmotic and activity coefficients follows from the equation of state of a mixture of hard spheres. Here we use the Mansoori-Car-nahan-Starling-Leland equation of state.47 In this case, the osmotic coefficient has the following form:47 HS _ (i + v + i?2) - 3iKyi + yin) - J?3y3 (17) A3 while the hard sphere contribution to the activity coefficient of the species i reads:24 + ^[3(1 - ad 4- fit + y (ai - ft - to - 1) Notations in equations 17 and 18 denote the follo wing: (18) àij = PiPj Jwij (ai ~ ajf P V psrai + ak. «i = — / T-A-, Pi Ja.ak ataj ik and as usual The osmotic and activity coefficients of the primitive electrolyte solutions are the summations of electrostatic and hard sphere contributions: Generally for the electrolyte, one calculates the mean activity coefficient of the solutions from the definition Yv++V- = f+yl-. 6. Acknowledgement M. L. and B. H. L. appreciate the support of the Ministry of Science and Education of Slovenia under grant P1-0201. 7. References 1. H. S. Harned, B. B. Owen, The Physical Chemistry of Electrolyte Solutions, Reinhold: New York, 1958. 2. R. A. Robinson, R. H. Stokes, Electrolyte Solutions, 2nd rev. ed., Dover Publications, Mineola, New York, 2002, pp. 432-455. 3. J. M. C. Barthel, H. Krienke, in: H. Baumgärtel, E. U. Franck, W. Grümbein (Ed.): Topics in Physical Chemistry, Vol. 5, Springer, New York, 1998. 4. J. O'M. Bockris, A. K. N. Reddy, Modern Electrochemistry, Vol. 1, Ionics, 2nd ed., Plenum Press, NY, 1998. 5. J. R. Loebe, M. D. Donohue, AIChE J. 1997, 43, 180-195. 6.G. Lamm, in: K. B. Lipkowitz, R. Larter, T. R. Cundari, in: Reviews in Computational Chemistry, vol. 19, Chapter 4, Wiley-VCH, John Wiley & Sons, Inc. 2003, pp. 147365. 7.L. L. Lee, Molecular Thermodynamics of Electrolyte Solutions, World Scientific, New Jersey 2008 8. J. P. Simonin, O. Bernard, L. Blum, J. Phys. Chem. B, 1998, 102, 4411-4417. 9. B. Hribar-Lee, V. Vlachy, J. Mol. Liq. 2005, 118, 163-169. 10. P. Debye, E. Hückel, Phys. Z. 1927, 24, 185-206. 11. J. P. Valleau, K. Cohen, D. N. Card, J. Chem. Phys. 1980, 72, 5942-5954. 12. E. Gutiérrez-Valladares, M. Luksic, B. Millán-Malo, B. Hri-bar-Lee, V. Vlachy, Cond. Matt. Phys. 2011, 14, 33003: 1-15. 13. K. Singer and C. W. Outhwaite, Statistical Mechanics, 1975, 2, 188-255. 14. M. M. Martinez, L. B. Bhuiyan, C. W. Outhwaite, J. Chem. Soc., Faraday Trans. 1990, 86, 3383-3390. 15. C. W. Outhwaite, M. Molero, L. B. Bhuiyan, J. Chem. Soc., Faraday Trans. 1991, 87, 3227-3230. 16. C. W. Outhwaite, M. Molero, L. B. Bhuiyan, J. Chem. Soc., Faraday Trans. 1993, 89, 1315-1320. 17. C. W. Outhwaite, Condens. Matter Phys. 2004, 7, 719-733. 18. M. Molero, C. W. Outhwaite, L. B. Bhuiyan, J. Chem. Soc., Faraday Trans. 1992, 88, 1541-1547. 19. C. W. Outhwaite, L. B. Bhuiyan, V. Vlachy, B. Hribar-Lee, J. Chem. Eng. Data, 2010, 55, 4248-4254. 20. J.-P. Hansen, I. R. McDonald, Theory of Simple Liquids, Elsevier, London, 2006. 21. L. Blum, Mol. Phys. 1975, 30, 1529-1535. 22. E. Waisman, J. L. Lebowitz, J. Chem. Phys. 1972, 56, 3086-3093, 3093-3099. 23. C. Sanchez-Castro, L. Blum, J. Phys. Chem. 1989, 93, 7478-7482. 24. W. Ebeling, K. Scherwinsky, Z. Phys. Chem. (Leipzig), 1983, 264, 1-14. 25. L. Blum, J. S. H0ye, J. Phys. Chem. 1977, 81, 1311-1316. 26. R. Triolo, J. R. Griera, L. Blum, J. Phys. Chem. 1976, 80, 1858-1861. 27. R. Triolo, L. Blum, M. A. Floriano, J. Phys. Chem. 1978, 82, 1368-1370. 28. H. R. Corti, J. Phys. Chem. 1987, 91, 686-689. 29. A. Kovalenko, F. Hirata, J. Chem. Phys. 1999, 110, 10095-10112. 30. G. Schmeer, A. Maurer, Phys. Chem. Chem. Phys. 2010, 12, 2407-2417. 31. J. C. Rasaiah, in: E. E. Kunhardt, L. G. Christophorou, L. H. Luessen (Ed.): The Liquid State and its Electrical Properties, NATO ASI Series B, Vol. 193, New York, Plenum Press, 1988. 32. V. Vlachy, Annu. Rev. Phys. Chem. 1999, 502, 145-165. 33. D. A. McQuarrie, Statistical Mechanics, Carper Collins, New York, 1973. 34. V. Vlachy, B. Hribar-Lee, L. B. Bhuiyan, J. Chem. Eng. Data, 2010, 55, 1855-1859. 35. L. Verlet, D. Levesque, Physica, 1962, 28, 1124-1142. 36. J.-P. Hansen, P. Vieillefosse, Phys. Rev. Lett. 1976, 37, 391394. 37. L. Belloni, Chem. Phys. 1985, 99, 43-54. 38. L. B. Bhuiyan, V. Vlachy, C. W. Outhwaite, Int. Rev. Phys. Chem. 2002, 21, 1-36. 39. Z. Abbas, E. Ahlberg, S. Nordholm, J. Phys. Chem. B, 2009, 113, 5905-5916. 40. M. P. Allen, D. J. Tildesley, Computer Simulations of Liquids, Oxford University, New York, 1989. 41. J. P. Valleau, L. K. Cohen, J. Chem. Phys. 1980, 72, 59355941. 42. B. Jamnik, V. Vlachy, J. Am. Chem. Soc. 1993, 115, 660666. 43. M. Luksic, B. Hribar-Lee., O. Pizio, Mol. Phys. 2011, 109, 893-904. 44. J. Rasaiah, Chem. Phys. Lett. 1970, 7, 260-264. 45. F. H. Stillinger, R. Lovett, J. Chem. Phys. 1968, 48, 38583868. 46. B. Hribar, H. Krienke, Yu. V. Kalyuzhnyi, V. Vlachy, J. Mol. Liquids, 1997, 75, 74, 277-289. 47. G. A. Mansoori, N. F. Carnahan, K. E. Starling, Jr. T. W. Le-land, J. Chem. Phys. 1971, 54, 1523-1525. 48. B. Widom, J. Chem. Phys. 1963, 39, 2808-2812. Povzetek Za študij termodinamičnih lastnosti elektrolitov, asimetričnih po velikosti in naboju, smo uporabili približek HNC. Poleg presežne notranje energije in osmoznega koeficienta smo v širokem območju koncentracij z uporabo Hansen-Vieil-lefosse-Bellonijeve enačbe (HVB) izračunali tudi srednji koeficient aktivnosti elektrolita in se na ta način izognili dolgotrajnemu računanju preko Gibbs-Duhemove zveze. Izračune smo primerjali z novimi rezultati velekanonskih simulacij Monte Carlo. Rezultati kažejo, da se srednji koeficienti aktivnosti ter osmozni koeficienti elektrolitov, kot jih napove približek HNC, dobro ujemajo z rezultati računalniških simulacij. Točnost rezultatov je odvisna predvsem od zasedenosti prostora v raztopini. Izkaže se, da pri višji zasedenosti prostora, dobimo nekoliko boljše rezultate s približkom MSA.