the ROZ Theory and GCMC Simulations - American Chemical Society

J. Phys. Chem. B 2001, 105, 4727-4734

4727

Equilibrium Properties of a Model Electrolyte Adsorbed in Quenched Disordered Charged Media: the ROZ Theory and GCMC Simulations Barbara Hribar,† Vojko Vlachy,*,† and Orest Pizio‡,§ Faculty of Chemistry and Chemical Technology, UniVersity of Ljubljana, Ljubljana 1001, SloVenia, and Department for Modeling of Physico-Chemical Processes, Maria Curie-Sklodowska UniVersity, 200-31 Lublin, Poland ReceiVed: January 26, 2001; In Final Form: March 8, 2001

We present a theoretical study of the quenched-annealed system consisting of an annealed +1: -1 electrolyte and a disordered quenched matrix modeled as a charge and size asymmetric +Z: -1 electrolyte, with Z ) 4 or 10. The annealed electrolyte is in thermodynamic equilibrium with an external reservoir of electrolyte at concentration cb such that the chemical potentials of the confined and external electrolytes are equal. Both the matrix and the adsorbed electrolyte are modeled as charged hard spheres in a dielectric continuum. The replica Ornstein-Zernike (ROZ) equations, supplemented by the hypernetted chain (HNC) approximation, are solved for this system. The effects of charge and size asymmetry of the matrix ions, their concentration, the prequenched conditions on the pair distribution functions, and the thermodynamic properties of the annealed electrolyte are investigated. The concentration of adsorbed electrolyte c1 is calculated as a function of the concentration of external electrolyte cb. For low concentration of external electrolyte, the average concentrations of electrolyte inside the matrix c1 is higher than cb. On the contrary, for higher values of cb, the electrolyte is desorbed from the matrix, c1 < cb. The theoretical results are supplemented by data from a grand canonical Monte Carlo computer simulation. It is shown that the ROZ/HNC theory yields results which are generally in good agreement with the simulations.

1. Introduction Study of the equilibrium distribution of an electrolyte between the porous phase and the bulk solution, which is a classical subject of physical chemistry, is fundamental to many chemical and biological processes. In the past decade, the fluids adsorbed in porous materials, pictured as a disordered set of obstacles, became the subject of numerous theoretical studies (for review see, e.g., refs 1 and 2). These advances have been made possible by the progress in statistical-mechanical theories and computer simulation techniques for continuum systems with quenched disorder.3-18 Theoretical research has been focused on a better understanding of the experimental results,19-23 indicating that properties of the adsorbed fluids differ substantially from those of the bulk fluids. Modeling the adsorbent as a disordered set of obstacles inherently provides the effects of pore connectivity. Moreover, the microporosity of the adsorbent and its surface area can be characterized by the theory with some degree of confidence. On the other hand, a more realistic modeling of the microporous adsorbents and their characterization can be reached using reverse Monte Carlo simulations.24 Hopefully, future modeling of the adsorbent would reach the level necessary to describe systems of practical interest. Quenched-annealed systems composed of particles interacting via the hard sphere or Lennard-Jones interactions have been studied quite extensively, while systems with Coulombic interactions received little attention so far. The study of partly quenched ionic fluids has been initiated by Chakraborty, Bratko, * To whom correspondence should be addressed. † University of Ljubljana. ‡ Maria Curie-Sklodowska University. § On sabbatical leave from Instituto de Quimica de la UNAM, Mexico.

and Chandler25-27 and more recently extended in several contributions from our group.28-31 In particular, we have applied the replica Ornstein-Zernike (ROZ) theory to study the adsorption of model electrolytes in disordered matrix media consisting of charged or neutral obstacles. A limited comparison with computer simulations indicated29,31 that the ROZ formalism yields at least semiquantitatively correct results for these systems. Encouraged by the agreement between the ROZ theory and computer simulations,31 we propose to extend our calculations to systems where the matrix medium is considered to be a quenched charge and size asymmetric +Z: -1 electrolyte. Of course, more sophisticated modeling of the matrix would be desirable for several aspects. Note that joining together species differing in size and in charge (using, for example, the chemical association mechanism), one can design heteronuclear chains or ramificated matrix structures with a distribution of excluded volume and/or an electric field as desired in advance. However, it seems more logical to test the accuracy of the theory for simpler systems first, before more refined theoretical models are studied. Therefore, the computer simulation data of this study provide an important benchmark for theoretical methodology. In the first example studied here, the matrix parameters were chosen to correspond to a +4: -1 electrolyte mimicking a ThCl4solution.32 In the second example, a matrix with polycations having 10 positive charges and small monovalent anions was studied. The matrix medium was assumed to be equilibrated under conditions (temperature T0 and dielectric permittivity 0) which may be different from the conditions observed (T1, ). It was assumed that the matrix can be brought from the prequenching conditions to the observation temperature without structural relaxation. The annealed electrolyte, being distributed

10.1021/jp010346b CCC: $20.00 © 2001 American Chemical Society Published on Web 04/28/2001

4728 J. Phys. Chem. B, Vol. 105, No. 20, 2001

Hribar et al.

