- Email: [email protected]
Computer simulations of physical adsorption: a historical review
Applied Surface Science 196 (2002) 3–12
Computer simulations of physical adsorption: a historical review William Steele* Department of Chemistry, Pennsylvania State University, University Park, PA 16802, USA
Abstract Computer simulations of two-dimensional (2D) fluids are discussed in connection with their use in determining the isotherms and thermodynamic properties of monolayers adsorbed on model surfaces. The use of powerful algorithms such as the Grand Canonical Monte Carlo (GCMC) to simulate a greatly expanded range of adsorption systems is then discussed in connection with computations of multilayer adsorption, adsorption on heterogeneous surfaces and in pores. # 2002 Published by Elsevier Science B.V. Keywords: Physical adsorption; Computer simulation; Isotherms; Surface models; Porous solids; Monolayers; Grand Canonical Monte Carlo
1. Introduction The development of algorithms [1] that allow one to utilize digital computers in the numerical evaluation of the fundamental equations of statistical mechanics that provide the link between molecular and macroscopic behavior has proved to be an extremely powerful research technique. Both thermodynamic and dynamic properties of fluids and solids are being computed with results that are reported in a very large number of research papers. Early studies tended to use the agreement between experiment and simulation as a test of the accuracy of the model intermolecular interactions employed and, later, as a way to generate realistic equations of state and transport data under conditions that are difficult or impossible to duplicate in the laboratory. However, much of the current work in this field is now devoted to systems that are too complex to characterize without the insights provided * Tel.: þ1-814-865-3711. E-mail address: [email protected] (W. Steele).
by simulations of models of such cases. Many important physical adsorption systems fall into this category if one views the adsorbates as fluid or solid phases in confined geometries produced by the solid adsorbent. This review comprises an attempt to follow the development of this field, starting with the studies of twodimensional (2D) hard-disc fluids that were published in the 1960s [2] and continuing with the development and application of powerful algorithms [3] for the computation of the thermodynamic and structural properties of films adsorbed on and in inert solid adsorbents [4].
2. In the beginning Roughly 40 years ago, the availability of digital computers for scientific research began to be adequate for the solution of many hitherto intractable problems. Among these was the calculation of the equation of state of simple fluids over a range of density up to and beyond that for the freezing point. The simplest
0169-4332/02/$ – see front matter # 2002 Published by Elsevier Science B.V. PII: S 0 1 6 9 - 4 3 3 2 ( 0 2 ) 0 0 0 3 8 - 7
4
W. Steele / Applied Surface Science 196 (2002) 3–12
system was that of hard spheres and, naturally, this was studied first as an example of numerical statistical mechanics. In addition to the problem of availability of adequate computer resources (which was impressively solved by the rapid advances in digital technology), several developments were needed for the successful simulation of model fluids. These included algorithms that would generate Gibbsian ensembles such as the canonical or grand canonical starting from a very small number of particles (by macroscopic standards) in a small box of volume that produced a fluid with a density of physical interest. The devices that allowed one to simulate macroscopic ensembles of interacting particles were and still are those of periodic boundary conditions and the minimum image convention. Briefly, systems with periodic boundaries were constructed by surrounding the initial small system with images of itself and allowing particles to pass in and out so that the influence of the walls is not important even though the number of particles remains invariant to passage out or in because an image particle passes in or out whenever passage through a wall occurs. This device works best when the system size is large enough to have the correlations between a particle and its nearest images be unimportant—which can be a problem when solid phases are studied. The minimum image convention requires that the total interaction of a particle with its neighbors be the sum of interactions with the actual neighbors plus those with the nearby images, especially when the particle of interest is near a wall or walls. Simulations of 2D fluids are facilitated by the small number of particles needed to give a comfortable system size. The calculation then requires that one evaluate many trajectories (molecular dynamics, MD) or random displacements (Monte Carlo, MC). For hard discs, an acceptable configuration requires that no particles overlap and that the number of trial configurations be large. The extension of this idea to non-circular particles (hard spherocylinders, chains of hard disks, etc.) is straightforward. The thermodynamic properties of ensembles of such systems can be obtained from the Gibbs isotherm, which relates the 2D pressure j at a given density G (particles per unit area) to m2D, the chemical potential of the 2D fluid [5]. However, the relationship to physical adsorption of a single layer of particles can be written as m2D ¼ mads mid , where mid is the chemical potential
of the adsorbed gas in its ideal or Henry’s law regime. One has Z x mads mid ad dj dx (1) ¼ 1 dx x kT kT 0 hp i mads mid ¼ ln KH G kT
(2)
where ad is the area per disc and KH the Henry’s law constant for the adsorption of the gas at pressure p (now in 3D) that is given by the well-known expression: Z A u KH ¼ ðe gsðzÞ=kT 1Þ dz (3) kT where ugs(z) is the gas–solid energy (i.e., the holding potential for the gas atoms on the surface) and A the surface area. The virial theorem [6] gives a generally valid expression for the 2D pressure j of a hard-disc fluid [7] that can be written as: j p ¼ 1 þ Gs2 ghd ðsÞ GkT 2
(4)
where s is the hard disc diameter and ghd(s) is the pair correlation function of the hard-disc fluid, extrapolated to the distance s where a pair of discs is in contact. This correlation function characterizes the structure of the fluid (or solid) and thus is an important factor that can be evaluated from the simulated configurations of the hard-disc fluid. Fig. 1 shows simulated 2D hard disc correlation functions for several densities [8]. In essence, g(r) shows the variation of local density in the neighborhood of a disc, normalized to 1 at large separation distances. Fig. 2 shows hard disc isotherms, simulated by MC and calculated from several theories [5]. It is evident that the scaled particle theory is the most accurate, although the isotherms calculated from the first six virial coefficients are good except at the highest densities shown. At high density, the widely used van der Waals’ isotherm is actually the worst of the theories shown here while the free area theory is the worst at low density. For hard bodies that are not circular, an extended form of the theory has been derived [9], exact when the particles have convex shapes and reasonably accurate for particles made up of overlapping discs (e.g. hard
W. Steele / Applied Surface Science 196 (2002) 3–12
5
Fig. 1. Simulated hard disc correlation functions for several densities. In essence, g(r) shows the local density in the neighborhood of a disc, normalized to 1 at large separation. The reduced densities G ¼ G=Gcp are 0.4 (solid line), 0.5 (dashed line) and 0.6 (short dashes), where Gcp is the density of the close-packed layer and r ¼ r=s.
diatomics or hard chains). One first replaces ghd(s) by gav, the pair correlation averaged over the orientations of a pair of particles in contact and then the contact length ps by S, the perimeter of a particle. Simulation studies have been reported for a number of 2D hard-body models, including diatomics (two overlapping circles) [10–12], triatomics (three overlapping circles with bond angles varying from 608 to 1808) [13] cyclic pentamers [14], disco-rectangles
[15,16], ellipses [17], hard rods [18] and oblate spheroids [19] (not strictly 2D because the particles are allowed to tilt out of the plane). Hard body systems are of interest by themselves, particularly for non-circular cases where the simulated configurations show orientational ordering similar to that found in 3D liquid crystals. As an illustration, snapshots of configurations for a fluid of hard ellipses with an axial ratio of 6 are shown in Fig. 3 for three
Fig. 2. Isotherms simulated for the hard disc system are compared with several theoretical expressions here. The dependence of reduced density x ¼ Gs2 upon the logarithm of the pressure times the Henry’s law constant is shown. (Note that the close-packed value for x is p/2 31=2 ¼ 0:907.).
6
W. Steele / Applied Surface Science 196 (2002) 3–12
Fig. 3. Snapshots of configurations of hard ellipses with an axial ratio of 6 are shown here for three values of Gpab corresponding to an isotropic phase (¼0.33), a nematic phase (¼0.60) and a solid phase (¼0.81). Here a and b are the axes of the ellipse.
densities [16]. Isotropic, nematic and solid phases are shown. Another possible use for 2D hard body simulations is as the basis for the perturbation theory of fluids that has been shown to give an excellent theory of 3D liquids. For example, monolayers of particles with Lennard-Jones interactions can be described as having the configurations (and thus, the pair correlation functions) of a hard-disc fluid with an appropriately chosen size parameter [20,21]. The attractive part of the Lennard-Jones interaction energies can then be evaluated for the 2D fluids by approximating the actual pair correlations by the hard particle correlation functions similar to now standard theories of the 3D liquid [22]. Note that the configurations of pure hard body systems are independent of temperature, whereas softer interactions (Lennard-Jones, multipolar, etc.) are significantly dependent upon T, so that the addition of a soft inter-
action as a perturbation to the hard-body fluid adds to the length of the calculation. Other variants of the hard disc model can be constructed by adding interactions such as those for dipole moments. (This model could be thought of as a simple representation of fluids such as the hydrogen halides where the dipole moments dominate the angle-dependence of the intermolecular interactions.) In fact, this system [23] is quasi-2D because the dipoles are allowed to orient out of plane. However, the temperature of the simulation shown in Fig. 4 is rather low ðT ¼ kTs3 =m3 ¼ 1=9Þ so that the interacting dipoles tend to form a 2D layer. Snapshots of configurations of 200 discs of this kind are shown here for two densities. In both cases, the particles form long chains of head-to-tail dipoles, but at the higher density, these chains begin to pack into close-packed arrays.
Fig. 4. Snapshots of the 2D configurations of hard disc þ dipole particles at densities r ¼ rs2 ¼ 0:5 (left) and 0.8 (right). The dipole orientations are shown by the lines through each disc.
W. Steele / Applied Surface Science 196 (2002) 3–12
3. The great step forward Although the work on 2D systems is described here as pioneering, it should be noted that much of the work on non-circular particles in 2D has been published rather recently. In face, the distinguishing feature of the studies discussed in the previous section is that they employ the canonical ensemble for their MC simulations and this is not really convenient for generating adsorption isotherms. Use of the hard particle results in the construction of perturbation theories for interacting soft particles was an interesting idea that was made obsolete by the widespread appearance of direct simulations of the 3D adsorption of non-hard particles such as those with Lennard-Jones interactions. Consequently, simulation algorithms came into wide use which allow for multilayer formation and for gas–solid potentials ugs(s, z) that vary with s, the particle position over a surface with an atomic structure that gives rise to irregularities or to periodic variations in the holding or gas–solid potential. Of course, simulations of such systems are of much more interest than the 2D models discussed to this point. The specific algorithm that is the basis for much of the modern work is the Grand Canonical Monte Carlo (GCMC). In it, ensembles are computer generated for systems of molecules at constant chemical potential, volume and temperature. Since adsorption is generally studied by measuring the pressure of the gas in equilibrium with the adsorbed particles at some fixed T, the chemical potential is obtainable from that of the (nearly) ideal gas. The system volume is that of the void space in which adsorption occurs and thus is constant. The GCMC thus has the number and potential energy of the (adsorbed) particles as fluctuating quantities whose average values are generated by evaluating long sequences of random displacements and, in addition, by varying the number of adsorbed particles by random deletion and insertion. The details of this procedure have been discussed in detail elsewhere [3,4] and need not be repeated here. It is sufficient to say that rules for the acceptance of the random changes are applied that cause the system to approach the statistical mechanical ensemble proposed by J. Willard Gibbs. Note that the presence of an adsorbing surface requires that the periodic boundary conditions in 3D are no longer practical. For the simple case of a planar or nearly planar adsorbing surface, periodic boundaries
7
make sense only in the two dimensions parallel to the surface. A boundary parallel to the adsorbing surface is then necessary, thus creating a parallel-walled slit pore. The effects of this boundary are minimized by assuming a gas-boundary wall interaction that is short-ranged (i.e., a hard wall) and setting the distance between the adsorbing and boundary walls to be as large as possible. Of course appropriate changes in the periodic boundaries are needed for systems such as cylindrical pores, etc. Simulations in such systems result in calculations of hNi and hUi, the average number and potential energy of the adsorbed (and non-adsorbed) particles in the adsorption volume V at a fixed temperature T and pressure p. In addition, the microscopic number and energy densities can be evaluated to give detailed information concerning the nature of the adsorption process as a function of pressure and temperature. Only the interaction energy of adsorbing particles with each other and with the (inert) solid adsorbent need be known to allow one to generate data that compares with experiment. At present, these interactions are not as well known as one might like. Nevertheless, a very large number of studies employing the GCMC and its variants (constant P rather than constant V, etc. [4]) have been published. This work has been reviewed elsewhere [24] so a lengthy repetition here is unnecessary. The results of an early GCMC simulation of the adsorption isotherm of argon on graphitized carbon black are compared with experiments in Fig. 5. The Lennard-Jones Ar–Ar interaction used in this study was rather well-known from work on the 3D gas. (Alterations in this energy due to the nearby presence of the solid adsorbent are present, but are often ignored.) The Ar–solid interaction was semi-empirical, with constants that are refined by comparison with the readily calculable Henry’s law constant. The gas– solid interactions for non-reactive gases on insulating solids are generally assumed to be sums over the atoms in the solid. (A different approach is needed for metallic substrates where the conduction electrons play a large role.) The evaluation of these sums requires a knowledge of the atom positions—straightforward for crystals if the exposed lattice plane is known, but not for amorphous or partially amorphous materials. In the case of graphitized carbon black, basal planes are exposed, for the most part, and the structure is well-known. If one takes the Ar–C atom potential to be a Lennard-Jones function with constants
8
W. Steele / Applied Surface Science 196 (2002) 3–12
Fig. 5. Isotherms for Ar on graphitized carbon black simulated at two reduced temperatures T ¼ kT=egg (egg is the Ar–Ar interaction well depthÞ ¼ 0:67 (squares) and 1.00. These are compared with an experimental isotherm measured at 78 K. y is the coverage in layers and p0 is the bulk vapor pressure.
(well-depth and size) that gives Henry’s law adsorption in agreement with high T, low coverage measurements, one obtains a simulated multilayer isotherm that agrees reasonably well with the low T data as indicated in Fig. 5. The rounded steps that are found are due to layer-by-layer adsorption at finite T on uniform graphite basal planes; the amount of rounding is most likely enhanced by edge effects associated with the finite extent of these planes in the experimental sample. If one accepts the site-wise summation models as adequate starting points for the moment, it is clear that both the location of these sites and the gas–site interactions are needed for an adsorption simulation to be realistic. One of the very large number of published adsorption simulations will be discussed here to illustrate how such computations can be used to gain a
deeper knowledge of this complicated process. The particular system was basically that of krypton in several straight-walled cylindrical pores [25]. The pore was constructed by first simulating a random close-packed cubic box of hard spheres [26]. Deleting the spheres from the desired cylindrical volume then produced a cylindrical cavity. Periodic boundary conditions were imposed on the ends of the cylinder and the system was transformed into an adsorbing medium by assuming Lennard-Jones interactions between each sphere and the atoms of the Kr adsorbate. Fig. 6 gives a view of the insides of two such cylinders. The first point to note is that the atom–atom interaction model will generate a heterogeneous adsorbing surface. Real adsorbents that approximate such a model are the MCM 21 zeolites that have been widely used for experimental adsorption studies over the past few years [27]. (Actually, these materials possess hexagonal pores rather than cylindrical, but this distinction is generally neglected. Also, their pore radii are roughly twice as large as those used in these simulations.) It is important to realize that the heterogeneity of such surfaces can be quite significant, depending upon the relative sizes of the adsorbate atom and the atoms of the silica adsorbents. To illustrate this point, adsorption energies defined as the minima in the gas–solid interaction potentials for Kr atoms at various positions on the walls were calculated for a net of about 3000 equally spaced points. Model walls with Lennard˚— Jones gas–solid atom sizes of 1.84, 2.30, 3.05 A with some deleted wall atoms—denoted by sc for ˚ without and with a larger small corrugation and 3.15 A
Fig. 6. Views down the axis of two cylindrical pores with walls ˚ made up of atoms. On the left, the atoms are small (1.84 A ˚ ); diameter) compared to the Kr atoms shown in the center (3.78 A ˚ ) and the wall has been on the right the atoms are larger (3.14 A made rougher by randomly deleting one-half of the wall atoms in the innermost two layers.
W. Steele / Applied Surface Science 196 (2002) 3–12
number of deleted wall atoms denoted by lc. The minima of the interaction energy and the positions of these minima have the expected dependence upon the Kr position. For instance, a Kr atom directly over a solid atom gives a smaller minimum at a location closer to the pore axis than one which is over the midpoint of a pair or triplet of wall atoms because such an adsorbed atom will sink into the center of the minicluster and ‘‘touch’’ two or three solid atoms rather than the single wall atom of an atop position. These
9
variations have been quantified by constructing histograms of the distributions of adsorption energies U and of the distances of these minima from the pore axis Rm. These histograms are shown in Fig. 7 for pores that have been modeled such that the average Rm is ˚ and the average Um is 8.8 kJ/mol for all 10 A systems. It is clear that these pore walls are significantly heterogeneous and that the heterogeneity increases as the wall atom/adsorbate size ratio increases or as the wall is roughened by random
Fig. 7. Histograms for the distributions of pore radius Rm (top) as defined by the positions of the adsorption energy minima U and for the values of these minima (bottom) are shown for all pores considered except that with the smallest wall atom size (and the least heterogeneity). The curves for wall atom size 3.15 are omitted here.
10
W. Steele / Applied Surface Science 196 (2002) 3–12
deletion of some of the adsorbing atoms in the wall. Furthermore, the physical smoothness as measured by R becomes significantly less as the wall atom size and/ or the corrugation produced by missing atoms increases. These model pores were used for simulations of the adsorption of Kr at 120 K. The isotherms and the gas–solid contribution to the adsorption ener-
gies at submonolayer coverages are shown in Fig. 8. Monolayer coverages for all systems except the large corrugation (lc) case are estimated to be 33 atoms/nm. In the lc pore, the large number of ‘‘holes’’ in the walls allows for significant adsorption but it is hard to define the ‘‘monolayer’’ for such a rough surface. It is clear in Fig. 8 (top) that a full (lc) pore contains more adsorbed
Fig. 8. Simulated isotherms and the gas–solid parts of the adsorption energy are shown on the top and the bottom, respectively. The isotherms were evaluated for adsorption up to pore-filling and both the adsorption and desorption branches of the curves are shown. The simulated energies are shown only for coverages up to the condensation points, which are slightly greater than monolayer.
W. Steele / Applied Surface Science 196 (2002) 3–12
atoms per unit length than those with smoother walls even when the average radius Rm is fixed. For submonolayer adsorption, the distributions of adsorption energies shown in Fig. 7 give rise to a nearly coverageindependent average gas–solid energy for the smallest wall atoms with the narrowest distributions of U. Much larger changes in this energy can be seen for larger wall atoms in pores both with and without imposed corrugation. (Note that the solid with size ˚ has adsorption properties almost parameter 3.15 A ˚ size parameter and identical to those with a 3.05 A a small additional wall corrugation.) The shapes of the isotherms at low coverages change significantly as the heterogeneity increases, going from curves with gradual almost linear rises for the surfaces lacking in strong sites to those with the rapid initial increases that are typically found for surfaces with strong gas–solid interactions. (Here, the use of the word site does not mean a preferred specific point for adsorption, but more generally, one of the 3000 points where the energies were calculated for use in the histograms.) Based on our understanding of the role of heterogeneity in physical adsorption [28], these changes in isotherm shape with changes in surface heterogeneity are as expected. However, it would be interesting to see if the adsorption energy distributions in the histograms of Fig. 7 could be used in the standard theories of adsorption on heterogeneous surfaces [27] to repro-
11
duce at least the monolayer portions of the isotherms shown here. Capillary condensation and hysteresis occur in all simulated isotherms shown here, with complete tube filling at 70–75 atoms/nm for all systems except the (lc) case. The positions and widths of the hysteresis loops vary somewhat with the changing surface heterogeneity. In spite of the presence of significant surface heterogeneity, the steps in these loops are essentially vertical on both the adsorption and desorption branches. Theory is unable to quantitatively explain this hysteresis at present, although it is likely that the considerable efforts being expended on a solution to this problem will give useful results soon. Finally, the dependence of the adsorbed phase densities upon distance from the pore axis is shown for coverages near monolayer completion in Fig. 9. Only for the two smoothest surfaces does one see anything like a single layer of atoms adsorbed near the wall. As the roughness increases, the atoms in the ‘‘monolayer’’ spread out over a range of distances, thus making the operational definition of monolayer completion more and more arbitrary and of less significance. We have devoted considerable space to discussion of the use of computer simulation not just to produce thermodynamic observables from a microscopic model but, in addition, to generate detailed data about
Fig. 9. The density of adsorbed atoms is shown as a function of distance from the pore axis for coverages close to monolayer completion. The formation of a layer at Ravg is expected for atoms adsorbed on a smooth wall.
12
W. Steele / Applied Surface Science 196 (2002) 3–12
the model surface properties that can enable critical tests of theory. It would seem that analyses which utilize the full power of the computer in this and other ways not discussed here will become an important aspect of our progress towards a full understanding of physical adsorption on or in complex solids.
[18] [19] [20] [21]
References [1] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, E. Teller, J. Chem. Phys. 21 (1953) 1087. [2] B.J. Alder, T.E. Wainwright, Phys. Rev. 127 (1962).; B.J. Alder, T.E. Wainwright, J. Chem. Phys. 33 (1960) 1439. [3] M.P. Allen, D.J. Tildesly, Computer Simulation of Liquids, Oxford Science Publications, 1987.; D. Frenkel, B. Smit, Understanding Molecular Simulations: From Algorithms to Applications, Academic Press, New York, 1996. [4] D. Nicholson, N.G. Parsonage, Computer Simulation and the Statistical Mechanics of Adsorption, Academic Press, New York, 1982. [5] W.A. Steele, The Interaction of Gases with Solid Surfaces, Pergamon Press, Oxford, 1974. [6] T.L. Hill, Statistical Mechanics, McGraw-Hill, New York, 1956. [7] T. Boublik, Mol. Phys. 29 (1975) 421. [8] D.G. Chae, F.H. Ree, T. Ree, J. Chem. Phys. 50 (1969) 1581. [9] M.J. Maeso, J.R. Solana, J. Chem. Phys. 102 (1995) 8562. [10] T. Boublik, Mol. Phys. 63 (1988) 685. [11] J. Talbot, D.J. Tildesley, J. Chem. Phys. 83 (1985) 6419. [12] C.W. Goulding, M. Rigby, Mol. Phys. 75 (1992) 623. [13] A. Bran˜ ka, K.W. Wojciechowski, Mol. Phys. 78 (1993) 1513. [14] I. Nezbeda, T. Boublik, O. Trnka, Czech. J. Phys. B 25 (1975) 119. [15] J. Viellard-Baron, J. Chem. Phys. 56 (1972) 4729. [16] J.A. Cuesta, D. Frenkel, Phys. Rev. A 42 (1990) 2126. [17] H. Schlacken, H.J. Mo¨ gel, P. Schiller, Mol. Phys. 93 (1998) 777.
[22] [23] [24]
[25]
[26]
[27]
[28]
D. Frenkel, R. Eppenga, Phys. Rev. A 31 (1985) 1776. P. Siders, Mol. Phys. 68 (1989) 1001. W.A. Steele, J. Chem. Phys. 65 (1976) 5256. A. Mulero, F. Cuadros, Langmuir 12 (1996) 3265; A. Mulero, F. Cuadros, Chem. Phys. 205 (1996) 379; A. Mulero, F. Cuadros, Chem. Phys. 156 (1991) 33; A. Mulero, F. Cuadros, Chem. Phys. 159 (1992) 89; A. Mulero, F. Cuadros, Chem. Phys. 160 (1992) 375; A. Mulero, F. Cuadros, Chem. Phys. 160 (1993) 375; A. Mulero, F. Cuadros, J. Phys. Chem. 99 (1995) 419; F. Cuadros, A. Mulero, L. Morala, V. Gomez-Serrano, Langmuir 17 (2001) 1576. J.D. Weeks, D. Chandler, H.C. Andersen, J. Chem. Phys. 54 (1971) 5237. J.J. Weis, Mol. Phys. 93 (1998) 361. W.A. Steele, in: J. Schwartz, C. Contescue (Eds.), Surfaces of Nanoparticles and Porous Materials, Marcel Dekker, New York, 1999. W.A. Steele, M.J. Bojan, Adv. Colloid Interf. Sci. 76–77 (1998) 153; B. McEnany, T.J. Mays, J. Rouque´ rol, F. Rodri uez-Reinoso, K.S.W. Sing, K.K. Unger (Eds.), Characterization of Porous Solids IV, Royal Society of Chemistry, 1997, p. 49.; M.D. LeVan (Ed.), Fundamentals of Adsorption, Kluwer Academic Publishers, Dordrecht, 1996. V.A. Bakaev, Surf. Sci. 215 (1988) 571; V.A. Bakaev, W.A. Steele, Langmuir 8 (1992) 148; A. Bonissent, B. Mutaftschieve, Phil. Mag. 35 (1977) 65. M. Grun, K.K. Unger, A. Matsumoto, K. Tsutsumi, Characterization of Porous Solids IV, Royal Society of Chemistry, 1997, p. 8.; P.J. Branton, P.G. Hall, K.S.W. Sing, Characterization of Porous Solids IV, Royal Society of Chemistry, 1997, p. 98.; P.I. Ravikovitch, G.L. Haller, A.V. Neimark, Adv. Colloid Interf. Sci. 76–77 (1998) 203. M. Jaroniec, R. Madey, Physical Adsorption on Heterogeneous Solids, Elsevier, Amsterdam, 1998.; W. Rudzin´ ski, D.H. Everett, Adsorption of Gases on Heterogeneous Surfaces, Academic Press, New York, 1992.
Computer simulations of physical adsorption: a historical review William Steele* Department of Chemistry, Pennsylvania State University, University Park, PA 16802, USA
Abstract Computer simulations of two-dimensional (2D) fluids are discussed in connection with their use in determining the isotherms and thermodynamic properties of monolayers adsorbed on model surfaces. The use of powerful algorithms such as the Grand Canonical Monte Carlo (GCMC) to simulate a greatly expanded range of adsorption systems is then discussed in connection with computations of multilayer adsorption, adsorption on heterogeneous surfaces and in pores. # 2002 Published by Elsevier Science B.V. Keywords: Physical adsorption; Computer simulation; Isotherms; Surface models; Porous solids; Monolayers; Grand Canonical Monte Carlo
1. Introduction The development of algorithms [1] that allow one to utilize digital computers in the numerical evaluation of the fundamental equations of statistical mechanics that provide the link between molecular and macroscopic behavior has proved to be an extremely powerful research technique. Both thermodynamic and dynamic properties of fluids and solids are being computed with results that are reported in a very large number of research papers. Early studies tended to use the agreement between experiment and simulation as a test of the accuracy of the model intermolecular interactions employed and, later, as a way to generate realistic equations of state and transport data under conditions that are difficult or impossible to duplicate in the laboratory. However, much of the current work in this field is now devoted to systems that are too complex to characterize without the insights provided * Tel.: þ1-814-865-3711. E-mail address: [email protected] (W. Steele).
by simulations of models of such cases. Many important physical adsorption systems fall into this category if one views the adsorbates as fluid or solid phases in confined geometries produced by the solid adsorbent. This review comprises an attempt to follow the development of this field, starting with the studies of twodimensional (2D) hard-disc fluids that were published in the 1960s [2] and continuing with the development and application of powerful algorithms [3] for the computation of the thermodynamic and structural properties of films adsorbed on and in inert solid adsorbents [4].
2. In the beginning Roughly 40 years ago, the availability of digital computers for scientific research began to be adequate for the solution of many hitherto intractable problems. Among these was the calculation of the equation of state of simple fluids over a range of density up to and beyond that for the freezing point. The simplest
0169-4332/02/$ – see front matter # 2002 Published by Elsevier Science B.V. PII: S 0 1 6 9 - 4 3 3 2 ( 0 2 ) 0 0 0 3 8 - 7
4
W. Steele / Applied Surface Science 196 (2002) 3–12
system was that of hard spheres and, naturally, this was studied first as an example of numerical statistical mechanics. In addition to the problem of availability of adequate computer resources (which was impressively solved by the rapid advances in digital technology), several developments were needed for the successful simulation of model fluids. These included algorithms that would generate Gibbsian ensembles such as the canonical or grand canonical starting from a very small number of particles (by macroscopic standards) in a small box of volume that produced a fluid with a density of physical interest. The devices that allowed one to simulate macroscopic ensembles of interacting particles were and still are those of periodic boundary conditions and the minimum image convention. Briefly, systems with periodic boundaries were constructed by surrounding the initial small system with images of itself and allowing particles to pass in and out so that the influence of the walls is not important even though the number of particles remains invariant to passage out or in because an image particle passes in or out whenever passage through a wall occurs. This device works best when the system size is large enough to have the correlations between a particle and its nearest images be unimportant—which can be a problem when solid phases are studied. The minimum image convention requires that the total interaction of a particle with its neighbors be the sum of interactions with the actual neighbors plus those with the nearby images, especially when the particle of interest is near a wall or walls. Simulations of 2D fluids are facilitated by the small number of particles needed to give a comfortable system size. The calculation then requires that one evaluate many trajectories (molecular dynamics, MD) or random displacements (Monte Carlo, MC). For hard discs, an acceptable configuration requires that no particles overlap and that the number of trial configurations be large. The extension of this idea to non-circular particles (hard spherocylinders, chains of hard disks, etc.) is straightforward. The thermodynamic properties of ensembles of such systems can be obtained from the Gibbs isotherm, which relates the 2D pressure j at a given density G (particles per unit area) to m2D, the chemical potential of the 2D fluid [5]. However, the relationship to physical adsorption of a single layer of particles can be written as m2D ¼ mads mid , where mid is the chemical potential
of the adsorbed gas in its ideal or Henry’s law regime. One has Z x mads mid ad dj dx (1) ¼ 1 dx x kT kT 0 hp i mads mid ¼ ln KH G kT
(2)
where ad is the area per disc and KH the Henry’s law constant for the adsorption of the gas at pressure p (now in 3D) that is given by the well-known expression: Z A u KH ¼ ðe gsðzÞ=kT 1Þ dz (3) kT where ugs(z) is the gas–solid energy (i.e., the holding potential for the gas atoms on the surface) and A the surface area. The virial theorem [6] gives a generally valid expression for the 2D pressure j of a hard-disc fluid [7] that can be written as: j p ¼ 1 þ Gs2 ghd ðsÞ GkT 2
(4)
where s is the hard disc diameter and ghd(s) is the pair correlation function of the hard-disc fluid, extrapolated to the distance s where a pair of discs is in contact. This correlation function characterizes the structure of the fluid (or solid) and thus is an important factor that can be evaluated from the simulated configurations of the hard-disc fluid. Fig. 1 shows simulated 2D hard disc correlation functions for several densities [8]. In essence, g(r) shows the variation of local density in the neighborhood of a disc, normalized to 1 at large separation distances. Fig. 2 shows hard disc isotherms, simulated by MC and calculated from several theories [5]. It is evident that the scaled particle theory is the most accurate, although the isotherms calculated from the first six virial coefficients are good except at the highest densities shown. At high density, the widely used van der Waals’ isotherm is actually the worst of the theories shown here while the free area theory is the worst at low density. For hard bodies that are not circular, an extended form of the theory has been derived [9], exact when the particles have convex shapes and reasonably accurate for particles made up of overlapping discs (e.g. hard
W. Steele / Applied Surface Science 196 (2002) 3–12
5
Fig. 1. Simulated hard disc correlation functions for several densities. In essence, g(r) shows the local density in the neighborhood of a disc, normalized to 1 at large separation. The reduced densities G ¼ G=Gcp are 0.4 (solid line), 0.5 (dashed line) and 0.6 (short dashes), where Gcp is the density of the close-packed layer and r ¼ r=s.
diatomics or hard chains). One first replaces ghd(s) by gav, the pair correlation averaged over the orientations of a pair of particles in contact and then the contact length ps by S, the perimeter of a particle. Simulation studies have been reported for a number of 2D hard-body models, including diatomics (two overlapping circles) [10–12], triatomics (three overlapping circles with bond angles varying from 608 to 1808) [13] cyclic pentamers [14], disco-rectangles
[15,16], ellipses [17], hard rods [18] and oblate spheroids [19] (not strictly 2D because the particles are allowed to tilt out of the plane). Hard body systems are of interest by themselves, particularly for non-circular cases where the simulated configurations show orientational ordering similar to that found in 3D liquid crystals. As an illustration, snapshots of configurations for a fluid of hard ellipses with an axial ratio of 6 are shown in Fig. 3 for three
Fig. 2. Isotherms simulated for the hard disc system are compared with several theoretical expressions here. The dependence of reduced density x ¼ Gs2 upon the logarithm of the pressure times the Henry’s law constant is shown. (Note that the close-packed value for x is p/2 31=2 ¼ 0:907.).
6
W. Steele / Applied Surface Science 196 (2002) 3–12
Fig. 3. Snapshots of configurations of hard ellipses with an axial ratio of 6 are shown here for three values of Gpab corresponding to an isotropic phase (¼0.33), a nematic phase (¼0.60) and a solid phase (¼0.81). Here a and b are the axes of the ellipse.
densities [16]. Isotropic, nematic and solid phases are shown. Another possible use for 2D hard body simulations is as the basis for the perturbation theory of fluids that has been shown to give an excellent theory of 3D liquids. For example, monolayers of particles with Lennard-Jones interactions can be described as having the configurations (and thus, the pair correlation functions) of a hard-disc fluid with an appropriately chosen size parameter [20,21]. The attractive part of the Lennard-Jones interaction energies can then be evaluated for the 2D fluids by approximating the actual pair correlations by the hard particle correlation functions similar to now standard theories of the 3D liquid [22]. Note that the configurations of pure hard body systems are independent of temperature, whereas softer interactions (Lennard-Jones, multipolar, etc.) are significantly dependent upon T, so that the addition of a soft inter-
action as a perturbation to the hard-body fluid adds to the length of the calculation. Other variants of the hard disc model can be constructed by adding interactions such as those for dipole moments. (This model could be thought of as a simple representation of fluids such as the hydrogen halides where the dipole moments dominate the angle-dependence of the intermolecular interactions.) In fact, this system [23] is quasi-2D because the dipoles are allowed to orient out of plane. However, the temperature of the simulation shown in Fig. 4 is rather low ðT ¼ kTs3 =m3 ¼ 1=9Þ so that the interacting dipoles tend to form a 2D layer. Snapshots of configurations of 200 discs of this kind are shown here for two densities. In both cases, the particles form long chains of head-to-tail dipoles, but at the higher density, these chains begin to pack into close-packed arrays.
Fig. 4. Snapshots of the 2D configurations of hard disc þ dipole particles at densities r ¼ rs2 ¼ 0:5 (left) and 0.8 (right). The dipole orientations are shown by the lines through each disc.
W. Steele / Applied Surface Science 196 (2002) 3–12
3. The great step forward Although the work on 2D systems is described here as pioneering, it should be noted that much of the work on non-circular particles in 2D has been published rather recently. In face, the distinguishing feature of the studies discussed in the previous section is that they employ the canonical ensemble for their MC simulations and this is not really convenient for generating adsorption isotherms. Use of the hard particle results in the construction of perturbation theories for interacting soft particles was an interesting idea that was made obsolete by the widespread appearance of direct simulations of the 3D adsorption of non-hard particles such as those with Lennard-Jones interactions. Consequently, simulation algorithms came into wide use which allow for multilayer formation and for gas–solid potentials ugs(s, z) that vary with s, the particle position over a surface with an atomic structure that gives rise to irregularities or to periodic variations in the holding or gas–solid potential. Of course, simulations of such systems are of much more interest than the 2D models discussed to this point. The specific algorithm that is the basis for much of the modern work is the Grand Canonical Monte Carlo (GCMC). In it, ensembles are computer generated for systems of molecules at constant chemical potential, volume and temperature. Since adsorption is generally studied by measuring the pressure of the gas in equilibrium with the adsorbed particles at some fixed T, the chemical potential is obtainable from that of the (nearly) ideal gas. The system volume is that of the void space in which adsorption occurs and thus is constant. The GCMC thus has the number and potential energy of the (adsorbed) particles as fluctuating quantities whose average values are generated by evaluating long sequences of random displacements and, in addition, by varying the number of adsorbed particles by random deletion and insertion. The details of this procedure have been discussed in detail elsewhere [3,4] and need not be repeated here. It is sufficient to say that rules for the acceptance of the random changes are applied that cause the system to approach the statistical mechanical ensemble proposed by J. Willard Gibbs. Note that the presence of an adsorbing surface requires that the periodic boundary conditions in 3D are no longer practical. For the simple case of a planar or nearly planar adsorbing surface, periodic boundaries
7
make sense only in the two dimensions parallel to the surface. A boundary parallel to the adsorbing surface is then necessary, thus creating a parallel-walled slit pore. The effects of this boundary are minimized by assuming a gas-boundary wall interaction that is short-ranged (i.e., a hard wall) and setting the distance between the adsorbing and boundary walls to be as large as possible. Of course appropriate changes in the periodic boundaries are needed for systems such as cylindrical pores, etc. Simulations in such systems result in calculations of hNi and hUi, the average number and potential energy of the adsorbed (and non-adsorbed) particles in the adsorption volume V at a fixed temperature T and pressure p. In addition, the microscopic number and energy densities can be evaluated to give detailed information concerning the nature of the adsorption process as a function of pressure and temperature. Only the interaction energy of adsorbing particles with each other and with the (inert) solid adsorbent need be known to allow one to generate data that compares with experiment. At present, these interactions are not as well known as one might like. Nevertheless, a very large number of studies employing the GCMC and its variants (constant P rather than constant V, etc. [4]) have been published. This work has been reviewed elsewhere [24] so a lengthy repetition here is unnecessary. The results of an early GCMC simulation of the adsorption isotherm of argon on graphitized carbon black are compared with experiments in Fig. 5. The Lennard-Jones Ar–Ar interaction used in this study was rather well-known from work on the 3D gas. (Alterations in this energy due to the nearby presence of the solid adsorbent are present, but are often ignored.) The Ar–solid interaction was semi-empirical, with constants that are refined by comparison with the readily calculable Henry’s law constant. The gas– solid interactions for non-reactive gases on insulating solids are generally assumed to be sums over the atoms in the solid. (A different approach is needed for metallic substrates where the conduction electrons play a large role.) The evaluation of these sums requires a knowledge of the atom positions—straightforward for crystals if the exposed lattice plane is known, but not for amorphous or partially amorphous materials. In the case of graphitized carbon black, basal planes are exposed, for the most part, and the structure is well-known. If one takes the Ar–C atom potential to be a Lennard-Jones function with constants
8
W. Steele / Applied Surface Science 196 (2002) 3–12
Fig. 5. Isotherms for Ar on graphitized carbon black simulated at two reduced temperatures T ¼ kT=egg (egg is the Ar–Ar interaction well depthÞ ¼ 0:67 (squares) and 1.00. These are compared with an experimental isotherm measured at 78 K. y is the coverage in layers and p0 is the bulk vapor pressure.
(well-depth and size) that gives Henry’s law adsorption in agreement with high T, low coverage measurements, one obtains a simulated multilayer isotherm that agrees reasonably well with the low T data as indicated in Fig. 5. The rounded steps that are found are due to layer-by-layer adsorption at finite T on uniform graphite basal planes; the amount of rounding is most likely enhanced by edge effects associated with the finite extent of these planes in the experimental sample. If one accepts the site-wise summation models as adequate starting points for the moment, it is clear that both the location of these sites and the gas–site interactions are needed for an adsorption simulation to be realistic. One of the very large number of published adsorption simulations will be discussed here to illustrate how such computations can be used to gain a
deeper knowledge of this complicated process. The particular system was basically that of krypton in several straight-walled cylindrical pores [25]. The pore was constructed by first simulating a random close-packed cubic box of hard spheres [26]. Deleting the spheres from the desired cylindrical volume then produced a cylindrical cavity. Periodic boundary conditions were imposed on the ends of the cylinder and the system was transformed into an adsorbing medium by assuming Lennard-Jones interactions between each sphere and the atoms of the Kr adsorbate. Fig. 6 gives a view of the insides of two such cylinders. The first point to note is that the atom–atom interaction model will generate a heterogeneous adsorbing surface. Real adsorbents that approximate such a model are the MCM 21 zeolites that have been widely used for experimental adsorption studies over the past few years [27]. (Actually, these materials possess hexagonal pores rather than cylindrical, but this distinction is generally neglected. Also, their pore radii are roughly twice as large as those used in these simulations.) It is important to realize that the heterogeneity of such surfaces can be quite significant, depending upon the relative sizes of the adsorbate atom and the atoms of the silica adsorbents. To illustrate this point, adsorption energies defined as the minima in the gas–solid interaction potentials for Kr atoms at various positions on the walls were calculated for a net of about 3000 equally spaced points. Model walls with Lennard˚— Jones gas–solid atom sizes of 1.84, 2.30, 3.05 A with some deleted wall atoms—denoted by sc for ˚ without and with a larger small corrugation and 3.15 A
Fig. 6. Views down the axis of two cylindrical pores with walls ˚ made up of atoms. On the left, the atoms are small (1.84 A ˚ ); diameter) compared to the Kr atoms shown in the center (3.78 A ˚ ) and the wall has been on the right the atoms are larger (3.14 A made rougher by randomly deleting one-half of the wall atoms in the innermost two layers.
W. Steele / Applied Surface Science 196 (2002) 3–12
number of deleted wall atoms denoted by lc. The minima of the interaction energy and the positions of these minima have the expected dependence upon the Kr position. For instance, a Kr atom directly over a solid atom gives a smaller minimum at a location closer to the pore axis than one which is over the midpoint of a pair or triplet of wall atoms because such an adsorbed atom will sink into the center of the minicluster and ‘‘touch’’ two or three solid atoms rather than the single wall atom of an atop position. These
9
variations have been quantified by constructing histograms of the distributions of adsorption energies U and of the distances of these minima from the pore axis Rm. These histograms are shown in Fig. 7 for pores that have been modeled such that the average Rm is ˚ and the average Um is 8.8 kJ/mol for all 10 A systems. It is clear that these pore walls are significantly heterogeneous and that the heterogeneity increases as the wall atom/adsorbate size ratio increases or as the wall is roughened by random
Fig. 7. Histograms for the distributions of pore radius Rm (top) as defined by the positions of the adsorption energy minima U and for the values of these minima (bottom) are shown for all pores considered except that with the smallest wall atom size (and the least heterogeneity). The curves for wall atom size 3.15 are omitted here.
10
W. Steele / Applied Surface Science 196 (2002) 3–12
deletion of some of the adsorbing atoms in the wall. Furthermore, the physical smoothness as measured by R becomes significantly less as the wall atom size and/ or the corrugation produced by missing atoms increases. These model pores were used for simulations of the adsorption of Kr at 120 K. The isotherms and the gas–solid contribution to the adsorption ener-
gies at submonolayer coverages are shown in Fig. 8. Monolayer coverages for all systems except the large corrugation (lc) case are estimated to be 33 atoms/nm. In the lc pore, the large number of ‘‘holes’’ in the walls allows for significant adsorption but it is hard to define the ‘‘monolayer’’ for such a rough surface. It is clear in Fig. 8 (top) that a full (lc) pore contains more adsorbed
Fig. 8. Simulated isotherms and the gas–solid parts of the adsorption energy are shown on the top and the bottom, respectively. The isotherms were evaluated for adsorption up to pore-filling and both the adsorption and desorption branches of the curves are shown. The simulated energies are shown only for coverages up to the condensation points, which are slightly greater than monolayer.
W. Steele / Applied Surface Science 196 (2002) 3–12
atoms per unit length than those with smoother walls even when the average radius Rm is fixed. For submonolayer adsorption, the distributions of adsorption energies shown in Fig. 7 give rise to a nearly coverageindependent average gas–solid energy for the smallest wall atoms with the narrowest distributions of U. Much larger changes in this energy can be seen for larger wall atoms in pores both with and without imposed corrugation. (Note that the solid with size ˚ has adsorption properties almost parameter 3.15 A ˚ size parameter and identical to those with a 3.05 A a small additional wall corrugation.) The shapes of the isotherms at low coverages change significantly as the heterogeneity increases, going from curves with gradual almost linear rises for the surfaces lacking in strong sites to those with the rapid initial increases that are typically found for surfaces with strong gas–solid interactions. (Here, the use of the word site does not mean a preferred specific point for adsorption, but more generally, one of the 3000 points where the energies were calculated for use in the histograms.) Based on our understanding of the role of heterogeneity in physical adsorption [28], these changes in isotherm shape with changes in surface heterogeneity are as expected. However, it would be interesting to see if the adsorption energy distributions in the histograms of Fig. 7 could be used in the standard theories of adsorption on heterogeneous surfaces [27] to repro-
11
duce at least the monolayer portions of the isotherms shown here. Capillary condensation and hysteresis occur in all simulated isotherms shown here, with complete tube filling at 70–75 atoms/nm for all systems except the (lc) case. The positions and widths of the hysteresis loops vary somewhat with the changing surface heterogeneity. In spite of the presence of significant surface heterogeneity, the steps in these loops are essentially vertical on both the adsorption and desorption branches. Theory is unable to quantitatively explain this hysteresis at present, although it is likely that the considerable efforts being expended on a solution to this problem will give useful results soon. Finally, the dependence of the adsorbed phase densities upon distance from the pore axis is shown for coverages near monolayer completion in Fig. 9. Only for the two smoothest surfaces does one see anything like a single layer of atoms adsorbed near the wall. As the roughness increases, the atoms in the ‘‘monolayer’’ spread out over a range of distances, thus making the operational definition of monolayer completion more and more arbitrary and of less significance. We have devoted considerable space to discussion of the use of computer simulation not just to produce thermodynamic observables from a microscopic model but, in addition, to generate detailed data about
Fig. 9. The density of adsorbed atoms is shown as a function of distance from the pore axis for coverages close to monolayer completion. The formation of a layer at Ravg is expected for atoms adsorbed on a smooth wall.
12
W. Steele / Applied Surface Science 196 (2002) 3–12
the model surface properties that can enable critical tests of theory. It would seem that analyses which utilize the full power of the computer in this and other ways not discussed here will become an important aspect of our progress towards a full understanding of physical adsorption on or in complex solids.
[18] [19] [20] [21]
References [1] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, E. Teller, J. Chem. Phys. 21 (1953) 1087. [2] B.J. Alder, T.E. Wainwright, Phys. Rev. 127 (1962).; B.J. Alder, T.E. Wainwright, J. Chem. Phys. 33 (1960) 1439. [3] M.P. Allen, D.J. Tildesly, Computer Simulation of Liquids, Oxford Science Publications, 1987.; D. Frenkel, B. Smit, Understanding Molecular Simulations: From Algorithms to Applications, Academic Press, New York, 1996. [4] D. Nicholson, N.G. Parsonage, Computer Simulation and the Statistical Mechanics of Adsorption, Academic Press, New York, 1982. [5] W.A. Steele, The Interaction of Gases with Solid Surfaces, Pergamon Press, Oxford, 1974. [6] T.L. Hill, Statistical Mechanics, McGraw-Hill, New York, 1956. [7] T. Boublik, Mol. Phys. 29 (1975) 421. [8] D.G. Chae, F.H. Ree, T. Ree, J. Chem. Phys. 50 (1969) 1581. [9] M.J. Maeso, J.R. Solana, J. Chem. Phys. 102 (1995) 8562. [10] T. Boublik, Mol. Phys. 63 (1988) 685. [11] J. Talbot, D.J. Tildesley, J. Chem. Phys. 83 (1985) 6419. [12] C.W. Goulding, M. Rigby, Mol. Phys. 75 (1992) 623. [13] A. Bran˜ ka, K.W. Wojciechowski, Mol. Phys. 78 (1993) 1513. [14] I. Nezbeda, T. Boublik, O. Trnka, Czech. J. Phys. B 25 (1975) 119. [15] J. Viellard-Baron, J. Chem. Phys. 56 (1972) 4729. [16] J.A. Cuesta, D. Frenkel, Phys. Rev. A 42 (1990) 2126. [17] H. Schlacken, H.J. Mo¨ gel, P. Schiller, Mol. Phys. 93 (1998) 777.
[22] [23] [24]
[25]
[26]
[27]
[28]
D. Frenkel, R. Eppenga, Phys. Rev. A 31 (1985) 1776. P. Siders, Mol. Phys. 68 (1989) 1001. W.A. Steele, J. Chem. Phys. 65 (1976) 5256. A. Mulero, F. Cuadros, Langmuir 12 (1996) 3265; A. Mulero, F. Cuadros, Chem. Phys. 205 (1996) 379; A. Mulero, F. Cuadros, Chem. Phys. 156 (1991) 33; A. Mulero, F. Cuadros, Chem. Phys. 159 (1992) 89; A. Mulero, F. Cuadros, Chem. Phys. 160 (1992) 375; A. Mulero, F. Cuadros, Chem. Phys. 160 (1993) 375; A. Mulero, F. Cuadros, J. Phys. Chem. 99 (1995) 419; F. Cuadros, A. Mulero, L. Morala, V. Gomez-Serrano, Langmuir 17 (2001) 1576. J.D. Weeks, D. Chandler, H.C. Andersen, J. Chem. Phys. 54 (1971) 5237. J.J. Weis, Mol. Phys. 93 (1998) 361. W.A. Steele, in: J. Schwartz, C. Contescue (Eds.), Surfaces of Nanoparticles and Porous Materials, Marcel Dekker, New York, 1999. W.A. Steele, M.J. Bojan, Adv. Colloid Interf. Sci. 76–77 (1998) 153; B. McEnany, T.J. Mays, J. Rouque´ rol, F. Rodri uez-Reinoso, K.S.W. Sing, K.K. Unger (Eds.), Characterization of Porous Solids IV, Royal Society of Chemistry, 1997, p. 49.; M.D. LeVan (Ed.), Fundamentals of Adsorption, Kluwer Academic Publishers, Dordrecht, 1996. V.A. Bakaev, Surf. Sci. 215 (1988) 571; V.A. Bakaev, W.A. Steele, Langmuir 8 (1992) 148; A. Bonissent, B. Mutaftschieve, Phil. Mag. 35 (1977) 65. M. Grun, K.K. Unger, A. Matsumoto, K. Tsutsumi, Characterization of Porous Solids IV, Royal Society of Chemistry, 1997, p. 8.; P.J. Branton, P.G. Hall, K.S.W. Sing, Characterization of Porous Solids IV, Royal Society of Chemistry, 1997, p. 98.; P.I. Ravikovitch, G.L. Haller, A.V. Neimark, Adv. Colloid Interf. Sci. 76–77 (1998) 203. M. Jaroniec, R. Madey, Physical Adsorption on Heterogeneous Solids, Elsevier, Amsterdam, 1998.; W. Rudzin´ ski, D.H. Everett, Adsorption of Gases on Heterogeneous Surfaces, Academic Press, New York, 1992.