292 Acta Chim. Slov. 2006, 53, 292–305 Scientific Paper Modelling Electrolyte Adsorption in Nanoporous Material Miha Lukšič a, Vojko Vlachy a*, and Orest Pizio b aFaculty of Chemistry and Chemical Technology,University of Ljubljana, Aškerčeva 5, 1000 Ljubljana, Slovenia bInstituto de Química de la UNAM, Circuito Exterior, Coyoacán 04510, México D.F. Received 11/5/2006 Dedicated to the memory of Prof. Dr. Davorin Dolar Abstract The structural and thermodynamic properties of a model electrolyte solution confined in disordered matrices with charged obstacles were studied by means of the grand canonical Monte Carlo simulation. A disordered nanoporous medium was modelled as i) an equi librium distribution of ions in a +1:-1 primitive model electrolyte; ii) a system of dipolar hard spheres; iii) a collection of chainlike molecules with alternating charge on the beads (polyampholyte); and iv) as a system of charged chainlike molecules (oligoelectrolyte) with the pertaining counterions. The confined electrolyte was assumed to be in thermodynamic equilibrium with the obstacles forming the nanoporous matrix and an external electrolyte of the same chemical composition. The solvent in all these cases was treated as a dielectric continuum. In the present study we were interested in effects of the distribution of charged obstacles on the mean activity coefficient of the confined electrolyte. The computer simulations were performed for a set of values of the model parameters such as the concentration of matrix ions and of the annealed electrolyte, pre-quenching conditions and the conditions of observation. The results confirmed our previous findings that the properties of an annealed electrolyte depend strongly on the conditions of observation (temperature and dielectric constant of solvent), as well as on the concentrations of all components. The effect of the matrix-charge distribution, investigated in this work, was found to be significant and more important for higher Coulomb couplings. Key words: nanoporous material, adsorption, electrolyte, Monte Carlo simulation 1. Introduction Understanding mechanisms governing electrolyte distribution between bulk solutions and nanoporous materials is of substantial interest for science as well as for technology. The latter interest stems from implications for separation and purification processes, heterogeneous catalysis, chemical sensing, and related areas (see, for example, ref. 1, 2, 3). In recent years such systems became the subject of numerous theoretical studies; particular attention has been paid to adsorbents containing charges, acknowledging the role that the Coulomb interaction plays in separation processes. The system of interest is heterogeneous on a molecular scale. One component (obstacles) is fixed in space (quenched component) and represents the nanoporous material, while the other component is fully mobile and thermally equilibrated within the matrix (annealed component). Such systems differ from regular two component mixtures; the statistical-mechanical average that is needed to obtain the thermodynamic properties of a confined fluid becomes a double ensemble average. First, at a given configuration of the matrix, the usual averaging over the annealed degrees of freedom takes place, while the second average is performed for all the possible values of the quenched variables. A decade ago we began a systematic investigation of partly-quenched ionic systems addressing the problem of screening of Coulombic forces by charged obstacles.4-11 In these papers the replica Ornstein–Zernike (ROZ) theory, developed by Madden, Glandt, Given, Stell and several others12-18 was generalized to study the adsorption of a model ionic fluid in a disordered medium with ionic and neutral obstacles (for a review, see ref. 19, 20). In addition to the replica theory, the grand canonical Monte Carlo method was utilized to examine the distribution of fluid between the porous and bulk phases.8,9,21-24 Computer simulations are perhaps the most useful tool in studying adsorption phenomena. The open ensemble Monte Carlo method is particularly well suited for this purpose for it unambiguously defines the equilibrium bulk phase. This approach, pioneered by Valleau and coworkers,25,26 has been successfully applied to partly-quenched systems.8-11 In several previous Lukšič et al. Modelling Electrolyte Adsorption in Nanoporous Materials Acta Chim. Slov. 2006, 53, 292–305 293 contributions from our group8-11 we studied electrolyte rejection from a charged matrix modelled by one possible equilibrium distribution of ions in the primitive model +1:-1 electrolyte. The matrix was assumed to be formed as follows: at a certain temperature T0 the electrolyte solution was subjected to a rapid quench. In addition, it was assumed that the spatial distribution of matrix ions (obstacles) did not change during the quench, and was therefore completely determined by the pre-quenching conditions given by T0 and the dielectric constant of the solvent (?) at that temperature. Within such a quenched electroneutral system of ionic obstacles, another electrolyte was distributed and allowed to thermally equilibrate at conditions (temperature, dielectric constant) that may have been different from those given by T0 and ?. Using computer simulations and the replica Ornstein–Zernike theory we were able to determine adsorption isotherms under a broad range of laboratory conditions. Note that the porous phase was always modelled as an electroneutral subsystem; in such situations the rejection of electrolyte is expected, and an unequal distribution of the fixed positive and negative charges in the matrix may be the only cause of the electrolyte intake. In the studies presented so far8-11 little attention has been paid to the modelling of the matrix phase. It seems reasonable to assume, however, that the spatial distribution of charged obstacles, in addition to other parameters, strongly influences the equilibrium distribution of electrolyte between the two phases. In the present study we utilized the computer simulation technique to make quantitative predictions about this effect. We used the grand canonical Monte Carlo method to examine the distribution of a charge symmetric +1:-1 electrolyte between the nanoporous phase containing an electroneutral distribution of ionic obstacles and bulk electrolyte solution. The parameters for the bulk electrolyte were chosen to mimic solutions of lithium chloride (LiCl) in solvents of various dielectric constant. We examined several different types of obstacle distribution, here called the matrix. In all the cases studied here these matrices were electroneutral. In one of the examples, the adsorbent was represented as a collection of fused charged hard spheres mimicking a dipolar fluid. In another case the matrix was modelled by an equilibrium distribution of chain-like oligomers modelling polyampholytes. In this example each oligomer with the number of monomer units N = 16 (or N = 32) was electroneutral having an equal number of negatively and positively charged beads. In yet another variant studied here, the matrix was assumed to be prepared from an equilibrium solution of charged oligomers (all the beads of the oligomer carrying a unit positive charge) with the respective counterions being distributed in solution. Both counterions and oligoions were assumed to be frozen in space, modelling a particular equilibrium charge distribution which can form in such a system. These results were compared with our previous (and new) calculations where the matrix material was represented as an electroneutral set of quenched charged hard spheres mimicking an equilibrium distribution of ions in the electrolyte solution. In summary, the purpose of this study was not to model a particular adsorbing material but rather to examine the effects caused by various distributions of charged obstacles. 2. Model system Our model system consists of two subsystems: the first is a quenched fluid with charged obstacles, which we call the matrix, while the second is a mobile (annealed) electrolyte solution. The latter equilibrates in the presence of matrix species. The annealed fluid is in thermodynamic equilibrium with an external (bulk) electrolyte of the same chemical composition. It is important to stress that within this model the matrix does not respond to the presence of the annealed fluid; this is how our system differs from regular mixtures. The notation used throughout the paper is that indices 0 and 1 depict matrix and annealed fluid species, respectively. a.) •0 b.) -\ Figure 1. (a) Dipolar obstacles (fused hard spheres of equal size bearing positive and negative charge); (b) chain-like obstacles with beads carrying alternating charges (N = 8), and (c) chain-like matrix particles with beads carrying positive charge neutralized by negative counterions (N = 8). In the present contribution we study the effects of different representations of the matrix subsystem. Firstly, i) the matrix is represented by positively and Lukšič et al. Modelling Electrolyte Adsorption in Nanoporous Materials 294 Acta Chim. Slov. 2006, 53, 292–305 o negatively charged hard spheres with charges ez+ and ez (z+ =1, z° = -l), such as to satisfy the electroneutrality condition. This is the so-called ”primitive model” of an electrolyte solution. Secondly, ii) the particles representing the matrix are fused hard spheres with positive and negative charges at their centres (see Figure 1a). The net charge of such a ”molecule” is zero; this is the model of a dipolar fluid. iii) The obstacles are chain-like molecules with beads carrying charges alternating in sign (polyampholyte), see Fig. 1b. In the final example studied here, iv) the obstacles are again presented as chain-like molecules, but with all the beads carrying charges of the same sign (Fig. 1c). In this case an equivalent number of quenched counterions is present to ensure the electroneutrality of the matrix subsystem. In all the cases treated we assume the matrix be formed as follows: at a certain temperature T0 and dielectric constant e, the respective fluid is subjected to a rapid quench. The distribution of particles in the matrix therefore corresponds to an equilibrium state of the fluid of concentration c0 at temperature T0. While the matrix subsystem is electroneutral in all the selected examples, the distribution of positive and negative charges is substantially different. In this study we are interested how the unequal distribution of charged obstacles affects the adsorbing power of the neutral adsorbent. In all the examples the size of the matrix particles are set equal to a°+ =a° =4.25 A and the interaction pair potential between the matrix particles (beads) is defined by: m oo Aiveoer r < ( ( [a\ + <7°)/2 (3) Most of the calculations in this paper apply to matrices formed under conditions characterized by the Bjerrum length /lB0=7.14 A \XB = eQ/4xL0dcBT). The structural and thermodynamic properties of the annealed electrolyte are examined for /lB0=7.14 A and 14.28 A at two different concentrations of matrix particles. The different values of /lBl used in the calculation may be ascribed to solvents having different dielectric permittivity Lx at temperature T1. 3. Grand Canonical Monte Carlo Method In studying adsorption processes it is natural to use the open ensemble method. In this method the configurational states are sampled at constant chemical potential µi and temperature T, while the concentration of ions inside the porous phase fluctuates. The simulation consists of two steps, the first one is canonical - the number of ions is fixed and a randomly chosen ion is moved into a new random position in the simulation cell. The attempted move is accepted with a probability fij equal to: fa = min[l,exp(—/3(t/j — Ui))] (4) where Ui is the configurational energy of state i (old position of the ion), Uj is the configurational energy of state j (attempted new position of the ion), and as usual p = \l(kBT) . In the grand canonical step the probabilities for insertion (and for deletion) of an electroneutral pair of singly charged ions to (from) the matrix are:25,26,27 fij = min[l, Yij] for adition fji = min[l, 11 Yij] for adition (5) After calculating the energies Ui and Uj , the transition probability fij from the state i with number of cations N* (and anions N~) to the state j, where jyt = N* +1 (and NJ = N~ +1), is: Ytj exp(B — /3(Uj — Ui)) (6) Uh I oo r< (aj +oj)/2 I S^7 r^(°i+^)/2 and (2) Parameter B = J3Juex+ \n(N+N_) in eq. 6 is related to the excess chemical potential (//ex = ju-ju ' ea ) of the bulk electrolyte with number concentration of cations p+ and anions p– and needs to be known in advance. I Lukšič et al. Modelling Electrolyte Adsorption in Nanoporous Materials Acta Chim. Slov. 2006, 53, 292–305 295 The computer simulation of this system involves in principle a double average. The first is the usual configurational averaging over the annealed degrees of freedom for a selected (equilibrium) distribution of obstacles. Next, it is necessary to perform averaging over the all possible distributions of matrix particles. In practice it turned out21-24,28-30 that only a small number of equilibrium matrix realizations is sufficient (see also, ref. 8-11) to obtain valid averages. We can ascribe this finding to the fact that for a sufficiently large matrix subsystem, the annealed particle visits all portions of the phase space relevant for the statistical average.22,28 To avoid effects due to the finite number of particles included in a simulation cell, the minimum image boundary conditions are used. We found in previous simulations, carried out at high couplings, that minimum image boundary conditions are perfectly acceptable under these conditions.11 4. Results 4.1 Numerical details Before presenting the results we wish to discuss some details of the Monte Carlo calculations. The first step was always the simulation of the matrix. In the first example i) the matrix was the +1:-1 electrolyte considered in the primitive model approximation. We applied the canonical Monte Carlo method to simulate an electroneutral system of equally-sized charged hard spheres (modelling a size symmetric +1:-1 electrolyte) at concentrations of c0=0.5 mol/L (and 1.0 mol/L), and for ?B,0=7.14 A. This value of ?B,0 corresponds to aqueous solutions at 298 K. The number of particles (sum of cations and anions) used in this calculation was 1000. The distribution of ions in the adsorbent was one of the possible equilibrium distributions (it was essentially the final one, taken after 10 × 106 attempted configurations) of the ions in the electrolyte solution formed under conditions T0 and ?. In the second example ii) the matrix was represented as a system of dipolar molecules (Fig. 1a), modelled as two fused hard spheres carrying positive and negative charges at their centres. The number of dipoles used in the simulation was in all cases equal to 500 and the machine calculations were performed for the monomer concentrations c0=0.5 mol/L (and 1.0 mol/L), and for two values of the parameter ?B,0 (?B,0=7.14 A and 14.28 A). In examples iii) and iv) 512 cations and 512 anions were arranged so as to form chainlike structures containing N=16 or N=32 monomeric units per chain. In example iii) the charges on monomer units in the chain alternated, while in case iv) all the beads on the chain carried positive charge. The anions (counterions) - their number was equal to the number of beads (monomer units) to satisfy the electroneutrality condition - were in the latter case distributed in the solution. Conditions of preparation corresponded to AB 0=7.14 A and c0=0.5 mol/L and 1.0 mol/L; in case iv) only to c0=1.0 mol/L. For iii) the simulation cell contained 32 chains with N=32 and 64 chains for N=16. On the other hand, in example iv) we treated simultaneously 16 oligomer chains with degree of polymerization N=32 and 512 related counterions in the basic Monte Carlo cell. In both cases, as also for ii) above, c0 denotes the concentration of monomer units (beads) composing the oligomer (or a dipole). In the examples ii-iv) the MOLSYM31 computer program was used to generate the equilibrium distributions. In each of these studies we needed around (20-50) × 103 passes (trial moves per particle) to equilibrate the Monte Carlo system. The simulations of charged oligomers with N=16 and N=32 were more demanding in this respect (50 × 103 passes). We used three types of moves for the chains in the simulation box. Chains were ”moved” by using a) translation, where the whole chain was translated by a particular random vector, b) rotation of the shorter end of the chain, and c) bend stretching. In the next step, the grand canonical Monte Carlo simulation was used to study the distribution of simple electrolyte between the bulk fluid and the matrix subsystem. The number of ions mimicking LiCl distributed within the matrix varied from ~ 60 at low to ~ 2000 at high concentrations. The annealed ions were first equilibrated over (1-3) × 106 states. After the equilibrium run, the production run of (6-50) × 107 attempted configurations was carried out to obtain canonical averages. 4.2 Characterization of the matrix As already mentioned, the matrix was represented by an equilibrium distribution of obstacles obtained by a sudden quench at T0 and s. The diameters of the individual (matrix) ions or beads were in all the examples studied here. In the first case, i) the matrix was represented by a restrictive primitive model +1:-1 electrolyte at concentrations c0=0.5 and 1.0 mol/ L (XB 0=7.14 A ). The distribution of charged obstacles could be characterized by the relevant pair distribution functions (PDFS), g®°(r/a), which for a solution with concentration c0=1.0 mol/L are shown in Figure 2. In ii) we picture the obstacles as fused hard spheres with positive and negative charges at their centres (see Figure 1a). The +,+ and +,- pair distribution functions for the concentration of monomers equal to c0=1.0 mol/L and for A,B0=7.14 A are shown in Figure 3. As seen from this figure, the solution exhibited very little structure beyond the first peak for this concentration. Lukšič et al. Modelling Electrolyte Adsorption in Nanoporous Materials 296 Acta Chim. Slov. 2006, 53, 292–305 r/o Figure 2. The ion-ion distribution functions gf(r/a) for the matrix formed from a +1:-1 primitive model electrolyte (case i). Empty symbols show PDF for unlike ions, full symbols present like-ion PDF. c0 = 1.0 mol/L and AB0= 7.14A. For the chain-like molecules (N=32 in our case) it is also important to have some informa tion about the global shape of the quenched matrix molecule. This information was provided from calculation of the end-to-end distances and radii of gyration during the Monte Carlo sim ulation of the matrix subsystem. For iii), where the beads on the electroneutral chain-like molecule carry alternating positive and negative charges (see Fig. 1b), the end-to-end distance (i?JL)1/2 was 39.16 A for c0=0.5 mol/L and 38.36 A for 1.0 mol/L solution. The corre sponding values of the radius of gyration {R-g)x 2 were 15.44 A and 15.19 A. The PDFS for the unlike charge (empty symbols; value at contact distance is not shovn and extends over 45.6) and like-charge beads (full symbols) of the polyampholyte matrix (c0=1 mol/L, A,B0=7.14 A and N=32) are shown on Figure 4. As for the case ii) little correlation between the beads of the chains is seen after the second peak for this conditions. 0 1 2 3 4 5 6 7 8 9 r/? Figure 4. The +,+ and +,- charge distribution functions gi0j0 (r /?) for the matrix formed from polyampholyte chains (case iii). Empty symbols show results for unlike charge PDF (value at contact distance extends over 45.6), full symbols present like-ion PDF. Concentration of monomers c0= 1.0 mol/L and ?B,0 = 7.14 A. In the last example iv), the chain-like molecules carry charge of the same sign (Fig. 1c) and the equivalent number of monovalent counterions is distributed in the matrix subsystem to make it electroneutral. The end-to-end distances (-ftfL)1/2in this case were 56.54 A for c0=0.5 mol/L and 51.95A for c0=1.0 mol/L solution, while {Rg)X 2values for these two con centrations were 20.89 A and 19.68 A respectively. These data show that charged chains iv) are considerably more extended than the electroneutral polyampholyte chains iii). Finally, in Fig ure 5 we present the oligoion-counterion and counterion-counterion pair distribution functions for this case. The graphs (empty symbols) indicate a relatively strong association of counteri ons with oligoions as normally found in linear polyelectrolyte solutions (see, for example, ref. 32, 33). r/? Figure 3. The +,+ and +,- charge distribution functions gi0j0(r /?) for the matrix formed from dipolar obstacles (case ii). Empty symbols show the PDF for unlike charge beads, full symbols present like-charge PDF at concentration of monomers c0= 1.0 mol/L and for ?B,0=7.14 A. r/o Figure 5. The counterion-counterions (full symbols) and oligoion-counterion (empty symbols) pair distribution functions gi0j0 (r /?) for the matrix formed from charged oligomers and respective number of counterions (case iv) at concentration of monomers c0= 1.0 mol/L and for ?B,0 = 7.14 A. 3 3 2.5 2.5 1.5 1.5 0.5 0.5 6 1.2 2.5 0.8 1.5 0.6 0.4 0.5 0.2 3 4 9 3 Lukšič et al. Modelling Electrolyte Adsorption in Nanoporous Materials Acta Chim. Slov. 2006, 53, 292–305 297 In general we may classify the model adsorbents studied here into two types. The first class is represented by dipolar fluid ii) and polyampholyte iii), where the individual particles (obsta cles) composing the matrix are neutral. The second class is constituted by +1:-1 electrolyte i) and polyelectrolyte solution iv), where the obstacles carry some net charge. The separation of the positively and negatively charged obstacles, possible in the latter two cases, has a decisive influence on the adsorption behaviour of the adsorbent. 4.3 Pair-distribution functions We start our discussion with the selected pair distribution functions between the annealed ions themselves, gy (rl a) and for the annealed ion-matrix particle (an ion or a bead) correlation pairs, gy°(rl a). For the annealed electrolyte confined in different matrices subscripts i and j denote Li (+) and CI (—) ions. Because of the small difference in the shapes of the PDFS for matrices with concentrations c0=0.5 and 1.0 mol/L, generally only the case for c0=1.0 mol/L is shown. Similarly, for matrices formed from chain-like molecules the difference in PDFS for N=16 and N=32 was small, and the majority of the results shown here apply to N=32. By the same reasoning, only the results where the matrices were prepared under conditions given by the Bjerrum length AB 0 = 7.14 A are presented. In all figures we use the reduced distance units r/a where a =4.25 A. i) Ionic matrices formed from the size symmetric +1:-1 electrolyte. First, in Figure 6a we show the pair distribution functions for annealed cations and anions modelling Li and CI ions in the 1.0 mol/L electroneutral ionic matrix prepared at AB 0=AB t= 7.14 A. The upper curves (with the values above unity) are for Li-Cl pair, and the lower ones (the values below unity) are for the Li-Li pair. Concentrations of the annealed electrolyte were c=0.1070 mol/L (empty symbols), and c^O.4996 mol/L (full symbols) in this example. In Figure 6b we show the same distributions at stronger electrostatic coupling; AB1 =14.28 A and c =0.1520 mol/L, with all other conditions unchanged (AB0=7.14 A). The ”overcharging effect”, reflected in the crossing of the Li-Cl and Li-Li distributions, noticed before in other studies,5,6,22 is clearly seen in this figure. Fluid-matrix pair distribution functions gx\ (CI ion-positive matrix charge) for c1=0.4996 mol/L, XB1= 7.14 A (empty symbols) and c=0.5295 mol/L, AB j = 14.28 A (full symbols), are displayed in panel (c) of the same figure. In summary, the results displayed in Figure 6 are similar to those published before; 5"8 they are shown here for the sake of completeness and to allow comparison with the distribution functions calculated for other choices of matrices, i.e. for examples ii-iv). 0 1 2 3 4 5 6 7 8 9 r/a 0 1 2 3 4 5 6 7 8 9 r/a 0 1 2 3 4 5 6 7 8 9 r/a Figure 6. (a) Annealed fluid ion-ion pair distribution functions g\\(r/a) and g\l_(r/a) for a model LiCl in 1.0 mol/L electroneutral ionic matrix prepared at AB0=7.14 A, AB1 = 7.14 A; cx =0.1070 mol/L (empty symbols), c1=0.4996 mol/L (full symbols). (b) As for case a) but for AB1=14.28 A and c1 =0.1520 mol/L. (c) Fluid-matrix pair distribution functions gx\{rlo) for this case: AB1 =7.14 A, c1 =0.4996 mol/L (empty symbols); AB1 =14.28 A cx= 0.5295 mol/L (full symbols). a = 4.25 A ii ) Matrix formed from dipolar particles as shown in Fig. 1a. Figures 7a-c show PDFS for adsorption of a model LiCl in the c0= 1.0 mol/L (note that here c0 is the concentration of monomer beads) dipolar matrix prepared at ?B,0 = 7.14 A. The actual conditions 3 2.5 1.5 0.5 3 2.5 1.5 0.5 3 2.5 1.5 0.5 Lukšič et al. Modelling Electrolyte Adsorption in Nanoporous Materials 298 Acta Chim. Slov. 2006, 53, 292–305 0 1 2 3 4 5 6 7 8 9 r/a 0 1 2 3 4 5 6 7 8 9 r/a 0 1 2 3 4 5 6 7 8 9 r/a 1 2 3 4 5 6 7 8 9 r/a r/a Figure 7. (a) Annealed fluid ion-ion pair distribution functions gl\(r/a) for a model LiCl in 1.0 mol/L dipolar matrix prepared at XB0 =7.14 A. Symbols denote: (•) AB1= 7.14 A, c=0.6085 mol/L; (o) XB1 = 14.28 A, c1=0.6083 mol/L; (?)AB1 = 28.56 A, cj= 0.5981 mol/L. (b) The same as in case a) for g"_(rl a). (c) Fluid-matrix pair distribution functions gl°+(r< c_ > ±^m (7) where < c+> and < c > are the average concentrations of cations and anions in the nanoporous phase (< c+ >=< c->=< c1>) and J+bulk is the mean activity coefficient of the bulk electrolyte. To obtain this quantity (T±bulk) for a given electrolyte concentration cbulk one has to run a separate grand canonical Monte Carlo simulation of the bulk electrolyte, or to calculate the quantity by some other theoretical method. In the present study, we used the hypernetted-chain theory to calculate the mean activity coefficient of the model LiCl electrolyte as a function of its concentration cbulk (see, for example, ref. 35). Eq. 7 defines the adsorption isotherm. These results are shown in Fig. 12, where the average concentration of annealed electrolyte, , is presented as a function of the mean activity, a\ = y±1cx.The isotherms for adsorption in different matrix representations are shown as follows: empty symbols denote adsorption in a 0.5 mol/L matrix, while full symbols denote adsorption at c0=1.0 mol/L. In all the cases shown here the matrix was prepared at ?B,0 =7.14 A. Further, squares represent adsorption at stronger Coulomb coupling, ?B,1 =14.28 A, while circles denote the results for ?B,1=7.14 A. As can be seen from Figures 12a-d, the amount of electrolyte adsorbed increases on increasing the mean activity of the adsorbing electrolyte, and it is lower for more concentrated matrices. This is not difficult to understand: for concentrated systems more space is occupied by the matrix particles and the exclusion volume effect, forcing electrolyte out of the matrix, plays a more important role. The difference between the amount of electrolyte adsorbed in 0.5 mol/L and 1.0 mol/L matrices is smaller at lower electrolyte activities. In addition, the adsorption is stronger for higher Bjerrum lengths. Electrostatic interactions between fluid ions, as well as between fluid and the matrix charges, in the case of ?B,1=14.28 A are twice as strong as for ?B,1 =7.14 A. It can be seen from Figures 12a-d that doubling of the Bjerrum length, and concomitant increase of the Coulomb interaction between ions of the system, has a more pronounced effect on the amount of electrolyte 0.4 0.1 0.2 0.3 0.4 0.5 0.6 O 0.4 a± (mol/L) 0.1 0.2 0.3 0.4 0.5 O.e a± (mol/L) 0.4 0.2 0.3 0.4 a± (mol/L) O 0.4 0.2 0.3 0.4 a± (mol/L) Figure 12. Adsorption isotherms (lines with symbols) for a LiCl model electrolyte confined in a) disordered electroneutral ionic matrix, b) dipolar matrix, c) chain-like matrix with beads carrying alternating charges (N=32) and d) chain-like matrix with beads carrying positive charge neutralized by negative counterions (N=32). From top to bottom: c 0= 0.5 mol/L, AB 1= 14.28 A (?); c0=1.0 mol/L, AB1=14.28 A (¦); c0=0.5 mol/L, AB1= 7.14 A, (?); c0=1.0 mol/L, AB1=7.14 A, (•).AB0=7.14 A. Symbols denote the GCMC simulation. 0.8 0.8 0.6 0.6 0.2 0.2 0 0 0 0 0.8 0.8 0.6 0.6 0.2 0.2 0 0 0 0.1 0.5 0.6 0 0.1 0.5 0.6 Lukšič et al. Modelling Electrolyte Adsorption in Nanoporous Materials 302 Acta Chim. Slov. 2006, 53, 292–305 adsorbed than, for example, doubling the concentration of matrix particles c0. Another way of presenting these results is often chosen in the electrochemical literature, i.e to calculate the so-called Donnan exclusion coefficient. This quantity is defined as: 27 Cbulk— < Ci > ,„. r =-----------------, (oj Cbulk where is the average concentration of model LiCl inside the matrix. Positive values of T indicate that electrolyte is excluded from the nanoporous material, which represents the basis of the desalination process, while negative values indicate electrolyte sorption. The results for T as a function of the bulk electrolyte concentration for various matrices are shown in Figure 13a-d. First we discuss Figures 13a and 13b showing T values for the case of adsorption in (a) c0=1.0 mol/L and (b) 0.5 mol/L matrices atAB1=7.14 A. It can be seen that in general the electrolyte was excluded (’rejected’) from the adsorbing material. An exception seems to be (cf panel a; triangles) matrices formed from chainlike ions (oligoelectrolyte) with related quenched counterions (case iv) at low values of cbulk. Also, exclusion of the model LiCl for the case of AB,1=7.14 A was greater for a dipolar matrix than for a chainlike ampholyte-type matrix in this range of observation. This difference was smaller for less concentrated matrices (panel b), which can be ascribed to the smaller exclusion volume effect. As expected, for all the examples studied here, the exclusion coefficient increased on increasing the electrolyte concentration. The same quantity for stronger electrostatic coupling (AB1=14.28 A) is shown Figures 13c-d. The situation is a little more complicated here; namely, at low electrolyte concentrations rejection was greater for polyampholyte matrices (stars and crosses) than for dipolar matrices (filled squares). This fact suggests that electrostatic interactions prevail over the exclusion volume term at smaller electrolyte concentrations and at higher couplings. The concentration of adsorbing electrolyte in case iii) was almost independent of chain length; it was only slightly smaller in the case of longer (N=32) chains. Interestingly, for adsorption in an electroneutral ionic matrix (case i , filled circles) intake (sorption) of invading electrolyte could be observed at higher electrostatic couplings (AB1 =14.28 A) and for low electrolyte concentrations as shown in panel (c). Sorption was stronger for oligoelectrolyte matrices (triangles) and could also be observed at lower concentrations for the case of AB,1=7.14 A (cf panel a). We can explain the sorption of the model electrolyte in examples i) and iv) by the fact that the adsorbing ions feel the separation of positive or negative charges and tend to compensate for them. In other words, an invading Li ion seeks to neutralize the negative matrix charge and oppositely, the invading Cl ion searches for ’pockets’ of positive charge. It is important that these ’pockets’ of positive and negative matrix charges be well separated; intake is expected to be bigger, for matrices where the separation between positive and negative charges is larger. This explains why the effect is the largest in case iv), where we have chains made of beads carrying exclusively positive charge and the corresponding negative counterions, occupying different domains of space. On the contrary, in the case of dipolar matrices this charge “separation” does not exist –an invading ion feels a dipole as, more or less, an electroneutral unit. The case of polyampholyte matrices lies between the cases ii) and iv). In summary, the strongest intake of electrolyte is expected in the case of oligoelectrolyte (charged oligoions) matrices, and at higher couplings also for the case of electroneutral ionic matrices. In all other cases, the electrolyte is rejected from the matrices: this effect is greater in more concentrated matrices and for a weaker Coulomb interaction (smaller ?B,1= 7.14 A). Hitherto, in all the cases discussed, the parameter defining the conditions of matrix prepara tion was fixed at ?B,0=7.14 A. We also carefully examined situations in which the matrix was prepared at higher Coulomb couplings, i.e. at ?B,0=14.28 A. We found that the conditions of matrix formation, in many cases, have only a minor impact on electrolyte adsorption. Of course, the effect of the matrix preparation is expected to be stronger when it yields distributions of obstacles with a more pronounced separation of the positive and negative charge within the matrix. In trying to understand the differences in adsorption in the various matrices from the structural point of view, we return to the radial distribution functions shown in subsection 4.3. It is enough to discuss two opposite cases; for the matrix containing polyampholyte chains iii) we have rejection of the electrolyte from the matrix (? > 0). In contrast to this a net intake of the electrolyte, i.e. ? < 0 is obtained for the example iv). The PDFS for these two cases should explain why we have rejection in one case and adsorption in another. While there is no crossing of the PDFS for Li-Li and Li-Cl in the case of polyampholyte matrices (not even at higher electrostatic interactions) (Fig. 8a-b), we observe this effect in the case of oligoelectrolyte matrices already at ?B,1=7.14 A. The ’charge inversion’ is more pronounced for stronger Coulomb coupling (Fig. 10a-b). Furthermore, PDFS for the annealed anion-positive charge of the matrix (or the annealed cation-negative charge of the matrix) for polyampholyte chains (case iii) are lower than unity for distances close Lukšič et al. Modelling Electrolyte Adsorption in Nanoporous Materials Acta Chim. Slov. 2006, 53, 292–305 303 0.4 c bulk 0.4 c bulk 0.4 c bulk "S^M^*"^ d 0.4 c bulk Figure 13. Donnan exclusion coefficient (T) for a LiCl model electrolyte confined in a disordered nanoporous material as a function of external bulk concentration. Figure a) c0 = 1.0 mol/L, XB1 = 7.14 A; b) c0= 0.5 mol/L, AB1=7.14 A; c) c0=1.0 mol/L, AB1=14.28 A; and d) c0=0.5 mol/L, AB1=14.28 A. Symbols denote: (•) disordered electroneutral ionic matrix; (¦) dipolar matrix; chain-like matrix with beads carrying alternating charges: (*) N=16; (+) N=32; and (A) chain like matrix with beads carrying positive charge neutralized by negative counterions (N=32). XB1=7.14 A. to the contact value. For the annealed anion it is not in a very favourable position (neither energetically, nor sterically) to make a close contact with a positive bead of the polyampholyte chain, since the latter is in contact with two negative beads. The same holds true for an annealed cation; the effect is more pronounced for weak Coulomb interactions and for adsorption in less concentrated matrices (see Fig. 9a-b). On the other hand, the fluid-matrix pair distribution functions, in the case of polyelectrolyte chains (iv) are higher than unity for the oppositely charged particles, which indicates attraction between the annealed ions and the matrix charges (Fig. 11a-b). The extremes in PDFS shown in Figures 10a-b and 11a-b indicate, on one hand, a) the accumulation of Cl ions around positively charged matrix chains, and b) the formation of triplets of Li ions with the quenched matrix anions, on the other. These two effects contribute toward increased adsorption, which is reflected in the negative values of exclusion coefficient ?, in contrast to the example (iii), where the coefficient ? is positive. 5. Conclusions Utilizing the grand canonical Monte Carlo calculations we examined the structural and thermodynamic properties of various partly quenched systems. The adsorbent-adsorbate system was composed of two subsystems: the so-called quenched phase (matrix) and an annealed fluid which under the conditions of observation was thermally equilibrated within quenched matrix particles (obstacles). The model for the annealed fluid was in all cases the same, i.e. the primitive model electrolyte, mimicking LiCl solution. It was assumed that thermal equilibrium was established with a surrounding bath of LiCl solution. The main purpose of this paper was to examine the adsorption of model electrolyte in matrices with various spatial distributions of the quenched charges. The distribution of annealed fluid within four different matrices was examined: i) a matrix formed by a rapid quench of some equilibrium distribution of a primitive model +1:-1 electrolyte at conditions 0.3 0.3 0.2 0.2 0.1 0.1 -0.1 -0.1 -0.2 -0.2 -0.3 -0.3 0.1 0.2 0.3 0.5 0.6 0.7 0.8 0.1 0.2 0.3 0.5 0.6 0.7 0.8 0.3 0.3 0.2 0.2 0.1 0.1 -0.1 -0.1 -0.2 -0.2 -0.3 -0.3 0.1 0.2 0.3 0.5 0.6 0.7 0.8 0.1 0.2 0.3 0.5 0.6 0.7 0.8 Lukšič et al. Modelling Electrolyte Adsorption in Nanoporous Materials 304 Acta Chim. Slov. 2006, 53, 292–305 T0 and ?, ii) a matrix consisting of fused hard spheres with opposite unit charge, iii) a matrix mimicking a quenched polyampholyte and represented as a collection of chainlike molecules with alternating unit charge on the beads, iv) a matrix composed of quenched charged chainlike molecules (oligoelectrolyte) and corresponding quenched counterions. The solvent in all the cases was treated as a dielectric continuum. Structural data for the adsorbed fluid, presented in the form of various pair distribution functions, show that the behaviour of adsorbed LiCl was different from that for the bulk electrolyte. The confinement effect was more pronounced for higher electrostatic interactions. Crossovers of like and unlike ion-ion PDFS were observed, indicating ’charge inversion’ for strong Coulomb interactions. The thermodynamic parameter of interest in this study was the excess chemical potential (mean activity coefficient) of the adsorbed electrolyte or, equivalently, the Donnan exclusion coefficient. The results indicate that electroneutral adsorbents with charged obstacles could either reject or adsorb the electrolyte; in other words, the Donnan exclusion coefficient could be positive or negative. Whether we have an intake of electrolyte or a rejection depends primarily on the spatial distribution of the charges in the matrix, and on the strength of the electrostatic interaction dictated by the dielectric constant of the solvent and also, to lesser degree, on the concentrations of the two phases. For small activities of the adsorbing electrolyte and for small concentration of the matrix particles (smaller exclusion volume effect), intake of LiCl was observed in cases i) and iv). On the other hand, a rejection of the electrolyte was observed for ii) and iii) over the whole range of observations, as also for higher activities of the bulk LiCl solution in the cases i) and iv). Clearly, electrolyte adsorption was strongly correlated with the spatial distribution of the matrix charges. In the examples i) and iv), the spatial separation of quenched positive and negative charges was compensated by an invading electrolyte and sorption of LiCl took place. The situation was more expressed in the latter case iv), since there was stronger separation of the positive (charged chains) and negative charges (counterions) than in case i). The charge separation was much smaller in the case of dipolar ii) or polyampholyte type of matrices iii), where individual obstacles are electroneutral, and consequently in these cases electrolyte rejection was observed. The exclusion volume effect dominated over electrostatic effects. In conclusion, we may say that the amount of electrolyte adsorbed depends strongly on the spatial distribution of the matrix charges, on the conditions of observation, given by the dielectric constant of the solvent, as well as by on the concentrations of all components present. Intake of electrolyte can be expected for matrices where separation between the domains of the positive and negative charge is large, and is amplified with increasing Coulomb coupling. 6. Acknowledgments This work was supported by the Slovenian Research Agency through Physical Chemistry Research Programme 0103-0201 and Research Project J1-6653. The authors wish to thank Dr. Jurij Reščič for the help with simulations of chainlike molecules. 7. References 1. N. K. Raman, M. T. Anderson, C. J. Brinker, Chem. Mater. 1996, 8, 1682–1701. 2. M. E. Davis, Nature 1993, 364, 391–393. 3. G. A. Ozin, C. Gil, Chem. Rev. 1989, 89, 1749–1764. 4. B. Hribar, O. Pizio, A. Trokhymchuk, V. Vlachy, J. Chem. Phys. 1997, 107, 6335–6341. 5. B. Hribar, O. Pizio, A. Trokhymchuk, V. Vlachy, J. Chem. Phys. 1998, 109, 2480–2489. 6. B. Hribar, V. Vlachy, O. Pizio, A. Trokhymchuk, J. Phys. Chem. B 1999, 103, 5361–5369. 7. B. Hribar, V. Vlachy, O. Pizio, J. Phys. Chem. B 2000, 104, 4479–4488. 8. B. Hribar, V. Vlachy, O. Pizio, J. Phys. Chem. B 2001, 105, 4727–4734. 9. B. Hribar, V. Vlachy, O. Pizio, Molec. Phys. 2002, 100, 3093–3103. 10. H. Dominguez, B. Hribar, V. Vlachy, O. Pizio, Physica A 2003, 324, 469–483. 11. V. Vlachy, H. Dominguez, O. Pizio, J. Phys. Chem. B, 2004, 108, 1046–1055. 12. W. G. Madden, E. G. Glandt, J. Stat. Phys. 1988, 51, 537–558. 13. W. G. Madden, J. Chem. Phys. 1992, 96, 5422–5432. 14. J. A. Given, G. Stell, J. Chem. Phys. 1992, 97, 4573–4574. 15. J. A Given, G. Stell, Physica A 1994, 209, 495–510. 16. J. A. Given, J. Chem. Phys. 1995, 102, 2934–2945. 17. H. Tatlipinar, G. Pastore, M. P. Tosi, Phil. Mag. Lett. 1993, 68, 357–361. 18. M. L. Rosinberg, G. Tarjus, G. Stell, J. Chem. Phys. 1994, 100, 5172–5177. 19. M. L. Rosinberg, in: C. Caccamo, J.-P. Hansen, G. Stell (Eds.): New Approaches to Problems in Liquid State Theory, Kluwer, Dordrecht, 1999, pp. 245–278. 20. O. Pizio, in: M. Borowko (Ed.): Computational Methods in Surface and Colloid Science, Marcel Dekker, New York, 2000, pp. 293–345. 21. R. D. Kaminsky, P. A. Monson, J. Chem. Phys. 1991, 95, 2936–2948. Lukšič et al. Modelling Electrolyte Adsorption in Nanoporous Materials Acta Chim. Slov. 2006, 53, 292–305 305 22. D. Bratko, A. K. Chakraborty, J. Chem. Phys. 1996, 104, 7700–7712. 23. D. M. Ford, E. D. Glandt, J. Chem. Phys. 1994, 100, 2391–2393. 24. E. Kierlik, M. L. Rosinberg, G. Tarjus, P. A. Monson, J. Chem. Phys. 1997, 106, 264–279. 25. G. M. Torrie, J. P. Valleau, J. Chem. Phys. 1980, 73, 5807–5816. 26. G. M. Torrie, J. P. Valleau, J. Phys. Chem. 1982, 86, 3251–3257. 27. B. Jamnik, V. Vlachy, J. Am. Chem. Soc. 1993, 115, 660–666. 28. D. Wu, K. Hui, D. Chandler, J. Chem. Phys. 1992, 96, 835–841. 29. K. S. Page, P. A. Monson, Phys. Rev. E 1996, 54, R29–R32. 30. L. Sarkisov, P. A. Monson, Phys. Rev. E 2000, 61, 7231–7234. 31. P. Linse, MOLSYM, Version 3.6.7 2000 (Lund University, Sweden). 32. V. Vlachy, D. Dolar, J. Chem. Phys. 1982, 76, 2010– 2014. 33. D. Bratko, V. Vlachy, Chem. Phys. Letters 1982, 90, 434–438. 34. G. Kortum, Treatise on Electrochemistry, Elsevier, Amsterdam, 1965. 35. V. Vlachy, T. Ichiye, A. D. J. Haymet, J. Am. Chem. Soc. 1991, 113, 1077–1082. Povzetek S pomočjo velekanonične Monte Carlo simulacije smo proučevali strukturne in termodinamske lastnosti raztopin modelnega elektrolita v neurejeni porozni snovi, ki vsebuje elektronevtralno kombinacijo nabojev. Porozno snov (adsorbent) smo si zamislili kot i) ravnotežno razporeditev ionov v +1: -1 elektrolitu pri neki koncentraciji c0, ii) kot sistem dipolov (delec sestavljata pozitivno in negativno nabiti kroglici), iii) kot zbir verig, sestavljenih iz togih kroglic z alternirajočim nabojem (poliamfolit), ter iv) kot sistem pozitivno nabitih verigastih molekul (oligoelektrolit) s pripadajočimi protiioni. Modelni elektrolit (adsorbat) je bil v termodinamskem ravnotežju z adsorbentom (delci adsorbenta so bili zamrznjeni v prostoru) ter zunanjo raztopino enake kemijske sestave. Topilo smo v vseh primerih obravnavali kot dielektrični kontinuum. Zanimalo nas je, kakšen vpliv ima razporeditev nabojev v porozni snovi na srednji koeficient aktivnosti adsorbiranega modelnega +1:-1 elektrolita. Izvedli smo računalniško simulacijo za niz parametrov, ki določajo stanje modela: spreminjali smo koncentracijo in razporeditev nabojev nanoporozne snovi ter koncentracijo modelnega elektrolita, pogoje priprave adsorbenta ter pogoje pri katerih sistem opazujemo. Računi so potrdili naša pretekla dognanja, da termodinamične lastnosti adsorbata močno zavisijo od temperature in dielektrične konstante topila, kakor tudi od koncentracije delcev adsorbenta ter adsorbata. Vpliv razporeditve nabojev adsorbenta na adsorpcijo elektrolita je pomemben, in še zlasti izrazit pri močnejših Coulombskih interakcijah. Lukšič et al. Modelling Electrolyte Adsorption in Nanoporous Materials