through the electroneutral matrix, was a charge symmetric electrolytesthe ionic radii were chosen to correspond to a KCl solution.32 The annealed electrolyte was assumed to be in equilibrium with an external reservoir (bulk) of the same electrolyte. The situation is in some aspects similar to the socalled Donnan equilibrium:33 in thermodynamic equilibrium, chemical potentials of mobile ionic species and the solvent within the microporous material phase must be equal to those in the bulk (external) electrolyte phase. In the study presented here, we are interested in the effects brought about by the presence of the matrix on the equilibrium properties of the +1: -1 electrolyte. Intuitively, one would expect for the mean activity coefficient of adsorbed electrolyte to be different from that in the equilibrium bulk electrolyte phase. The electrolyte adsorbed in the model matrix with charged obstacles was examined as a function of various parameters such as the concentration and charge asymmetry +Z: -1 of the matrix ions and the prequenching conditions. The major part of the computation applied to the ROZ equations, supplemented with the hypernetted chain closure (ROZ/HNC), which has been successfully applied to partly quenched systems in our previous studies.29-31 So far, a fairly detailed study of the structure of model electrolytes, as measured through the pair distribution functions between various species, has been presented. In the present paper, however, the emphasis was put on thermodynamic properties, more precisely, on the mean activity coefficient of the adsorbed electrolyte. In parallel, the grand canonical Monte Carlo (GCMC) method was used to study the equilibrium partitioning of the model electrolyte between the matrix and bulk solution. The outline of this paper is the following. After this introduction, we present the details of the model (section 2). The theoretical procedures used in this study are outlined in section 3. In section 4, the numerical results for structural and thermodynamic properties are presented. The conclusions are summarized in the final section (section 5).

respectively. The matrix ions are immersed in a dielectric continuum with the dielectric constant 0. In this case, the interaction potential reads

2. Model System

3. Theoretical Methods

The model system consists of two electroneutral subsystems. The first is a quenched ionic fluid, which is called the matrix, and the second is an annealed ionic fluid which thermally equilibrates in the presence of matrix species. The matrix does not respond to the presence of the annealed fluid. The notation used in this paper is as follows: the superscripts 0 and 1 correspond to the matrix and the annealed fluid species, respectively. The matrix, represented as an electroneutral system of positively and negatively charged hard spheres, is formed as described before:27 at a certain temperature T0, the electrolyte solution is subjected to a rapid quench. It is assumed that the structure of the solution does not change during this procedure; the structure of the matrix corresponds to an equilibrium state of an ionic fluid of concentration c0 at temperature T0. This assumption allows us to use integral equation methods to calculate the pair distribution functions for matrix ions. In this calculation, we choose to examine two different matrix models. In the first example, the matrix is formed from a +4: -1 electrolyte with parameters mimicking ThCl4 solution.32 The charges of the ions in this example are ez0+ and ez0-, where z0+ ) + 4 and z0- ) - 1. The diameters are chosen to be σ0+ ) 12.1 Å and σ0+ ) 3.62 Å. In the second example, we studied the matrix mimicking a +10: -1 electrolyte. The parameters for this case are z0+ ) + 10, z0- ) - 1, σ0+ ) 30 Å, and σ0+ ) 3.62 Å and do not correspond to any specific micellar system. The number concentration of matrix particles is F0+ and F0-,

Integral Equation Theory. The structure and thermodynamics for the model described above has been obtained by solving the replica Ornstein-Zernike equations. The equation for the matrix subsystem reads

U00 ij (r) )

{

r < (σ0+ + σ0-)/2 ∞ 2 0 0 e zi zj /0r r g (σ0+ + σ0-)/2

}

(1)

where i and j assume values + and -. The model for the annealed electrolyte is similar to that of the matrix; it corresponds to an electroneutral system of charged hard spheres with charges ez1+ and ez1- (z1+ ) 1, z1- ) -1) and diameters σ1+ and σ1-. The sizes of the ions were chosen to mimic KCl electrolyte solution: σ1+ ) 3.45 Å and σ1- ) 3.62 Å. The total ionic number concentration of annealed species is equal to F1 ) F1+ + F1-. The adsorbing ionic fluid is an electrolyte solution containing ions and solvent. The latter is, in analogy with the matrix subsystem, modeled as a dielectric continuum with constant . Adsorption may be studied at the temperature T1, which is different from the temperature T0 where the matrix was equilibrated. In analogy with previous studies,27-31 we define the quenching parameter Q ) e0T0/T1. The fluid-fluid U11 ij (r) and the fluid-matrix U10 ij (r) interaction potentials are defined as

U11 ij (r) )

{

}

(2)

U10 ij (r) )

{

}

(3)

and

r < (σ1i + σ1j )/2 ∞ , e2z1i z1j /r r g (σ1i + σ1j )/2

r < (σ1i + σ0j )/2 ∞ , e2z1i z0j /r r g (σ1i + σ0j )/2

respectively. The calculations presented in this paper apply to T1 ) 298 K and ) 78.54.

H00 - C00 ) C00 X F0H00

(4)

where F0 is a 2 × 2 diagonal matrix with diagonal elements F0+ and F0- while the symbol X denotes convolution. The correlation functions in eq 4 are 2 × 2 matrices with the elements 00 00 00 f00 - -(r), f+ +(r), and f+ -(r) ) f- +(r), where f stands for h or c. The ROZ equations for the fluid-matrix and fluid-fluid correlation are

H10 - C10 ) C10 X F0H00 + C11 X F1H10 - C12 X F1H10 H11 - C11 ) C10 X F0H01 + C11 X F1H11 - C12 X F1H21 H12 - C12 ) C10 X F0H01 + C11 X F1H12 + C12 X F1H11 2C12 X F1H21 (5) where F1 is a 2 × 2 diagonal matrix with diagonal elements F1+ and F1-. The correlation functions are 2 × 2 matrices with elements f+ +(r), f+ -(r), f- +(r), and f- -(r), each of functions describing annealed fluid species, H11and C11, consisting of connected and blocking part5-8. The connected parts of the functions correspond to correlations between particles within

