Adsorption of Simple Gases in MCM-41 Materials: The Role of Surface

Langmuir , 2006, 22 (1), pp 194–202 ... Molecular Modeling and Adsorption Properties of Ordered Silica-Templated CMK .... Langmuir 0 (proofing), ...
0 downloads 0 Views 614KB Size
Download PDF
194

Langmuir 2006, 22, 194-202

Adsorption of Simple Gases in MCM-41 Materials: The Role of Surface Roughness Benoit Coasne,*,†,‡ Francisco R. Hung,†,⊥ Roland J.-M. Pellenq,§ Flor R. Siperstein,| and Keith E. Gubbins† Center for High Performance Simulation and Department of Chemical and Biomolecular Engineering, North Carolina State UniVersity, Raleigh, North Carolina 27695-7905, Laboratoire de Physicochimie de la Matie` re Condense´ e, CNRS and UniVersite´ Montpellier 2, Place Euge` ne Bataillon, 34095 Montpellier Cedex 05, France, Departament d’Enginyeria Quı´mica, UniVersitat RoVira i Virgili, Campus Sescelades AV. dels Paı¨sos Catalans, 26, 43007 Tarragona, Spain, and Centre de Recherche sur la Matie` re Condense´ e et Nanosciences, CNRS Campus de Luminy, Case 913, 13288 Marseille Cedex 09, France ReceiVed June 22, 2005. In Final Form: October 19, 2005 This paper reports the development and testing of atomistic models of silica MCM-41 pores. Model A is a regular cylindrical pore having a constant section. Model B has a surface disorder that reproduces the morphological features of a pore obtained from an on-lattice simulation that mimics the synthesis process of MCM-41 materials. Both models are generated using a similar procedure, which consists of carving the pore out of an atomistic silica block. The differences between the two models are analyzed in terms of small angle neutron scattering spectra as well as adsorption isotherms and isosteric heat curves for Ar at 87 K and Xe at 195 K. As expected for capillary condensation in regular nanopores, the Ar and Xe adsorption/desorption cycles for model A exhibit a large hysteresis loop having a symmetrical shape, i.e., with parallel adsorption and desorption branches. The features of the adsorption isotherms for model B strongly depart from those observed for model A. Both the Ar and Xe adsorption branches for model B correspond to a quasicontinuous pore filling that involves coexistence within the pore of liquid bridges and gas nanobubbles. As in the case of model A, the Ar adsorption isotherm for model B exhibits a significant hysteresis loop; however, the shape of the loop is asymmetrical with a desorption branch much steeper than the adsorption branch. In contrast, the adsorption/desorption cycle for Xe in model B is quasicontinuous and quasireversible. Comparison with adsorption and neutron scattering experiments suggests that model B is too rough at the molecular scale but reproduces reasonably the surface disorder of real MCM-41 at larger length scales. In contrast, model A is smooth at small length scales in agreement with experiments but seems to be too ordered at larger length scales.

1. Introduction The siliceous mesoporous materials MCM-411 and SBA-152 attract a great deal of attention because of their possible uses as adsorbents or catalytic supports for gas adsorption, phase separation, catalysis, preparation of nanostructured materials, drug delivery, etc.3-5 These materials are obtained by a template mechanism involving the formation of surfactant or block copolymers micelles in a mixture composed of a solvent and a silica source. After polymerization of the silica and removal of the organic micelles, one obtains a material made up of a hexagonal array of cylindrical pores. The pore size distribution is narrow with an average value that can be varied from 1.5 up to 20 nm, depending on the synthesis conditions.5 Properties of MCM-41 and SBA-15 materials have been extensively studied * To whom correspondence should be addressed. E-mail: bcoasne@ lpmc.univ-montp2.fr. Phone: +33 4 67 14 33 78. Fax: +33 4 67 14 42 90. † North Carolina State University. ‡ Laboratoire de Physicochimie de la Matie ` re Condense´e (UMR 5617CNRS). § Centre de Recherche sur la Matie ` re Condense´e et Nanosciences (UPR 7251-CNRS). | Universitat Rovira i Virgili. ⊥ Present address: Department of Chemical and Biological Engineering, University of Wisconsin, Madison, WI 53706-1691. (1) Beck, J. S.; Vartulli, J. C.; Roth, W. J.; Leonowicz, M. E.; Kresge, C. T.; Schmitt, K. D.; Chu, C. T.-W.; Olson, D. H.; Sheppard, E. W.; McCullen, S. B.; Higgins, J. B.; Schlenker, J. L. J. Am. Chem. Soc. 1992, 114, 10834-10843. (2) Zhao, D.; Feng, J.; Huo, Q.; Melosh, N.; Fredrickson, G. H.; Chmelka, B. F.; Stucky, G. D. Science 1998, 279, 548-552. (3) Corma, A. Chem. ReV. 1997, 97, 2373-2419. (4) Ciesla, U.; Schu¨th, F. Microp. Mesop. Mater. 1999, 27, 131-149. (5) Soler-Illia, G. J. de A. A.; Sanchez, C.; Lebeau, B.; Patarin, J. Chem. ReV. 2002, 102, 4093-4138.

by combining transmission electronic microscopy (TEM), X-ray diffraction, and adsorption experiments.3,4 However, some uncertainties remain regarding the surface chemistry (presence of impurities, defects, etc.) as well as the surface texture (microporosity, surface roughness) of these materials.6-8 For instance, it has been shown that the cylindrical mesopores in SBA-15 materials are connected through transverse microporous channels.9-11 On the other hand, it is generally agreed that MCM41 is made up of unconnected mesopores.10,12 From a fundamental point of view, MCM-41 and SBA-15 silica pores are considered as model materials to investigate the effect of anisotropic nanoconfinement on the thermodynamic properties of fluids. As a result, many experimental, simulation, and theoretical works have been reported on gas adsorption in these materials (for a review, see ref 13). In most simulation and theoretical studies, the silica MCM-41 and SBA-15 materials (6) Sonwane, C. G.; Bhatia, S. K.; Calos, N. J. Langmuir 1999, 15 (5), 46034612. (7) Berenguer-Murcia,A.; Garcia-Martinez, J.; Cazorla-Amoros, D.; MartinezAlonso, A.; Tascon, J. M. D.; Linares-Solano, A. In Studies in Surface Science and Catalysis 144; Rodriguez-Reinoso, F., McEnaney, B., Rouquerol, J., Unger, K., Eds.; Elsevier Science: Amsterdam, 2002; pp 83-90. (8) Edler, K. J.; Reynolds, P. A.; White, J. W. J. Phys. Chem. B 1998, 102, 3676-3683. (9) Imperor-Clerc, M.; Davidson, P.; Davidson, A. J. Am. Chem. Soc. 2000, 122, 11925-11933. (10) Ryoo, R.; Ko, C. H.; Kruk, M.; Antochshuk, V.; Jaroniec, M. J. Phys. Chem. B 2000, 104, 11465-11471. (11) Galarneau, A.; Cambon, H.; Di Renzo, F.; Fajula, F. Langmuir 2001, 17, 8328-8335. (12) Jun, S.; Joo, S. H.; Ryoo, R.; Kruk, M.; Jaroniec, M.; Liu, Z.; Ohsuna, T.; Terasaki, O. J. Am. Chem. Soc. 2000, 122, 10712-10713. (13) Gelb, L. D.; Gubbins, K. E.; Radhakrishnan R.; Sliwinska-Bartkowiak, M. Rep. Prog. Phys. 1999, 62, 1573-1659.

10.1021/la051676g CCC: $33.50 © 2006 American Chemical Society Published on Web 11/24/2005

Adsorption of Gases in MCM-41 Materials

are modeled as pores with a regular cylindrical geometry. On the other hand, recent simulation studies14-20 have investigated the behavior of pure fluids and mixtures confined in atomistic pore models of these materials. The effect of surface roughness and morphological defects, such as constrictions, on adsorption and capillary condensation of gases was also recently addressed.18,19 However, it is not clear whether the atomistic silica pore models used in the simulation studies mentioned above are realistic representations of the morphology of MCM-41 pores, due to the lack of conclusive experiments regarding their surface roughness and morphology. Some experimental results7,21 attribute some surface roughness to the pore walls, whereas others6,8 suggest a smooth surface at molecular scales (3-7 Å) with roughness and other morphological defects (constrictions and tortuosity) at larger length scales (20-50 Å). One possible approach to modeling porous solids is to mimic the synthesis process of the real material, a strategy used in the past15 to develop realistic models for Vycor and controlled pore glasses. Recently, Siperstein and Gubbins23,24 reported on-lattice Monte Carlo simulations of surfactant-solvent-silica systems that were able to reproduce the formation of hexagonal structure, resembling the arrangement of silica MCM-41 pores.23,24 Moreover, a simulation study of gas adsorption in the coarse-grained structure obtained by the authors was in fair agreement with experiments.23 The motivation for this work is two-fold. We report the preparation of an atomistic silica mesopore that keeps the morphological features of the MCM-41 lattice model developed by Siperstein and Gubbins. The porous material is generated by carving out of a silica block the skeleton of a MCM-41 pore that was obtained by lattice Monte Carlo simulations. In doing this “downscaling” process, atomic details can be included, more accurate potentials can be used for gas adsorption, and the effect of structural defects on the adsorbate phase behavior can be assessed. In the second part of this work, we present grand canonical Monte Carlo (GCMC) simulations of Ar and Xe adsorption in our model of MCM-41 pores. We discuss the effect of the pore shape by comparing the adsorption isotherm and the isosteric heat curves with those obtained for an atomistic, regular cylindrical nanopore having a similar pore diameter. When possible, comparison with experimental data is made. The remainder of the paper is organized as follows. In section 2, we present the numerical method used to generate atomistic silica MCM-41 pores as well as the GCMC technique. The simulation results are discussed in section 3. Section 4 contains concluding remarks and suggestions for future work.

2. Computational Details 2.1. Preparation of Realistic MCM-41 Silica Pores. Atomistic pores were generated using the method proposed by Pellenq and Levitz to prepare numerical Vycor samples.25 Coasne and Pellenq showed that this technique can be used to prepare pores (14) Maddox, M. W.; Olivier, J. P.; Gubbins, K. E. Langmuir 1997, 13, 17371745. (15) Gelb, L. D. Mol. Phys. 2002, 100, 2049-2057. (16) Coasne, B.; Grosman, A.; Ortega, C.; Pellenq, R. J.-M. In Studies in Surface Science and Catalysis 144; Rodriguez-Reinoso, F., McEnaney, B., Rouquerol, J., Unger, K., Eds.; Elsevier Science: Amsterdam, 2002; pp 35-42. (17) He, Y.; Seaton, N. A. Langmuir 2003, 19, 10132-10138. (18) Coasne, B.; Pellenq, R. J.-M. J. Chem. Phys. 2004, 120, 2913-2922. (19) Coasne, B.; Pellenq, R. J.-M. J. Chem. Phys. 2004, 121, 3767-3774. (20) Coasne, B.; et al. Langmuir 2005, in preparation. (21) Fenelonov, V. B.; Derevyankin, A. Y.; Kirik, S. D.; Solovyov, L. A.; Shmakov, A. N.; Bonardet, J.-L.; Gedeon, A.; Romannikov, V. N. Microp. Mesop. Mater. 2001, 44-45, 33-40. (22) Gelb, L. D.; Gubbins, K. E. Langmuir 1998, 14, 2097-2111 (23) Siperstein, F. R.; Gubbins, K. E. Mol. Simul. 2001, 27, 339-352. (24) Siperstein, F. R.; Gubbins, K. E. Langmuir 2003, 19, 2049-2057. (25) Pellenq, R. J.-M.; Levitz, P. E. Mol. Phys. 2002, 100, 2059-2077.

Langmuir, Vol. 22, No. 1, 2006 195

Figure 1. Plane and transverse section views of an atomistic cylindrical pore having a constant radius R0 ) 3.2 nm (model A). White and gray spheres are oxygen and silicon atoms, respectively. Black spheres are hydrogen atoms that delimit the pore surface. Two simulation boxes aligned along the pore axis have been represented for visualization purpose.

of various morphologies and/or topologies, such as cylindrical, hexagonal, ellipsoidal, and constricted pores.16,18,19 Recently, He and Seaton used a similar procedure to generate cylindrical silica pores.17 In the present work, two different models of MCM41 pores were generated. Model A consists of a regular cylindrical pore having a constant section of radius R0 ) 3.2 nm. Model B has the same average pore size but exhibits an important surface disorder that reproduces the morphological features of an onlattice simulation mimicking the synthesis process of MCM-41 pores. Model A. To generate the first MCM-41 nanopore (model A), we carved out of an atomistic silica block a regular cylindrical pore of radius R0 ) 3.2 nm. The initial silica matrix is made of 13 × 13 × 15 unit cells of nonporous cristoballite, corresponding to a simulation box of 9.3 nm × 9.3 nm × 10.7 nm. To mimic the pore surface in a realistic way, we first removed the Si atoms that are in an incomplete tetrahedral environment. In a second step, we removed all oxygen atoms that are nonbonded. This procedure ensures that the remaining silicon atoms have no dangling bonds and the remaining oxygen atoms have at least one saturated bond with a Si atom. Then, the electroneutrality of the simulation box was ensured by saturating all oxygen dangling bonds with hydrogen atoms. The latter are placed in the pore void, perpendicularly to the pore surface, at a distance of 1 Å from the closest unsaturated oxygen atom. It has been shown15 that the density of OH groups obtained using such a procedure is close to that obtained experimentally for porous silica glasses (7-8 OH per nm2).26 Then, we displace slightly and randomly all of the O, Si, and H atoms in order to mimic an amorphous silica surface (the maximum displacement in each direction x, y, and z is 0.7 Å). The resulting atomistic porous material is shown in Figure 1. Following the previous work by Pellenq and Levitz,25 we relaxed the porous structure obtained after randomization of the surface by performing Monte Carlo simulations in the canonical NVT ensemble. The interactions between atoms i and j separated by a distance rij were calculated using two-body potentials of the BKS type27

Uij )

( )

qiqj rij f6,ij(rij)C6,ij + Aij exp rij Fij r6

(1)

ij

The potential function given by eq 1 is the sum of the Coulomb interaction, the short-range repulsion, and the dispersion interaction. The Coulomb energy was computed without the use of Ewald summation. Indeed, Puibasset and Pellenq28 have shown (26) Landmesser, H.; Kosslick, H.; Storek, W.; Frick, R. Solid State Ionics 1997, 101-103, 271-277. (27) van Beest, B. W. H.; Kramer, G. J.; van Santen, R. A. Phys. ReV. Lett. 1990, 64, 1955-1958. (28) Puibasset, J.; Pellenq, R. J.-M. Phys. Chem. Chem. Phys. 2004, 6, 19331937.

196 Langmuir, Vol. 22, No. 1, 2006

Coasne et al.

Table 1. Potential Parameters Used for the Surface Relaxation of the MCM-41 Silica Porea partial charges (e) +2.4 +0.6 -1.2 -1.2

Si H O (SiO) O (OH) Buckingham potential

A (eV)

F (Å)

C (eV‚Å6)

1.2-

-O (ref 27) O1.2--O1.2- (ref 27) H0.6+-O1.2- (ref 31)

18003.7572 1388.7730 311.97000

0.20520 0.36232 0.25000

133.5381 175.0000 0.000000

Morse potential

D (eV)

R (Å-1)

r0 (Å)

7.05250

3.17490

0.94285

Si2.4+

-O

0.6+

H

1.2-

(ref 31)

a

Oxygen atoms forming a OH bond at the pore surface and oxygen atoms in the bulk region have the same properties in this model. All interactions were determined using two-body potentials of the BKS type,27 except the interaction between H and O atoms forming a hydroxyl group that was modeled using a Morse potential.30,31 Note that the following interactions are purely Coulombic: Si/Si, Si/H, and H/H. This set of parameters is taken from refs 27 and 31.

that the correction due to the long-range contribution is very small for a large simulation box (we note that the box size in the present work is even larger than that in ref 28). The f6,ij(rij) term in eq 1 is a damping function that avoids the divergence of the dispersive energy when rij tends to 0 (this divergence comes from the mathematical expansion in power of 1/r2n that becomes nonvalid for small distances r).29 The function f6,ij(rij) has the following form: 6

f6,ij(rij) ) 1 -

∑ k)0

(rij/Fij)k k!

( )

exp -

rij

Fij

(2)

where Fij is characteristic of the overlap at short distances of the wave functions of the two interacting atoms (see eq 1). The damping function in eq 2 varies steeply from 0 up to 1 as the distance rij between the two atoms increases. In the relaxation process, all interactions were determined using eq 1, except the interaction between H and O atoms forming a hydroxyl group that was modeled using a Morse potential30,31

Uij ) D(1 - exp[-R(rij - r0)])2 - D

(3)

The values of parameters for eqs 1 and 3 as well as the partial charges used in the simulation of the relaxation process are reported in Table 1. To accelerate the simulation run, only atoms located at a distance smaller than 15 Å from the H atoms at the pore surface were allowed to move. The canonical Monte Carlo simulations were performed at a temperature of 2500 °C. This high temperature (the melting point of pure silica is ∼1700 °C) was used to compensate the fact that Monte Carlo simulations, which allow the displacement of one atom at a time only, are not efficient in permitting bond breaking between Si-O and O-H atoms. Indeed, chemical bonding involves energy barriers of a few eV, which are difficult to overcome using Monte Carlo simulations. As a result, we decided to use a quite large temperature to allow the particles to move apart from each other by giving more thermal energy to the system. The temperature of 2500 °C is of the same (29) Tang, K. T.; Toennies, J. P. J. Chem. Phys. 1984, 80, 3726-3741. (30) Gale, J. D. Philos. Mag. B 1996, 73, 3-19. (31) de Leeuw, N. H.; Manon Higgins, F.; Parker, S. C. J. Phys. Chem. B 1999, 103, 1270-1277.

Figure 2. On-lattice silica structure reproducing the hexagonal array of independent cylindrical pores. 2 × 2 × 1 simulation boxes are shown. For visualisation purposes, the porous matrix has been rotated by 45° along the z direction to be oriented perpendicularly to the axis of the pores (adapted from ref 23).

order of magnitude as that used in previous simulations of amorphous silica.32 The maximum displacement allowed for the particle moves was initially set to δ ) 0.5 Å. This rather large step in particle displacement allows the stability of the porous structure to be checked by sampling configurations where chemical bonds are broken (the equilibrium distances of the OH and SiO bonds are 0.9 and 1.6 Å, respectively). The structure was relaxed using this maximum displacement until the relative change in energy was less than 0.1%. Afterward, we refined the equilibration process by gradually reducing the parameter δ (each time the simulation was performed until ∆E/E was less than 0.1%). The relaxation simulation was considered complete when the previous criterion over the change in energy was reached with δ ) 0.05 Å. The length of the simulation was such that at least 105 Monte Carlo steps per atom were performed. Model B. In a recent work, Siperstein and Gubbins23,24 studied the phase behavior of surfactant-solvent-silica containing systems using on-lattice Monte Carlo simulations. Depending on the concentration of the system, the authors observed different silica structures that reproduce the series of mesoporous M41 materials.5 Of particular interest for the present work, it was found that the system forms a hexagonal phase mimicking the MCM-41 porous structure, for a surfactant/silica concentration ratio equal to 0.2 (see Figure 2). To estimate the pore size distribution and the silica/void interface in a rigorous way, the authors superimposed to the model material a grid having a spacing of 0.1 lattice unit. The largest sphere that could be inserted without overlapping with any of the beads representing the onlattice silica structure was calculated for each grid site. The number N(r) of spheres having a radius r that were inserted using such a procedure provides an accurate estimate of the pore size distribution of the porous material. Assuming a lattice unit of 0.5 nm, the lattice porous structure shown in Figure 2 was found to have a pore length of 28 nm and a mean pore radius of about 3.2 nm.23 The choice of such a lattice spacing was motivated by the fact that it corresponds to a reasonable approximation for the size of the bead in the on-lattice simulations (a few SiO2 tetrahedron units per bead). Following the general procedure described above for model A, we generate the atomistic silica nanopore model B by carving out of a silica block the porous matrix that corresponds to the lattice structure shown in Figure (32) Binder, K. Comput. Phys. Comm. 1999, 121-122, 168-175 and references therein.

Adsorption of Gases in MCM-41 Materials

Figure 3. Plane and transverse section views of an atomistic MCM41 pore corresponding to model B. This configuration was obtained after relaxation of the structure obtained by carving out of a silica matrix one of the pores shown in the lattice structure in Figure 2. White and gray spheres are oxygen and silicon atoms, respectively. Black spheres are hydrogen atoms that delimit the pore surface.

2. The initial silica matrix is made of 19 × 19 × 39 unit cells of nonporous cristoballite, corresponding to a simulation box of 13.5 nm × 13.5 nm × 27.8 nm. This method enables us to design fully atomistic pores that keep the morphological features of the original template mesoporous material obtained by Siperstein and Gubbins. We isolated one of the pores shown in Figure 2 instead of the complete hexagonal array of pores. The porous material that we obtained using such a procedure was relaxed using NVT Monte Carlo simulations, as for the case of model A. The resulting atomistic porous material (model B) is shown in Figure 3. Characterization of the Porous Materials. The pore size distribution of the MCM-41 pore model B was estimated as the distribution of the distances between hydrogen atoms at the pore surface and the center of the pore (the latter was estimated as the center of mass of the distribution of hydrogen atoms). This pore size distribution is shown in Figure 4 and has a mean value of 3.2 nm and a dispersion of (1 nm. We note that this mean value differs a little from the most probable pore size R ) 3.5 nm, due to the slightly asymmetrical shape of the pore size distribution (which is not exactly Gaussian). To characterize the surface roughness, we also calculated the average pore radius as a function of the position on the pore axis Z (Figure 4). It was found that the average pore dispersion along the pore axis is about (0.5 nm. We prepared a second porous sample of model B using another pore obtained from the on-lattice structure. The two atomistic samples were found to have very similar pore size distributions and surface roughnesses. Consequently, the nanopore shown in Figure 3 was considered representative of the porous solid. Following the work by Pellenq et al.,33 small angle neutron scattering spectra (SANS) were calculated for each porous material considered in this work. Such spectra I(Q) provide crucial information regarding the surface properties of the material in (33) Pellenq, R. J.-M.; Rodts, S.; Pasquier, V.; Delville, A.; Levitz, P. Adsorption 2000, 6, 241-249.

Langmuir, Vol. 22, No. 1, 2006 197

Figure 4. (Top) Pore size distribution for model B estimated from the distance R between hydrogen atoms at the pore surface and the pore center. (Bottom) Average pore size R(Z) as a function of the position on the pore axis Z for model B. R(Z) was estimated from the average distance to the pore center of hydrogen atoms located at a position Z.

Figure 5. Small angle neutron scattering spectra for pore model A (black line) and pore model B (grey line). The spectrum for pore model A has been multiplied by 10 for the sake of clarity. The Bragg peaks of cubic cristobalite are in very good agreement with the experimental X-ray peaks: 1.52, 2.49, 2.90, 3.03, 3.83, and 4.55 Å-1. The two straight lines correspond to an algebraic decay over the range 0.05-0.55 Å-1 with the exponent 4.0 for model A and 3.6 for model B (see text).

the Porod range, i.e., QD > 1 (D is the width of the pore).34 SANS spectra are shown in Figure 5 for models A and B. Peaks obtained at large momentum transfer (Q > 1 Å-1) are in very good agreement with experimental X-ray data on the structure of cristoballite: 1.52, 2.49, 2.90, 3.03, 3.83, and 4.55 Å-1. Spectra were fitted using an algebraic decay law I(Q) ) Q-x over the range 0.05-0.55 Å-1 in order to determine the Porod exponent that characterizes the surface roughness of the porous solids at length scales ∼12-50 Å. x was found to be equal to 4.0 ( 0.1 in the case of model A, as expected for a cylindrical pore having a smooth pore/void interface that exhibits only atomistic surface roughness. In contrast, the exponent x for pore model B is equal to 3.6 ( 0.1, which is close to the value found for disordered porous silica samples such as Vycor35,36 or silica gel.6,37 Such (34) Ramsay, J. D. F. AdV. Colloid Interface Sci. 1998, 76-77, 13-37. (35) Levitz, P.; Ehret, G.; Sinha, S. K.; Drake, J. M. J. Chem. Phys. 1991, 95, 6151-6161.

198 Langmuir, Vol. 22, No. 1, 2006

a value is in agreement with SAXS or SANS experimental investigations of the surface disorder for MCM-41 materials.6,8 Edler et al.8 reported that the Porod exponent of MCM-41 samples obtained from the ordinary preparation is in the range 3-3.5. This result was confirmed later by Sonwane et al.,6 who showed that the Porod exponent for MCM-41 samples of different pore sizes is in the range 3.4-3.6 at length scales 20-50 Å. 2.2. Grand Canonical Monte Carlo Simulation of Gas Adsorption. To test the atomistic MCM-41 pores that we built, we performed GCMC simulations of Ar and Xe adsorption and compared our results with experimental data. The GCMC technique is a stochastic method that simulates a system having a constant volume V (the pore with the adsorbed phase), in equilibrium with an infinite fictitious reservoir of particles imposing its chemical potential µ and its temperature T.38-40 The absolute adsorption isotherm is given by the ensemble average of the number of adsorbed atoms versus the pressure of the gas reservoir P (the latter is obtained from the chemical potential according to the bulk equation of state for an ideal gas). Periodic boundary conditions were used along the pore axis, so that the system is equivalent to a pore of infinite length. Interactions between Ar or Xe atoms and substrate atoms were calculated using the PN-TraZ potential as reported for rare gas adsorption in zeolite41 or in porous silica glass.25 The intermolecular energy is written as the sum of the dispersion interaction with a repulsive short-range contribution and an induction term due to the interaction of the adsorbed atom with the local field created by the partial charges of the atoms in the substrate. Values of the interaction parameters as well as details of the intermolecular potential functions can be found in the previous work by Pellenq and Levitz.25 We calculate the adsorbate/substrate interaction by using an energy grid;42 the potential energy is calculated at each corner of each elementary cube (about 1 Å3). An accurate estimation of the energy is then obtained by a linear interpolation of the grid values. Such a procedure enables simulation of adsorption in mesoporous media of complex morphology and/or topology without a direct summation over matrix species in the course of GCMC runs.15,19,22,43 The Ar/Ar and Xe/Xe interactions were calculated using Lennard-Jones potentials with Ar ) 120 K, σAr ) 0.34 nm, and Xe ) 211 K, σXe ) 0.41 nm.44 The bulk saturating vapor pressure P0 for both Ar and Xe was determined according to the equation of state proposed by Kofke for Lennard-Jones fluids.45 The system was first allowed to equilibrate in the course of a first GCMC run. Afterward, the number of particles in the system and the isosteric heat of adsorption were averaged using 105 Monte Carlo steps per particle.

3. Results and Discussion In this section, we report molecular simulations of Ar adsorption at 87 K and Xe adsorption at 195 K in pore models A and B. Results are analyzed by comparing for both pores the adsorption isotherm and isosteric heat curves with experimental data. Filling (36) Mitropoulos, A. C.; Haynes, J. M.; Richardson, R. M.; Kanellopoulos, N. K. Phys. ReV. B 1995, 52, 10035-10042. (37) Drake, J. M.; Levitz, P.; Klafter, J. Isr. J. Chem. 1991, 31, 135-145. (38) Nicholson, D.; Parsonage, N. G. Computer Simulation and the Statistical Mechanics of Adsorption; Academic Press: New York, 1982. (39) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, 1987. (40) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed.; Academic Press: London, 2002. (41) Pellenq, R. J.-M.; Nicholson, D. J. Phys. Chem. 1994, 98, 13339-13349. (42) Pellenq, R. J.-M.; Nicholson, D. Langmuir 1995, 11, 1626-1635. (43) Pellenq, R. J.-M.; Levitz, P. E. Mol. Simul. 2001, 27, 353-370. (44) Streett, W. B.; Staveley, L. A. K. J. Chem. Phys. 1967, 47, 2449-2454. (45) Kofke, D. A. J. Chem. Phys. 1993, 105, 4149-4162.

Coasne et al.

Figure 6. (Top) Ar adsorption isotherms at 87 K obtained by GCMC simulations: (circles) pore model A, (squares) pore model B. (triangles) are experiments for MCM-41 pores with R ) 3.25 nm by Kruk and Jaroniec.53 Open and closed symbols are the adsorption and desorption data, respectively. P0 is the bulk saturation pressure of the gas. Adsorbed amounts have been normalized to the number of Ar atoms, N0, when the pore is completely filled. (Bottom) Xe adsorption isotherms at 195 K obtained by GCMC simulations: (circles) pore model A, (squares) pore model B. Open and closed symbols are the adsorption and desorption data, respectively. P0 is the bulk saturation pressure of the gas. Adsorbed amounts have been normalized to the number of Xe atoms, N0, when the pore is completely filled.

and emptying of each porous material are also investigated by examining configurations of the adsorbed atoms at different pressures. 3.1. Adsorption Isotherms. The Ar adsorption isotherms obtained at 87 K for the atomistic MCM-41 pore models A and B are shown in Figure 6. We also report the Xe adsorption isotherms at 195 K for both pore models. Adsorbed amounts have been normalized to the total number of atoms, N0, when the pores are filled. The uncertainty in N/N0 was estimated for each pressure from the statistical error 1/xNr for a Gaussian distribution of Nr values centered around its average value. The maximum uncertainty, which is found for the low-pressure range, is at most a few %. Of course, such an estimate does not include the error in the condensation/evaporation pressures, which is related to the use of a finite number of attempted configurations in the simulations (see the discussion below). The Ar and Xe adsorption isotherms for pore model A conform to the typical behavior observed for adsorption/condensation of gas in regular cylindrical pores;13,18,19,46 the adsorbed amount increases continuously in the multilayer adsorption regime until a jump occurs due to capillary condensation of the fluid within the pore. The condensation pressures are P ) 0.69P0 for Ar at 87 K and P ) 0.70P0 for Xe at 195 K. The evaporation pressures, P ) 0.27P0 for Ar and P ) 0.55P0 for Xe, are much lower than the condensation pressures, so that a large hysteresis loop is observed in both cases. The adsorption isotherms for model B exhibit very different features from those for model A. Adsorption/condensation of Ar and Xe in model B is a quasicontinuous mechanism that is composed of reversible paths and two small jumps in the adsorbed (46) Coasne, B.; Gubbins, K. E.; Pellenq, R. J.-M. Part. Part. Syst. Char. 2004, 21, 149-160.

Adsorption of Gases in MCM-41 Materials

amount. These two discontinuities, which occur at P ∼ 0.44P0 and P ∼ 0.51P0 for Ar and P ∼ 0.47P0, and P ∼ 0.56P0 for Xe, correspond to partial condensations of the fluid within the pore (see below). In contrast to the condensation mechanism, it is found that the evaporation mechanism for model B depends on the nature of the adsorbate. In the case of Ar at 87 K, the desorption for this pore is similar to that for model A; due to the infinite pore length, the confined fluid cannot evaporate at equilibrium and the pore empties in a one-step process through a cavitation effect.19,46-51 This spontaneous nucleation of the gas phase within the pore occurs at the same pressure for pore models A and B, i.e., P ∼ 0.25P0. Such a result suggests that the spinodal limit of Ar confined at 87 K does not depend much on the surface disorder for pore sizes about R0 ) 3.2 nm. The asymmetrical shape of the hysteresis loop for Ar in model B strongly departs from the symmetrical shape observed for model A. This result shows that asymmetrical hysteresis loops, which are usually observed for porous materials with connected pores, can be observed also for unconnected materials with, however, important surface disorder. This result is in agreement with our previous works18,19,46 where capillary hysteresis loops for unconnected pores with constrictions were found to be asymmetrical. The evaporation of Xe confined in pore model B is very different from that observed for Ar as desorption for Xe resembles the adsorption mechanism, i.e., pore emptying consists of a quasicontinuous process composed of continuous and discontinuous paths. The hysteresis loops that accompany the jumps in the adorbed amounts are much smaller than that observed in the case of the regular cylindrical pore. We note that the experimental observation of such two-step filling mechanisms may be difficult as the hysteresis loops are very narrow and the jumps in the adsorbed amounts are small. Also, it is not clear whether the hysteresis loops found in the simulations for Xe in model B are composed of long-life metastable states. One may expect such states along the adsorption/desorption branches to disappear if longer simulations or better sampling techniques were used; in this hypothetical situation, the hysteresis loops shown in Figure 6 would vanish and the condensation/evaporation processes would be truly continuous. Filling and emptying of both Ar and Xe in model B occur at much smaller pressures than those for model A. Again, this result is in agreement with previous works18,19,46 showing that an important surface disorder (ellipsoidal or constricted pores) induces a decrease in the capillary condensation pressure of confined fluids. It is also found that, at a given pressure, the adsorbed amount prior to capillary condensation is much larger for model B than for model A. This result agrees with our previous works,18,52 where it was shown that the surface roughness leads to an increase in the adsorbed amount in the multilayer adsorption regime. The results for Ar adsorption in pore models A and B are compared in Figure 6 with the experimental data obtained by Kruk and Jaroniec53 for MCM-41 pores with R0 ) 3.25 nm. The adsorption isotherm for model A is much closer to the experimental curve than that for model B. The condensation pressure for model B largely underestimates the experimental transition pressure. Also, the adsorbed amounts for this model (47) Sarkisov, L.; Monson, P. A. Langmuir 2001, 17, 7600-7604. (48) Ravikovitch, P. I.; Neimark, A. V. Langmuir 2002, 18, 1550-1560. (49) Vishnyakov, A.; Neimark,A. V. Langmuir 2003, 19, 3240-3247. (50) Van Der Voort, P.; Ravikovitch, P. I.; De Jong, K. P.; Benjelloun, M.; Van Bavel, E.; Janssen, A. H.; Neimark, A. V.; Weckhuysen, B. M.; Vansant, E. F. J. Phys. Chem B 2002, 106, 5873-5877. (51) Ravikovitch, P. I.; Neimark, A. V. Langmuir 2002, 18, 9830-9837. (52) Pellenq, R. J. M.; Coasne, B.; Levitz, P. E. In Adsorption and transport at the nanoscale; Quirke, N., Nicholson, D., Eds.; Taylor and Francis: London, 2005. (53) Kruk, M.; Jaroniec, M. Chem. Mater. 2000, 12, 222-230.

Langmuir, Vol. 22, No. 1, 2006 199

Figure 7. Transverse section view of Ar atoms adsorbed at 87 K in the MCM-41 pore model A (R0 ) 3.2 nm) at P ) 0.69P0. Black spheres are hydrogen atoms at the pore surface and gray spheres are Ar atoms. The two extremities of the pore are connected to each other, due to the use of periodic boundary condition along the pore axis. The black lines delimit the regions of low (gas) and high (adsorbed film) density.

overestimate the experimental values prior to capillary condensation. As will be discussed below, it seems that the atomic surface roughness in model B is too great. On the other hand, the condensation pressure for model A is close to the experiments. Also, the Ar adsorption isotherm for model A reproduces quite well the experimental adsorbed amounts at low pressure, which shows the ability of the PN-TraZ to model accurately the interactions between the adsorbate and the silica substrate. However, model A also overestimates at the onset of capillary condensation the adsorbed amounts that are observed in the experiment. A possible explanation for this result is the use in the simulation of pairwise potentials to estimate the adsorbateadsorbate and adsorbate-substrate interactions, which account only in an effective way for many body interactions in the system. Effective potentials, which are usually derived for bulk phases, can lead to an inaccurate estimate of the interaction energy for 2D systems such as the interface between the gas phase in the center of the pore and the adsorbed film at the pore surface.54 Fernandez-Alonso et al.55 have shown that three-body interactions are responsible for ∼15% of the total energy for Ar in silicalite. The authors also found that pairwise effective potentials overestimate the interaction energy calculated using potentials including many body terms. This deviation, which was observed for small zeolite channels, is expected to be less for the bigger MCM-41 pores; however, even a small difference in the estimate of the total energy may account for the disagreement between the adsorbed amounts for model A and the experimental data at the onset of capillary condensation. 3.2. Pore Filling and Emptying. A typical configuration of Ar atoms adsorbed at 87 K in pore model A is shown in Figure 7 for P ) 0.69P0 (at the onset of capillary condensation). The pore surface is covered by a molecularly thick cylindrical film until the capillary condensation pressure is reached. As expected for a regular nanopore, condensation occurs through an irreversible transition between the partially filled and the completely filled configurations. Configurations of Xe atoms adsorbed at 195 K (not shown) reveal that the filling of this pore model is identical to that for Ar at 87 K. As can be seen in Figure 8, adsorption and condensation in pore model B are significantly different from those for pore model A. At P ) 0.44P0, the confined fluid exhibits a region of low density that coexists with a region of high density; the latter region forms a liquidlike bridge within the pore that separates gas nanobubbles. We note that this result was corroborated by examining the density profiles along the pore axis.56 Such a (54) Nicholson, D. Surf. Sci. 1985, 151, 553-569. (55) Fernando-Alonso, F.; Pellenq, R. J.-M.; Nicholson, D. Mol. Phys. 1995, 86, 1021.

200 Langmuir, Vol. 22, No. 1, 2006

Figure 8. Transverse section views of Ar adsorbed at 87 K in the MCM-41 pore model B (R ) 3.2 nm) at P ) 0.44P0 (top) and P ) 0.46P0 (bottom). Black spheres are hydrogen atoms at the pore surface and gray spheres are Ar atoms. The two extremities of the pore are connected to each other, due to the use of periodic boundary condition along the pore axis. The black lines delimit the regions of low density that are not filled (nanobubbles).

gas-liquid coexistence has been already reported for adsorption in chemically heterogeneous slit pores,57,58 cylindrical pores with constrictions19,46,47,49,59 and disordered porous materials.25,60,61 As the pressure increases slightly above P ) 0.46P0, partial irreversible condensation occurs within the pore as revealed by the jump in the adsorbed amount (see Figure 6). Only a much smaller gas bubble remains within the pore after this sudden uptake in Ar atoms (see the second snapshot in Figure 8). The filling of the remaining porosity consists of the deformation of the last gas bubble until its irreversible condensation occurs at a pressure P ) 0.52P0. As shown in Figure 9, the pore filling for Xe at 195 K is similar to that for Ar at 87 K. The fluid first condenses at P ∼ 0.47P0 in the center part of the pore, leading to the formation of a liquid bridge. The pore filling ends with the irreversible condensation of the gas nano-bubble at P ∼ 0.56P0. As mentioned above, pore emptying for model B depends on the adsorbate. Ar in this pore empties entirely at P ) 0.25P0 through cavitation, as in the case of pore model A. We already noted that this effect is due to the infinite length of the pores used in this work, which prevents evaporation to occur at equilibrium. In contrast, desorption for Xe confined at 195 K in pore model B is quasi-continuous (i.e., a combination of reversible and irreversible paths) and resembles the pore filling mechanism. These results show that the strength of the adsorbate/adsorbent interactions can have a significant effect on the way evaporation occurs. 3.3. Isosteric Heat of Adsorption. Ar and Xe isosteric heats of adsorption, Qst, are shown in Figure 10 for pore models A and B. These curves are typical of adsorption on heterogeneous surfaces as they decrease down to a plateau value that is close (56) Coasne, B.; Hung, F. R.; Siperstein, F. R.; Gubbins, K. E. Ann. Chim. Sci. Mater. 2005, 30, 375-383. (57) Bock, H.; Schoen, M. Phys. ReV. E 1999, 59, 4122-4136. (58) Bock, H.; Diestler, D. J.; Schoen, M. J. Phys.: Condens. Matter 2001, 13, 4697-4714. (59) Libby, B.; Monson, P. A. Langmuir 2004, 20, 4289-4294. (60) Detcheverry, F.; Kierlik, E.; Rosinberg, M. L.; Tarjus, G. Phys. ReV. E 2003, 68, 061504. (61) Detcheverry, F.; Kierlik, E.; Rosinberg, M. L.; Tarjus, G. Langmuir 2004, 20, 8006-8014.

Coasne et al.

Figure 9. Transverse section views of Xe adsorbed at 195 K in the MCM-41 pore model B (R ) 3.2 nm) at P ) 0.47P0 (top) and P ) 0.56P0 (bottom). Black spheres are hydrogen atoms at the pore surface and gray spheres are Xe atoms. The two extremities of the pore are connected to each other, due to the use of periodic boundary condition along the pore axis. The black lines delimit the regions of low density that are not filled (nano-bubbles).

Figure 10. (Top) Isosteric heat of adsorption of Ar at 87 K obtained by GCMC simulations: (circles) pore model A, (squares) pore model B. (closed triangles) are experiments for MCM-41 pores with D ) 4.5 nm by Neimark et al.62 Adsorbed amounts in abscissa have been normalized to the number of Ar atoms, N0, when the pore is completely filled. (Bottom) Isosteric heat of adsorption of Xe at 195 K obtained by GCMC simulations: (circles) pore model A, (squares) pore model B. (closed triangles) are experiments by Burgess et al.63 for silica glass Vycor having an average pore diameter of 4.5 nm. Adsorbed amounts in abscissa have been normalized to the number of Xe atoms, N0, when the pore is completely filled.

to the heat of liquefaction of bulk Ar, ∼6.5 kJ/mol, and bulk Xe, 12.5 kJ/mol. We also report in Figure 10 the experimental heat of adsorption curve obtained by Neimark et al.62 for Ar at 87 K in MCM-41 pores of diameter about 4.5 nm. Our results are also compared with the experimental data by Burgess et al.63 for Xe adsorption at 195 K in a silica Vycor glass having an average (62) Neimark, A. V.; Ravikovitch, P. I.; Gru¨n, M.; Schu¨th, F.; Unger, K. K. J. Colloid Interface Sci. 1998, 207, 159-169.

Adsorption of Gases in MCM-41 Materials

Langmuir, Vol. 22, No. 1, 2006 201

are smooth at the molecular scale but rough over the range 2050 and 80-250 Å. Despite the fact that pore model B was carved out of an atomistic block and that the structure was relaxed using an atomistic potential, the observed heterogeneity using Ar or Xe as a probe may still be influenced by the lattice structure used to generate the original material, because these probes are smaller than one lattice unit. It is worth noting that lattice simulations are useful when the properties of interest are at scales larger than one lattice spacing. On the other hand, lattice simulations cannot be expected to provide reliable information when looking at phenomena that occurs on length scales smaller than one lattice unit.

4. Conclusion

Figure 11. (Left) Fluid-wall (open symbols) and fluid-fluid (closed symbols) contributions to the isosteric heat of adsorption of Ar at 87 K in silica pores. Circles and squares are data for MCM-41 pore models A and B, respectively. (Bottom) Fluid-wall (open symbols) and fluid-fluid (closed symbols) contributions to the isosteric heat of adsorption of Xe at 195 K in silica pores. Circles and squares are data for MCM-41 pore models A and B, respectively.

pore diameter of 4.5 nm. Both the Ar and Xe isosteric heats for model A match better the experiments than those for model B. At low pore filling fractions, i.e., N/N0 < 0.2 for Ar and N/N0 < 0.3 for Xe, Qst for model B is significantly larger than that for model A. In agreement with previous works,43,52 the overestimated isosteric heat for the rough pore is consistent with the larger adsorbed amounts observed for this system compared to the smooth pore (Figure 6). The energetic differences between the two atomistic pore models can be understood by examining in Figure 11 the wall/ fluid and fluid/fluid contributions of the total isosteric heat of adsorption. The Ar/Ar and Xe/Xe contributions are very similar for both pore models A and B over the entire range of filling fractions. In contrast, the wall/Ar and wall/Xe contributions for model B are larger than those for model A; this difference, which is due to the surface roughness of model B, is responsible for the large isosteric heat of adsorption Qst for this pore (see Figure 10). The total isosteric heat of adsorption for model A reproduces reasonably the experimental curves (in the previous section, it was also found that model A matches better the experimental adsorption isotherm for Ar in MCM-41 pores). Comparison with the experimental data for Ar and Xe suggests that the surface disorder/roughness for model B is somewhat too great, i.e., contains too many high energy sites. The latter are of course related to the presence of single lattice unit cubes that are reminiscent of the initial lattice structure. On the other hand, we recall that the Porod exponent for model B, x ) 3.6, is much closer to the experimental values reported in the literature than that for model A, x ) 4.0. These results suggest that model B has approximately the correct surface disorder at length scales given by the Porod range 12-50 Å, whereas its atomistic surface roughness is too great. In contrast, it seems that model A has the correct atomistic roughness but does not reproduce the surface disorder of the experimental materials at larger length scales. This conclusion is compatible with the experimental investigation by Sonwane et al.,6 where it was shown that MCM-41 samples (63) Burgess, C. G. V.; Everett, D. H.; Nuttall, S. Langmuir 1990, 6, 17341738.

In this work, we report the preparation of two different atomistic models of MCM-41 silica mesopores having the same pore size. Model A consists of a regular cylindrical pore having a constant section of radius R0 ) 3.2 nm. Model B has an important surface disorder that reproduces the morphological features of an onlattice simulation mimicking the synthesis process of MCM-41 pores. Both models are obtained by carving the pore out of a silica block. The pore/void interface is then modeled in a realistic way by (i) saturation of the oxygen dangling bonds with hydrogen atoms and (ii) by amorphization of the surface using both randomization and relaxation techniques. Differences between the two pore models are discussed in terms of SANS spectra, adsorption isotherms, and isosteric heat curves. Many different features are found between the Ar or Xe adsorption isotherms for model A and those for model B. For instance, the adsorption branch for model B corresponds to a quasi-continuous pore filling involving coexistence within the pore of liquid bridges and gas nano-bubbles, while the discontinuous adsorption branch for model A conforms to the classical picture of capillary condensation in regular nanopores. It is also found that the features of the adsorption/desorption cycles for model B depend on the adsorbate, i.e., Ar or Xe. A large hysteresis loop with asymmetrical adsorption and desorption branches is observed in the case of Ar adsorption in model B; this shape departs from the symmetrical hysteresis loop obtained for the regular nanopore corresponding to model A. On the other hand, the adsorption/desorption isotherm for Xe in model B is a quasicontinuous and quasireversible mechanism, which is composed of both reversible/continuous and slightly irreversible/discontinuous paths. To interpret the differences between both pore models, we have examined for each system the isosteric heat of adsorption and its fluid/wall and fluid/fluid contributions. It has been found that, due to the atomic surface roughness for model B, the Ar and Xe isosteric heats for this sample are larger than those for model A. Comparison with experimental data shows that model A matches better the adsorption isotherm and isosteric heat curves than model B. However, it cannot be concluded that model A is closer to realistic MCM-41 samples than model B, since the Porod exponent for the latter model is much closer to the experimental value found by several authors. Such a result indicates that lattice simulations are able to reproduce properties on length scales that are larger than one lattice unit. Our conclusions can be summarized as follows. Model B seems to be too rough at the molecular scale, but reproduces reasonably the surface disorder of real MCM-41 at length scales given by the Porod range 12-50 Å. In contrast, model A is smooth at small length scales in agreement with experiments, but does not exhibit the appropriate surface disorder at larger length scales. We are currently trying to improve our model of MCM-41 porous materials by starting with model B obtained after the use of b-spline64 to model the pore surface in the lattice structure. The b-spline approach serves as a

202 Langmuir, Vol. 22, No. 1, 2006

reasonable interpolation to obtain a pore surface with subatomistic resolution < 0.1 nm; the control points in this b-spline model are the lattice surface coordinates, which allows us to retain the surface undulations and irregularities from the mimetic simulation.65 We should also consider the use of more efficient relaxation techniques to reconstruct the surface of the pore. This may allow us to obtain porous samples combining the correct atomistic roughness (similar to that for model A) and the correct surface disorder at larger length scales (similar to that for model B). Once a satisfactory model of MCM-41 materials is obtained, we will extend our work to realistic models of SBA-15 pores, i.e., silica pores which are linked through porous interconnections.11,66 Numerical models of these materials, which are obtained from on-lattice Monte Carlo simulations, are available.67 In future works, we will also investigate freezing of simple adsorbates such as krypton in these realistic models of mesoporous silica materials. Such a study will allow us to determine whether the disordered pore surface prevents crystallization of the confined fluid.68,69 Finally, we will also use our realistic models to determine the effect of the surface disorder on the adsorption and desorption (64) Piegl, L.; Tiller, W. The NURBS Book; Springer: New York, 1995. (65) Bhattacharya, S.; Coasne, B.; Hung, F. R.; Gubbins, K. E. In Proceedings of the 7th International Symposium on the Characterization of Porous Solids. Stud. Surf. Sci. Catal. 2005, accepted. (66) Galarneau, A.; Cambon, H.; Di Renzo, F.; Ryoo, R.; Choi, M.; Fajula, F. New J. Chem. 2003, 27, 73-79. (67) Bhattacharya, S.; Gubbins, K. E. J. Chem. Phys. 2005, 123, 134907.

Coasne et al.

scanning curves. Such curves, which are obtained by reversing upon adsorption or desorption the direction of change in the pressure, are of fundamental importance to test the validity of the independent domain theory for capillary condensation hysteresis.70,71 In parallel with experiments, this work would allow us to clarify the origin of capillary hysteresis loops that are usually observed for nanoporous materials. Acknowledgment. We thank the Department of Energy (Grant No. DE-FG02-98ER14847) and the Petroleum Research Fund of the American Chemical Society for support of this work. F.R.S. is grateful to the Ramo´n y Cajal program from the Ministerio de Ciencia y Tecnologı´a and the support of the project CTQ2004-03346/PPQ. This research was performed using supercomputing resources from San Diego Supercomputer Center (NSF/MRAC-CHE050047S), the National Energy Research Scientific Computing Center (DOE-DE-FGO2-98ER14847) and the High Performance Computing Center at North Carolina State University. LA051676G (68) Hung, F. R.; Coasne, B.; Gubbins, K. E.; Siperstein, F. R.; SliwinskaBartkowiak, M. In Proceedings of the 7th International Symposium on the Characterization of Porous Solids. Stud. Surf. Sci. Catal. 2005, accepted. (69) Hung, F. R.; Coasne, B.; Gubbins, K. E.; Sliwinska-Bartkowiak, M. to be published. (70) Everett, D. H. The Solid-Gas Interface; Flood, E. A., Ed.; Marcel Dekker: New York, 1967; Vol. 2, 1055-1113. (71) Coasne, B.; Gubbins, K. E.; Pellenq, R. J.-M. Phys. ReV. B 2005, 72, 024304.