Quenched Disordered Charged Media

J. Phys. Chem. B, Vol. 105, No. 20, 2001 4729

the same replica, whereas the blocking parts contain correlation between particles of different replicas. The particles belonging to different replicas are correlated merely to the presence of the matrix, and they do not interact directly. Functions H12 and C12 represent the blocking contribution to the pair correlation functions and to the direct correlation functions of the annealed fluid (fluid-fluid functions) H11 and C11, respectively. To solve the set of integral equations represented by eq 5, we need the so-called closure relations. In the present calculation, the ROZ equations were supplemented by the hypernetted chain (HNC) closure. This approximation was found to be very successful in description of bulk symmetric electrolytes34 and moderately successful for electrolyte solutions asymmetric in charge and size35. The HNC closure relations for the replica OZ system of equations, successfully used in several previous studies29-31, read

Cmn(r) ) exp[-βUmn(r) + γmn(r)] - 1 - γmn(r) C12(r) ) exp[γ12(r)] - 1 - γ12(r)

(6)

where γmn ) Hmn - Cmn, with the superscripts m and n assuming values of 0 and 1, and Umn are the matrices of interparticle interactions for different components. Because of the long-range Coulomb forces, eq 5 has to be renormalized before its numerical solution is possible. The renormalization procedure for partly quenched systems was explained in refs 29-31 and will not be repeated here. The ROZ equations, rewritten in a renormalized form, are presented below: 00 00 00 0 00 00 0 00 H00 (s) - C(s) ) C(s) X F (H(s) + q ) + Φ X F H(s) (7) 0 00 10 0 00 10 10 00 H10 (s) - C(s) ) C(s) X F (H(s) + q ) + Φ X F H(s) + 1 10 11 1 10 1 10 10 12 10 C11 (s) X F (H(s) + q ) + Φ X F H(s) - C(s) X F (H(s) + q ) 0 01 10 0 01 11 10 01 H11 (s) - C(s) ) C(s) X F (H(s) + q ) + Φ X F H(s) + 1 11 11 1 11 1 21 11 12 21 C11 (s)X F (H(s) + q ) + Φ X F H(s) - C(s) X F (H(s) + q ) (8) 0 01 10 0 01 12 10 01 H12 (s) - C(s) ) C(s) X F (H(s) + q ) + Φ X F H(s) + 1 12 11 1 12 1 11 12 12 C11 (s) X F (H(s) + q ) + Φ X F H(s) + C(s) X F H 1 21 2C12 (s) X F H

The set of ROZ/HNC equations was solved numerically by direct iteration over n ) 214 points and for the integration step ∆r ) 0.05 Å. This large number of integration points was needed because the interactions are very long range, and at the same time, pair distribution functions are rapidly varying near the multivalent matrix ions. Monte Carlo Simulations. The methodology of computer simulation of a quenched-annealed system has been explained in some detail in previous papers.25,27,31 First, the configurational properties of an adsorbed fluid are obtained as ensemble averages for a chosen matrix configuration, and then, the procedure is repeated for other matrix realizations in order to obtain the matrix average. In principle, this second average extends over all possible configuration states of the matrix. In practice, however, a small number of equilibrium matrix realizations is usually generated.16,25,31 An exception is the work of Bratko and Chakraborty,27 which applies to the infinitely dilute annealed fluid. Recent simulation studies of the phase transition behavior of a Lennard-Jones fluid in a disordered porous structure indicate16,36 that some details of the phase

diagram may be sensitive to the number of matrix realizations of the porous material. We believe that this problem does not apply to the system studied here; for values of the parameters used in this work, the partly quenched system in question is far from conditions where the results may be realization-sensitive.36 In the first stage of the simulation procedure we used the canonical Monte Carlo method to simulate the matrix, i.e., the +4: -1 or +10: -1 electrolyte at T0. The number of matrix particles used in the simulations varied from 250 to 550, depending on the matrix parameters. The matrix subsystem was equilibrated over 9 × 106 steps and the final configuration chosen as a first realization of the matrix. In the second stage the grand canonical Monte Carlo procedure was used to calculate the distribution of mobile ions between the porous and bulk phase. The grand canonical method has the advantage that by sampling at constant chemical potential, the relevant bulk electrolyte phase is defined unambiguously.37,38 In the GCMC simulations, the chemical potential of the solution is fixed by the external reservoir of electrolyte at temperature T1 and mean activity cbγ(,b. The concentration of the adsorbed electrolyte c1 follows directly from the simulation, while the mean activity coefficient γ(,1 is obtained from the equilibrium condition given by eq 9. To implement the GCMC simulation method for a partly quenched electrolyte system, one needs first to obtain the mean activity coefficients, γ(,b, of the bulk electrolyte phase at a chosen concentration cb. This quantity can be obtained by a separate GCMC simulation of a bulk electrolyte37 or by some other theoretical method. In the present work, we use the hypernetted-chain theory34,35 to calculate the mean activity coefficient of the model +1: -1 electrolyte as a function of the concentration cb. Note that by comparison with computer simulation, this theory was found to be very accurate for +1: -1 electrolytes.34,39 The initial number of annealed fluid ions to be distributed within the matrix was calculated according to the external electrolyte concentration cb, and it varied from 50 (at low concentration) to 600 (at high concentration). The Ewald summation method40 was applied to treat the Coulomb interactions adequately. The ions, distributed within the quenched asymmetric electrolyte representing the matrix, were first equilibrated over the 3-5 × 106 states. After the equilibration run, a production run consisting of 3-5 × 107 steps was needed to obtain reasonably good statistics. The simulation steps consisted of (i) attempts to either create or annihilate an ionic pair at randomly selected coordinates in the simulation cell and (ii) displacement of ions common to the canonical Monte Carlo simulation. The overall acceptance during simulation in the productive stage of all the systems studied here was from 30% to 50%. A new simulation run was started by changing the representative equilibrium distribution of the quenched fluid (no annealed fluid ions were present during this part of the procedure) by attempting an additional 9 × 106 configurations. The final configuration of this run was chosen as the next representative distribution of the matrix particles. Such a run provides a statistically independent matrix configuration, which, however, is characterized by the same statistical average in terms of the pair distribution function of matrix species. Next, the electrolyte ions were distributed in the “new” matrix, and with the matrix particles rigidly fixed in their positions, the grand canonical simulation for the annealed electrolyte was repeated. The average over “all matrix realizations” actually involves only two or three matrix configurations. As before,31 we did not observe


Hribar et al.

Figure 1. (a) Annealed fluid ion-ion pair distribution functions 11 g11 + -(r) and g+ +(r) for a model +1: -1 electrolyte solution in a +4: -1 matrix. (b) Fluid-matrix pair distribution functions g10 + +(r) and 0 g10 - +(r) for this case. The calculations apply to c+ ) 0.125 M, the quenching parameter Q ) 1.2, and c1 ) 0.0336 M. The symbols represent GCMC results, and the ROZ-HNC calculations are shown by lines.

significant changes of statistical averages as a function of the number of matrix representations. As a result of the computer simulation, the pair distribution functions and the average concentration of adsorbed electrolyte c1 were obtained. Note that for the charge symmetric electrolyte studied here, c1 ) F+1/NA, where NA is the Avogadro number, the average concentration of mobile cations c1+ (K+) and anions (Cl-) equals the average concentration of the annealed electrolyte c1. The mean activity coefficient of the adsorbed electrolyte γ(,1 is obtained from the thermodynamic relation

γ(,1 cb ) γ(,b c1

(9)

Finally, we would like to mention that only the total pair correlation functions were evaluated in this simulation. In other words, no attempt was made to calculate separately the blocking or connected contribution to the ion-ion distribution functions. Meroni et al.,10 in their study of a partly quenched hard sphere fluid, report certain practical difficulties in connection with the evaluation of the blocking contribution to the pair distribution function. For a partly quenched system with Coulomb interactions, this may require a nontrivial simulation procedure. We hope to address this issue in one of our future studies. 4. Results Pair Distribution Functions: We consider first results for the structure of the annealed +1: -1 electrolyte in the +4: -1 matrix. The pair distribution functions (pdfs), g11 ij (r), for this example are shown in Figures 1 and 2. Figure 1a shows the annealed electrolyte ion-ion pdfs at concentration c1 ) 0.0336 mol dm-3. The concentration of the matrix cations is

Figure 2. (a) Same pdfs as those for Figure 1a but for c1 ) 0.01 M and for two different values of the quenching parameter: Q ) 0.8 (solid lines) and 1.7 (dotted lines). (b) Fluid-matrix pair distribution functions; the parameters and notation are the same as those for Figure 2a.

c0+ ) 0.125 mol dm-3, and the quenching parameter Q ) 1.2. The fluid-matrix distribution functions are shown in Figure 1b. The agreement between the ROZ/HNC calculation (lines) and the GCMC results (symbols) shown in Figure 1a and 1b is very good. It is of interest to discuss these results in view of those obtained previously for a symmetric +1: -1 matrix.29 The most noticeable differences are in the shapes of the fluidmatrix pair distribution functions. In Figures 1b and 2b, we see 10 the so-called charge inversion; g10 + +(r) and g- +(r) functions cross at a certain distance from the multivalent matrix ion. In our view, this effect reflects the strong interaction between the multivalent (+4) matrix ions and the monovalent fluid counterions. In contrast to this, the annealed fluid-matrix pair distribution functions g10 ij (r), calculated for the +1: -1 electrolyte in the +1: -1 matrix, exhibit monotonic behavior (cf. Figure 3 of ref 29). The annealed fluid pdfs shown in Figures 1a and 2a agree qualitatively with those obtained for a symmetric +1: -1 matrix.29 Along with the concentrations of electrolyte and matrix ions, the quenching parameter Q represents another important parameter of the theory. The effects of this parameter on the structure of the adsorbed electrolyte is shown in Figure 2. These results apply to the highest matrix concentration studied in this paper, i.e., to c0+ ) 0.125 mol dm-3, while the annealed electrolyte concentration for this example is c1 ) 0.01 mol dm-3. In Figure 2a, we see the effect obtained by increasing the value of the quenching parameter from Q ) 0.8 (solid line) to 1.7 (dotted line) on the fluid ion-ion pdfs. Note that the distribution functions obtained at Q ) 1.7 exhibit nonmonotonic behavior. Altogether, the results presented in Figure 2a and 2b confirm previous findings27-31 that prequenching conditions play an important role in shaping the structure of these systems. The results for the +10: -1 matrix are shown in the next three figures. First, in Figure 3a, we show the annealed fluid ion-ion pdfs at c1 ) 0.01 mol dm-3 and for c0+ ) 0.009 11


Figure 3. (a) Annealed fluid ion-ion pair distribution functions 11 g11 + -(r) and g+ +(r) for a model +1: -1 electrolyte in a +10: -1 matrix. (b) Fluid-matrix pair distribution functions g10 + +(r) and 0 g10 (r). The calculations apply to c ) 0.009 11 M, c ) 0.01 M, and 1 -+ + Q ) 1.2. The symbols represent GCMC results, and the ROZ-HNC calculations are shown by lines.

Figure 4. (a) Annealed fluid ion-ion pair distribution functions 11 g11 + -(r) and g+ +(r) in the ROZ/HNC approximation for a model +1: -1 electrolyte in a +10: -1 matrix; c1 ) 0.01 M, c0+ ) 0.02 M. The values of the quenching parameter are Q ) 0.8 (solid lines) and 1.7 (dotted lines). (b) Fluid-matrix pair distribution functions; the parameters and notation are the same as those for Figure 4a.

mol dm-3. In Figure 3b, the ion-matrix distribution functions are presented for these values of the parameters; the value of


Figure 5. Same as that for Figure 4, but with c0+ ) 0.05 M and for two different fluid concentrations: c1 ) 0.001 M (solid lines), and 0.1 M (dotted lines); Q ) 1.2.

the quenching parameter is 1.2. Again, the ROZ/HNC pair distribution functions shown by lines are in good agreement with the computer simulations (symbols). In contrast to the previous case (+4: -1 example, Figure 1a), crossover behavior 11 of the functions for like and unlike ions g11 ++(r) and g-+(r) is observed in Figure 3a. This effect is caused by the strong and long-range interaction between the matrix and annealed fluid ions. Charge inversion, documented by crossing of the like and unlike annealed fluid-matrix pair distribution functions (see also Figures 1b and 2b), is shown in Figure 3b. More HNC/ROZ predictions for the -10: +1 matrix are presented in Figure 4, where the influence of the quenching parameter is shown; note the similarity with Figure 2. Again, we see that the value of the quenching parameter strongly affects the annealed fluid distribution functions. Finally, in Figure 5, we compare the structure of adsorbed electrolyte at two different concentrations c1. The concentration of matrix c0+ ) 0.05 mol dm-3(ηm ) 0.433) is relatively high in this case. The solid lines in Figure 5 apply to electrolyte concentration c1 ) 0.001 mol dm-3, and the dotted lines to c1 ) 0.1 mol dm-3. It is interesting that a 100-fold increase of the electrolyte concentration results in a relatively small change of the fluid structure in this concentration range. Thermodynamic Properties. In this subsection, we discuss the thermodynamic properties which can be calculated once the distribution functions are known. The excess internal energy of a charged fluid inside the matrix is given by

βEex/N1 )

1

11 x1i F1j ∫ dr g11 ∑ ∑ ij (r)Uij (r) + 2 i)+,- j)+,-

∑ ∑

i)+,- j)+,-

x1i F0j

∫ dr g10ij (r)U10ij (r)

(10)

where x1i ) F1i /F1 is the fraction of particles of species i in the annealed fluid. The excess internal energies for some of the


Hribar et al.

TABLE 1: Reduced Excess Internal Energies, -Eex/N1kBT1, for a +1: -1 Electrolyte Absorbed in a +4: -1 Matrix, as Predicted by the ROZ/HNC Theory and GCMC Smulations

TABLE 2: Reduced Excess Internal Energies, -Eex/N1kBT1, for a +1: -1 Electrolyte Absorbed in a +10: -1 Matrix, as Predicted by the ROZ/HNC Theory and GCMC Simulations

-Eex/N1kBT

-Eex/N1kBT

c0+ [M]

Q

c1 [M]

ROZ/HNC

MC

c0+ [M]

Q

c1 [M]

ROZ/HNC

MC

0.009 11

1.2

0.0001 0.001 0.01 0.1 0.0001 0.001 0.0355 0.0675 0.128 0.0001 0.001 0.01 0.1 0.0001 0.001 0.0336 0.0631 0.119 0.0001 0.001 0.01 0.1

0.555 0.513 0.432 0.410 0.929 0.882 0.688 0.650 0.626 0.875 0.853 0.795 0.711 1.187 1.137 0.917 0.862 0.812 1.497 1.441 1.248 0.935

0.702 ( 0.002 0.661 ( 0.002 0.633 ( 0.002 0.936 ( 0.002 0.866 ( 0.002 0.814 ( 0.002 -

0.009 11

1.2

0.02

0.8

0.0001 0.001 0.01 0.0313 0.0601 0.1 0.0001 0.001 0.01 0.1 0.0001 0.001 0.01 0.1 0.0001 0.001 0.01 0.1 0.0001 0.001 0.01 0.1

0.715 0.677 0.577 0.518 0.496 0.492 0.603 0.621 0.582 0.545 0.833 0.822 0.722 0.600 1.131 1.046 0.868 0.652 0.958 1.068 0.989 0.808

0.544 ( 0.004 0.507 ( 0.004 0.491 ( 0.004 -

0.05

0.125

0.8

1.2

0.125

1.7

examples studied in this paper are collected in Table 1 (+4: -1 matrix) and in Table 2 (+10: -1 matrix). In several cases, the ROZ calculations are supplemented by the GCMC results, and the agreement between the two types of calculations is reasonably good. Closer inspection of these results indicates that the agreement is generally better for higher concentrations of the annealed fluid. To calculate the equilibrium distribution of electrolyte between the porous phase and bulk electrolyte, we need the individual activity coefficients of all mobile species. For bulk electrolytes, there are several ways of calculating the activity coefficient, and one route, valid within the hypernetted-chain approximation, has been used on several occasions41,42

ln γi,b ) -

∑

j)+,-

Fjc(s)ij(0) + 0.5

∑ Fj ∫ dr [hij(hij - cij)]

j)+,-

(11) where c(s)(0) denotes the Fourier transform of the short-range part of the direct correlation function at k ) 0. We used this relation to calculate the mean activity coefficient γ(,b of the bulk +1: -1 electrolyte at the various values of cb, needed to implement the GCMC calculations. Recently eq 11 was rederived to be applicable within the replica formalism31

ln γi,1 ) 0.5

10 11 12 F0j c(s)ij (0) - ∑ F1j [c(s)ij (0) - c(s)ij (0)] + ∑ j)+,j)+,-

10 10 1 Fj0 ∫ dr h10 ∑ ij (hij - cij ) + 0.5 ∑ Fj ∑ ∫ dr j)+,j)+,j)+,11 11 12 12 12 [h11 ij (hij - cij ) - hij (hij - cij )] (12)

First, let us discuss the predictions of the ROZ/HNC theory for the equilibrium distribution of ions between the bulk electrolyte and charged matrix. These data are presented in Table 3 and also in Figures 6 and 7. The first two columns of Table 3 characterize the matrix. In the third column of this Table, the concentration of bulk electrolyte (cb) with the mean activity coefficient equal to γ(,b (column 4) is given. Note that these values of γ(,b (obtained by eq 11) were used as input to the

1.2

1.7

0.05

1.2

TABLE 3: Equilibrium Concentration c1 and the Mean Activity Coefficient γ(,1 of an Adsorbed +1: -1 Electrolyte as Obtained by the GCMC Simulationsa z0+

c0+ [M]

cb [M]

γ(,b

c1 [M]

γ(,1

4 4 4 4 4 4 4 10 10 10

0.125 0.125 0.125 0.125 0.05 0.05 0.05 0.009 11 0.009 11 0.009 11

0.031 25 0.0625 0.125 0.25 0.031 25 0.0625 0.125 0.01 0.031 25 0.0625

0.845 0.803 0.759 0.718 0.845 0.803 0.759 0.901 0.845 0.803

0.0336 0.0631 0.119 0.224 0.0355 0.0675 0.128 0.0106 0.0312 0.0600

0.786 0.796 0.800 0.802 0.743 0.743 0.739 0.850 0.846 0.836

a Q ) 1.2 in this calculation. The mean activity coefficients of bulk electrolyte γ(,b, used as input in the GCMC simulations, are calculated from the OZ/HNC theory34 and eq 11. Numerical uncertainties in c1 and in γ(,1 are estimated to be from 1% to 2%

GCMC simulations of the partly quenched system. In the fifth column, the concentration of the adsorbed electrolyte c1, in equilibrium with the bulk electrolyte, is given. Finally, in the last column, we present γ(,1, which is the mean activity coefficient of the adsorbed electrolyte. The results for c1 were obtained directly from the GCMC simulations, and γ(,1 was then calculated from eq 9. The GCMC results presented in Table 3 indicate that for low concentrations of the external electrolyte cb, the electrolyte is adsorbed in the matrix; i.e., c1 is higher than the corresponding equilibrium concentration of the bulk electrolyte. This effect reflects the attraction between the charges on the matrix particles and electrolyte ions. At higher external electrolyte concentrations cb, the electrolyte is “excluded” from the matrix; in other words, its equilibrium concentration in the matrix c1 is lower than in the bulk. This seem to be a combination of the electrostatic and volume-exclusion effects. Alternatively, the equilibrium distribution of ions between the bulk and porous phase can be presented as a ratio between the two activity coefficients Y ) γ(,1/γ(,b ) cb/c1 (see eq 9). The ROZ/HNC results for Y (Figures 6 and 7) were obtained as follows: (i) for a given concentration cb, we calculated γ(,b


Figure 6. (a) Ratio of the mean activity coefficients Y ) γ(,1/γ(,b for a model +1: -1 electrolyte in a +4: -1 matrix as a function of cb. The matrix concentrations are c0+ ) 0.05 M (dashed line) and c0+ ) 0.125 M (dotted line); Q ) 1.2. The GCMC results for c0+ ) 0.125 M are shown by circles, and the results for c0+ ) 0.05 M by crosses. (b) The ratio of the mean activity coefficients Y as a function of cb for Q ) 0.8 (solid line), 1.2 (dashed line), and 1.7 (dotted line); c0+ ) 0.125 M.

of the bulk electrolyte using eq 11. In the next step, (ii) we applied eq 9 and, by a trial and error procedure, determined the mean activity coefficient γ(,1 and concentration c1 of the adsorbed electrolyte. In this step, the mean activity γ(,1 was calculated from eq 12. Note that for Y < 1 there is an increased concentration of the adsorbing electrolyte in the matrix relative to the bulk, c1 > cb. For Y > 1, the situation is reversed, c1 < cb, and in this case, the electrolyte desorbs from the matrix. The comparison between the ROZ/HNC and GCMC results for Y is presented in the next two figures. Figure 6a shows the results for the +4: -1 matrix. The GCMC data shown by circles should be compared with the dotted line; these results apply to c0+ ) 0.125 M. In the same figure, the result shown by crosses corresponds to GCMC simulations at c0+ ) 0.05 M. In this case, the ROZ/HNC calculations are shown by a broken line. All the results are obtained for the quenching parameter Q ) 1.2. As we see, the general trends of the adsorption behavior are successfully captured by the solution of the ROZ/HNC approximation in conjunction with eq 12. In Figure 6b, the effect of the prequenching conditions on electrolyte adsorption is examined. The concentration of the adsorbed electrolyte increases by increasing values of the quenching parameter Q. The calculations presented in Figure 7 apply to a +10: -1 matrix. The GCMC simulations of partly quenched systems necessarily involve a large number of particle and are therefore time-consuming. For this reason, we present fewer results for this highly asymmetric case. A limited comparison between the computer simulations and the ROZ/HNC theory in this case again suggests that the theory is able to reproduce the main features of the computer results.


Figure 7. (a) Ratio of the mean activity coefficients Y ) γ(,1/γ(,b for a +1: -1 electrolyte in a +10: -1 matrix as a function of cb. The calculations apply to c0+ ) 0.009 11 M (solid line), 0.02 M (dashed line), and 0.05 M (dotted line); the quenching parameter is Q ) 1.2. The symbols represent the GCMC results for c0+ ) 0.00911 M. (b) The ratio of the mean activity coefficients Y as a function of cb for Q ) 0.8 (solid line), 1.2 (dashed line), and 1.7 (dotted line); c0+ ) 0.02 M.

5. Conclusions The equilibrium distribution of a model electrolyte between a microporous material and the bulk solution is of much interest for basic science as well as technology. The purpose of the present study was 2-fold: first, we were interested in the mean activity of the adsorbed electrolyte, γ(,1, relative to the value of the same quantity in the equilibrium bulk electrolyte. This important quantity was studied as a function of the matrix and the external electrolyte concentration and for different values of the quenching parameter. Alternatively, the results could be presented in the form of the adsorption isotherm. The second goal of this research was to apply the GCMC method to a partly quenched ionic system and to test critically the predictions of the ROZ/HNC theory. In the course of this study, two different types of matrices were examined. In the first case, the microporous material was depicted as a quenched charge and size asymmetric +4: -1 electrolyte. In the second case, the charge asymmetry was even higher, +10: -1. In both cases, the invading electrolyte was modeled as a +1: -1 electrolyte, with the parameters mimicking KCl solutions. As in previous studies,27-31 the solvent was treated as a dielectric continuum. The electrolyte adsorbed in the matrix was assumed to be in equilibrium with an external electrolyte of concentration cb and the mean activity cbγ(,b. Adsorption of the model electrolyte was studied by the ROZ/ HNC theory and GCMC computer simulations, varying the matrix and electrolyte concentrations and the quenching parameter. We found the concentration of the adsorbed electrolyte c1 to be higher than the equilibrium bulk electrolyte concentration cb if the concentrations of matrix particles and the bulk

4734 J. Phys. Chem. B, Vol. 105, No. 20, 2001 electrolyte were sufficiently low. For higher bulk electrolyte and/or matrix concentrations, the concentration of the adsorbed electrolyte are less than those in the bulk; in other words, c1 < cb. One important conclusion of this study is that the ROZ/ HNC theory, which is computationally much more efficient than the GCMC method, yields fair agreement with the machine calculations. As mentioned in the Introduction, the partly quenched system studied here resembles what is often called the Donnan equilibrium.33 The Donnan equilibrium is established by the distribution of a simple electrolyte between an aqueous electrolyte phase (bulk electrolyte) and a polyelectrolyte-electrolyte mixture in water, when the two phases are separated by a membrane not permeable to polyelectrolyte molecules. This system has been studied by numerous investigators, and modern theoretical techniques were recently applied to this classical problem.42-44 In these papers, the polyelectrolyte solution was modeled as highly asymmetric electrolyte such as, for example, the matrix fluid (cf. Section 2) in the present paper. We wish to stress here that the equilibrium systems examined before42-44 differ in one important aspect to the partly quenched system studied in our present work. The difference is that in the classical Donnan case, all the small ions are mobile and can travel through the semipermeable membrane. This is not the case in our system, where the matrix counterions neutralizing the macroions are quenched and form an integral part of the matrix. It is not surprising, therefore, that the results obtained in the present study differ qualitatively from those for the classical Donnan system. In the latter case, namely, a rejection (“negative adsorption”) of the electrolyte on the side of the polyelectrolyte solution is obtained at low electrolyte concentration (see, for example, ref 45), a result which is exactly the opposite to what was obtained for the partly quenched systems studied above. The present work permits several important extensions which will be the subject of our future studies. First, it is of considerable interest to explore electrolyte adsorption in matrices of more complicated structure than those merely of charged hard spheres differing in size and charge. Second, a problem of practical importance is how to influence electrolyte adsorption in templated porous materials by varying the template size and density.17 The methodological tools used in this work permit solution of these problems with reasonable accuracy. Acknowledgment. O.P. was supported by the CONACyT of Mexico under Research Project 25301-E. B.H. and V.V. acknowledge the support of the Slovene Ministry of Science and Technology. In addition, O.P. is grateful to the academic staff of the Department for Modeling of Physico-Chemical Processes at Maria Curie-Sklodowska University for hospitality during his sabbatical leave, in particular to Professors S. Sokolowski, A. Patrykiejew, and M. Borowko. References and Notes (1) Rosinberg, M. L. In New Approaches to Problems in Liquid State Theory; Caccamo, C., Hansen, J.-P. Stell, G., Eds.; Kluwer: Dordrecht, The Netherlands, 1999.

Hribar et al. (2) Pizio, O. Adsorption in Random Porous Media. In Computational Methods in Surface and Colloid Science; Borowko, M., Ed.; Marcel Dekker: New York, 2000. (3) Madden, W. G.; Glandt, E. D. J. Stat. Phys. 1988, 51, 537. (4) Madden, W. G. J. Chem. Phys. 1992, 96, 5422; 1995, 102, 5572. (5) Given, J. A.; Stell, G. J. Chem. Phys. 1992, 97, 4573. (6) Given, J. A.; Stell, G. Phys. A 1994, 209, 495. (7) Tatlipinar, H.; Pastore, G.; Tosi, M. P. Philos. Mag. Lett. 1993, 68, 357. (8) Lomba, E.; Given, J. A.; Stell, G.; Weis, J. J.; Levesque, D. Phys. ReV. E 1993, 48, 233. (9) Rosinberg, M. L.; Tarjus, G.; Stell, G. J. Chem. Phys. 1994, 100, 5172. (10) Meroni, A.; Levesque, D.; Weis, J. J. J. Chem. Phys. 1996, 105, 1101. (11) Vega, C.; Kaminsky, R. D.; Monson, P. A. J. Chem. Phys. 1993, 99, 3003. (12) Ford, D. M.; Glandt, E. D. J. Chem. Phys. 1994, 100, 2392. (13) Trokhymchuk, A. D.; Pizio, O.; Holovko, M. F.; Sokolowski, S. J. Phys. Chem. 1996, 100, 17004. (14) Kierlik, E.; Rosinberg, M. L.; Tarjus, G.; Monson, P. A. J. Chem. Phys. 1997, 106, 264. (15) Paschinger, E.; Kahl, G. Phys. ReV. E 2000, 61, 5330. (16) Sarkisov, L.; Monson, P. A. Phys. ReV. E 2000, 61, 7231. (17) Zhang, L.; Van Tassel, P. R. Mol. Phys. 2000, 98, 1521. (18) Fernaud, M.-J.; Lomba, E.; Lee, L. J. Chem. Phys. 2000, 112, 810. (19) Wong, A. P. Y.; Chan, M. H. W. Phys. ReV. Lett. 1990, 65, 2567. (20) Frisken, B. J.; Cannell, D. S. Phys. ReV. Lett. 1992, 69, 632. (21) Goh, M. C.; Goldburg, W. I.; Knobler, C. M. Phys. ReV. Lett. 1987, 58, 1008. (22) Maher, J. V.; Goldburg, W. I.; Pohl, D. W.; Lanz, M. Phys. ReV. Lett. 1984, 53, 60. (23) Wong, A. P. Y.; Kim, S. B.; Goldburg, W. I.; Chan, M. H. W. Phys. ReV. Lett. 1993, 70, 954. (24) Thomson, K. T.; Gubbins, K. E. Langmuir 2000, 16, 5761. (25) Chakraborty, A. K.; Bratko, D.; Chandler, D. J. Chem. Phys. 1994, 100, 1528. (26) Bratko, D.; Chakraborty, A. K. Phys. ReV. E 1995, 51, 5805. (27) Bratko, D.; Chakraborty, A. K. J. Chem. Phys. 1996, 104, 7700. (28) Hribar, B.; Pizio, O.; Trokhymchuk, A.; Vlachy, V. J. Chem. Phys. 1997, 107, 6335. (29) Hribar, B.; Pizio, O.; Trokhymchuk, A.; Vlachy, V. J. Chem. Phys. 1998, 109, 2480. (30) Hribar, B.; Vlachy, V.; Pizio, O.; Trokhymchuk, A. J. Phys. Chem. B 1999, 103, 5361. (31) Hribar, B.; Vlachy, V.; Pizio, O. J. Phys. Chem. B 2000, 104, 4479. (32) Simonin, J. P.; Bernard, O.; Blum, L. J. Phys. Chem. B 1998, 102, 4411. (33) Guggenheim, E. A. Thermodynamics, 5th ed; North-Holland: Amsterdam, 1967; pp 305-307. (34) Rasaiah, J. C. In The Liquid state and its Electrical Properties; Kunhardt, E. E., Christophorou, L. G., Luessen, H. L., Eds.; N. A. T. O. ASI Series B 193; Plenum: New York, 1988. (35) Vlachy, V. Annu. ReV. Phys. Chem. 1999, 50, 145. (36) Page, K. S.; Monson, P. A. Phys. ReV. E 1996, 54, R29. (37) Torrie, G. M.; Valleau, J. P. J. Chem. Phys. 1980, 73, 5807. (38) Valleau, J. P.; Cohen, L. K. J. Chem. Phys. 1980, 72, 5935. (39) Vlachy, V.; Ichiye, T.; Haymet, A. D. J. J. Am. Chem. Soc. 1991, 113, 1077. (40) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Oxford University: New York, 1989. (41) Belloni, L. Chem. Phys. 1985, 99, 43. (42) Vlachy, V.; Prausnitz, J. M. J. Phys. Chem. 1992, 96, 6465. (43) Svensson, B.; Akesson, T.; Woodward, C. E. J. Chem. Phys. 1991, 95, 2717. (44) Zhou, Y.; Stell, G. J. Chem. Phys. 1988, 89, 7010. (45) Vlachy, V.; Haymet, A. D. J. Aust. J. Chem. 1990, 43, 1961.

the ROZ Theory and GCMC Simulations - American Chemical Society

Recommend Documents