Adsorption of Spherical Molecules in Probing the Surface Topography

Langmuir 2002, 18, 3963-3979

3963

Adsorption of Spherical Molecules in Probing the Surface Topography: 2. Model of Conditional Probabilities Fre´de´ric Villie´ras,*,† Vadim Sh. Mamleev,‡ David Nicholson,§ and Jean-Maurice Cases† Laboratoire Environnement et Mine´ ralurgie, INPL and CNRS UMR 7569, BP 40, 54 501 Vandoeuvre les Nancy Cedex, France, Kazakh-American University, Satpaev Str., 18 a, Almaty 480013, Kazakhstan, and Department of Chemistry, Imperial College of Science, Technology and Medicine, Exhibition Road, SW7 2AY, London, United Kingdom Received November 5, 2001. In Final Form: December 18, 2001 Adsorption on a heterogeneous surface is considered within the framework of the lattice model taking into account a correlation in arrangement of nearest adsorption sites. The expression for an adsorption isotherm is derived in the mean field approximation. From an experimental isotherm, the proposed algorithm computes fractions of different adsorption sites, reduced adsorption energies, and a matrix expressing conditional probabilities of detection of different sites in a neighborhood of sites of a chosen kind. The elements of that matrix are supposed to be expressed via the fractions of the sites. In particular, in the case of the well-known model of patchwise heterogeneity, this matrix has diagonal form. In the latter case the algorithm allows one to estimate energies of lateral interactions for admolecules localized on different patches. The patchwise model is used for treatment of an isotherm of argon adsorption on muscovite. As established, the lateral interactions cannot be ignored even at very low degrees of surface coverage. Perhaps this fact is explained by grouping of the strongest adsorption sites into linear chains such as edges, growth steps, spiral dislocations, etc. However, the most reasonable explanation is accumulation of the adsorbate molecules around needle-shaped contacts between adsorbent particles. It follows from general considerations that adsorption sites in a neighborhood of such contacts must be the strongest ones and, hence, must be occupied first and foremost, i.e., at lowest pressures. In a vicinity of such contacts the adsorbate molecules can form rings that can be considered, in turn, as a one-dimensional lattice with coordination number 2. The next step of the theory should consist in determination of energies of lateral interactions by using standard systems with known crystallography.

1. Introduction The thermodynamic properties of adsorbed films are unambiguously determined by the properties of both surfaces and adsorbate molecules. Therefore it is legitimate to approach the problem of analyzing van der Waals relief of surfaces with thermodynamic information concerning adsorption. For many years the solution of a such problem seemed extremely promising for the study of surface physics;1-5 therefore this problem often attracted significant creative forces. At the same time we are not aware of any work in which information obtained in such a manner has been strictly interpreted at a molecular level. One of the main questions related to such an analysis is a choice of a priori information necessary for stable computations of the statistical distribution, F(), of adsorption sites (centers) with respect to their binding energy, , with an adsorbate. This choice is complicated * To whom correspondence should be addressed. E-mail: [email protected]. † Laboratoire Environnement et Mine ´ ralurgie, INPL and CNRS UMR 7569. ‡ Kazakh-American University. § Imperial College of Science, Technology and Medicine. (1) Roginskii, S. Z.. Adsorption and Catalysis on Heterogeneous Surfaces; Academy of Sciences of the USSR: Moscow, 1948 (in Russian). (2) Jaroniec, M.; Madey, R. Physical adsorption on Heterogeneous Solids; Elsevier: Amsterdam, 1988. (3) Rudzin˜ski, W.; Everett, D. H. Adsorption of Gases on Heterogeneous Surfaces; Academic Press: London, 1991. (4) Jaroniec, M.; Bra¨uer, P. Surf. Sci. Rep. 1986, 6, 65. (5) House, W. A. Colloid Sci. 1983, 4, 1.

due to lack of a strict theory of formation of surfaces. At present the Monte Carlo simulation6,7 can hardly give an unambiguous answer on how a real surface is structured, even if the corresponding adsorbent is crystalline. Somehow this answer is hidden in experimental thermodynamic information. If any short-range ordering (SRO) on a surface does not exist, it is natural to assume that the function F() is a smooth one and it can be evaluated8-11 by the known algorithms of regularization.12 If SRO exists, one may subdivide adsorption sites into groups with similar properties and search the discrete sets {Fi}, {i}. At present any of these approaches can give rise to many doubts. However, earlier13 we showed that for surfaces of crystals, like muscovite, both of them give similar results. Hence, in principle, SRO can be established by handling of an experiment. Final proofs of a discrete nature of energy distributions for monocrystals were not obtained as yet. However, formation of surfaces of crystals in the form of a manifold of cleavage planes is accepted in classic theory.14,15 Apparently singular crystal faces make up energy-optimal (6) Bakaev, V. A. Surf. Sci. 1988, 198, 571. (7) Bakaev, V. A. Langmuir 1992, 8, 1372, 1379. (8) Merz, P. J. Comput. Phys. 1980, 38, 64. (9) Von Szombathely, M.; Bra¨uer, P.; Jaroniec, M. J. Comput. Chem. 1992, 13, 17. (10) Koopal, L. K. Vos, C. H. W. Langmuir 1993, 9, 2593. (11) Jagiello, J. Langmuir 1994, 10, 2778. (12) Tikhonov, A. N.; Arsenin, V. Ya. The Methods of Solution of Incorrect Problems; Nauka: Moscow, 1979 (in Russian). (13) Mamleev, V. Sh.; Villie´ras, F.; Cases, J. M. submitted to Langmuir

10.1021/la0116387 CCC: $22.00 © 2002 American Chemical Society Published on Web 04/19/2002

3964

Langmuir, Vol. 18, No. 10, 2002

surface.15 Predominant existence of basal planes on surfaces of carbon black and some halides was proven in the adsorption experiments.16 In this paper we shall accept discreteness of an energy distribution a priori, although it does not exclude additional difficulties. Among all variants of equations describing adsorption on different sites, only the simplest Langmuir local isotherm, which ignores lateral interactions, minimally touches on the question about a surface structure. It is noteworthy that because of lack of ponderable arguments to the advantage of stricter models, the Langmuir isotherm was used in publications most often.1-5 At the same time as coverage approaches a monolayer, lateral interactions of admolecules clearly influence thermodynamics of adsorption; therefore a fine resolution of F(), particularly in the region of low adsorption energy, seems to be impossible without a correction for lateral interactions. Perhaps, the main reason of failures in analysis of surface heterogeneity is connected with this circumstance. When solving the problem of F() evaluation, some corrections for lateral effects were made on the basis of the models of the pathchwise and random topography.4,5 However, as a rule, even individual cleavage planes contain different kinds of adsorption sites. Thus, these models are rather rough. In the past decade the models with intermediate topographies were used in solving the direct problem, i.e., for a generation of overall isotherms with artificially assigned topography.17-20 For example, Zgrablich and Riccardo et al.18 investigated the influence of the change of topographies within the framework of the, so-called, generalized Gaussian model.21 Different topographies are described in terms of the pair correlation function in arrangement of sites with different energies. This model seems semiempirical. Indeed, it very hard to imagine how this model could be used for description of adsorption on crystalline surfaces, where mutual dispositions of adsorption sites are strictly correlated but the sites have different adsorption energies. Generally, it seems strange that despite weakness of the quantitative theory of physical adsorption, different authors digress from the simplest model of a surface consisting of a set of cleavage planes with different fractions (“concentrations”). Against other approaches that we know, the ideas of Tovbin22,23 seem most promising for investigation of adsorption on crystals. He selected the variant of a crystalline surface as paramount for a practical experimental examination. In this paper we shall adhere to his approach. In the second section we shall derive the equation for an integral isotherm of localized adsorption with known correlation in arrangement of nearest adsorption sites. The theory of localized adsorption on heterogeneous surfaces can be developed by Bogolubov’s method24 in the (14) Laudise, R. A. The Growth of Single Crystals; Prentice Hall: Englewood Cliffs, NJ, 1970. (15) Dunning, W. J. In The Solid-Gas Interface,1; Flood, E. A., Ed.; Marcel Dekker: New York, 1967. (16) Thomy, A.; Duval, X.; Regnier, J. Surf. Sci. Rep. 1981, 1, 1. (17) Steele, W. A. Langmuir 1999, 15, 6083. (18) Riccardo, J. A..; Pereyra, V. D.; Rezzano, J. L.; Rodrigues, D., Zrgablich, G. Surf. Sci. 1988, 204, 286. (19) Riccardo, J. A.; Chade, M. A.; Pereyra, V. D.; Zgrablich, G. Langmuir 1992, 8, 1518. (20) Zgrablich, G.; Zuppa, C.; Ciacera, M.; Riccardo, J. A.; Steele, W. A. Surf. Sci. 1996, 356, 257. (21) Ripa, P.; Zgrablich, G. J. Phys. Chem. 1975, 79, 2118. (22) Tovbin, Yu. K. Teor. Eksp. Khim. 1982, 18, 417 (in Russian). (23) Tovbin, Yu. K. Theory of Physical Chemistry Processes at GasSolid Interface; Nauka: Moscow, 1990 (in Russian). (24) Tyablikov, S. V. Methods of the Quantum Theory of Magnetism; Nauka: Moscow, 1982 (in Russ.).

Villie´ ras et al.

mean field22 approximation. A more sophisticated approach23 allows one to describe localized adsorption in the quasichemical approximation. Mathematically the accuracy of the quasichemical approximation is higher. However, because of some uncertainties in physical information concerning real systems, an advantage of this approximation in comparison to the mean field approximation becomes questionable. In this paper we use the mean field (molecular-field) approximation as a compromise between accuracy and simplicity, perhaps, optimally at least at the current stage of development of the theory. We repeat the results of Tovbin22 again by the standard method25,26 in order to avoid complicated terminology accepted in quantum statistics22 that, perhaps, hinders in understanding of this approach. In addition to these results we shall obtain an expression for the grand potential predicting points of two-dimensional (2D) phase transitions of the first kind in adsorbed films. Possibility of multiple 2D transitions for the patchwise model seems trivial. However, Hill25 predicted the theoretical possibility of multiple 2D transitions in the case of random topography as well, which later on was confirmed by the Monte Carlo simulation.27 Steele17 mentioned 2D phase transitions for noble gases when investigating the heterogeneous surface of carbon black. Zgrablich et al. detected clusterization20 of adsorbate while simulating surfaces with intermediate topographies by the Monte Carlo method; however they did not study the conditions of phase transitions. Heterogeneity leads to lowering of critical temperature of 2D condensation. That is why for adsorption of noble gases on heterogeneous surfaces 2D transitions are not observed at the temperature of liquid nitrogen. For adsorption of Ar on crystals of TiO2, the 2D condensation was not detected28 even up to the temperatures 15-20 K. Nevertheless, investigation of the critical temperature of 2D transitions is of theoretical interest, since it could testify in favor of a choice of a correct model of heterogeneity. However, a detailed examination of this question remains beyond the present article, since for the system under consideration13 (Ar-muscovite at 78 K) adsorption occurs in the supercritical region. In the third section we briefly consider integral equations that can be formulated according to the cluster models of adsorption. We try to confront the model of Riccardo et al.19 with the model used. The difference is that in the one case the correlation function is used as a weighting coefficient, whereas in another case it enters into the kernel of the integral equation. The fourth section is devoted to the method of evaluation of discrete energy distributions by means of a new algorithm “validating admolecules’ fractional energies” (VADFRE). Development of this algorithm is the main purpose of the presented work. The reader can regard the above-mentioned sections as preparatory information. In the fifth and sixth sections, the results of VADFRE are shown in handling artificial isotherms and the isotherm of argon adsorption on muscovite. 2. Lattice Model of Adsorption on a Heterogeneous Surface in the Mean Field Approximation Imagine a surface carrying I kinds of distinguishable adsorption sites. The number of each kind of sites Bi (25) Hill, T. L. J. Chem. Phys. 1949, 17, 762. (26) Nicholson, D. J. Chem. Soc., Faraday Trans. 1 1975, 71, 238. (27) Gordon, R. J. Chem. Phys. 1968, 48, 1408. (28) Morrison, J. A.; Los, J. M.; Drain, L. E. Trans. Faraday Soc. 1951, 47, 1023.

Surface Topography

Langmuir, Vol. 18, No. 10, 2002 3965

function Ξ is represented by the series29 B

∑ exp(Nµ/kT)QN

Ξ)

(2)

N)0

where QN is the canonical partition function, k is Boltzmann’s constant, T is a temperature, and µ is the chemical potential of the gaseous adsorbate in a bulk phase. The function QN, in turn, is expressed as

QN ) Figure 1. The model of a regular surface with different adjacent adsorption sites.

(1 e i e I) is great enough, i.e., Bi . 1. Ni molecules of an adsorbate occupy sites of the ith kind. The numbers for Ni are also very large Ni . 1. Our problem is to find the dependencies of the partial occupancies θi(P) ) Ni(P)/Bi upon the adsorbate pressure P in a bulk phase. Thereupon, one can readily find an overall degree of a monolayer coverage that is expressed as the sum I

Θ(P) )

I

Ω ) -kT ln Ξ At the thermodynamic limit B f ∞, the partition function can be restricted by its maximum term

Ξ ) exp(Nµ/kT)QN

where B ) ∑i)1IBi is the total amount of adsorption sites (centers) on the surface and Fi is the fraction of sites of the ith kind in the total adsorption capacity of the monolayer. Suppose that an adsorbate molecule being localized in any site of the ith kind on the surface interacts only with molecules localized in the nearest sites of the jth kind (see Figure 1). We note that in the mean field approximation the models, which take into account long-distance interactions and only nearest interactions, in a phenomenological aspect, are equivalent,22 and so consideration of only the nearest interactions does not violate generality. Since the adsorbed molecules interact, in addition to the manifold {Fi} one should know also the manifold {fij} in which the element fij means a conditional probability to detect a site of the jth kind near a site of the ith kind. The numbers fij are connected with the numbers of pairs of adjacent sites Bij. On the strength of the following balance

(1 + ∆ij)Bij ∑ j)1

one may write

fij ) Bij(1 + ∆ij)/ZiBi

(1)

I

fij ) 1 ∑ j)1 where ∆ij is the Kronecker delta: ∆ij ) 1 (when i ) j); ∆ij ) 0 (when i * j). Note that the numbers Bij are necessary only to illustrate the derivation of an expression for the adsorption isotherm, whereas in particular calculations only the matrix {fij} is needed. The elements fij (1 e i e I, 1 e j e I) are determined from real crystallography of the surface (see Figure 1). Let us derive the equation for the adsorption isotherm from the Gibbs grand ensemble in which the partition

(4)

where N corresponds to the maximum term in series (2). In the grand ensemble the numbers Ni (1 e i e I) are independent, i.e., ∂N/∂Ni ) 1. Thus, determination of the maximum term in Ξ leads to the follow system of equations:

∂ ln Ξ/∂N ) ∂ ln Ξ/∂Ni ) ∂ ln QN/∂N + µ/kT ) ∂ ln QN/∂Ni + µ/kT ) 0 (5) 1eieI According to the lattice model, adsorbed molecules can be disposed only in the neighborhood of adsorption sites (or close to centers of potential “dimples”). A molecule sunk into such a “dimple” makes small oscillations. After summing eq 3 over the vibrational levels, one can rewrite QN in the form25,26 I

QN )

I

(3)

where gn is the degeneracy factor and En is a value of the full internal energy of the adsorption system at the nth energy level. The summation in eq 3 is carried out over all energy levels of the system. The natural variables in the grand ensemble are µ, T, and B, and the characteristic function is the grand potential that is expressed as

θi(P)Bi/B ) ∑ Fiθi(P) ∑ i)1 i)1

ZiBi )

∑n gn exp(-En/kT)

q˜ iN ∑n (∏ i)1

(n)

i

)gn* exp(-Φn/kT)

(6)

where q˜ i is the partition function for the vibrational energy of admolecules localized in the adsorption sites of the ith kind, gn* is a new degeneracy factor equal to the number of ways in which N molecules can be arranged on the surface under a constant potential energy Φn of their interactions with the surface and with each other. The sum q˜ i is expressed30 as 3

[1 - exp(-hνip/kT)]-1 ∏ p)1

q˜ i ) exp(-v,i/kT)

where h is Planck’s constant; νip (p ) 1, 2, 3) are the vibration frequencies along the normal coordinates x, y, z; v,i ) (1/2)h(νi1 + νi2 + νi3) is a zero point energy. (29) Hill, T. L. Statistical Mechanics, Principles and Selected Applications; McGraw-Hill: New York, 1956. (30) Lopatkin, A. A. Theoretical Fundamentals of Physical Adsorption; Moscow State University: Moscow, 1983 (in Russian).

3966

Langmuir, Vol. 18, No. 10, 2002

Villie´ ras et al.

In the existing models, as a rule, the polarization of admolecules is supposed to be independent of their molecular surroundings. In such a case the overall potential of lateral interactions is described as a sum of pair interactions. Thus, one can write the total potential energy of interactions of the admolecules with the surface and with each other as I

Φ)-

Nii - ∑ Nijwij ∑ i)1 N

we ignored correlation caused by the interaction between admolecules. In such approximation the degeneracy factor g*({Ni},{Nij}) is actually replaced by g*({Ni}). Note that, according to eq 9, the energy does not vary at permutations of Ni molecules within the manifold of Bi sites. For each group of Ni molecules the Bi!/Ni!(Bi - Ni)! variants of such permutations are possible; hence, their total number for the whole surface gives

(7)

I

g*({Ni}) )

ij

where Nij and wij are the number of pairs and the pair interaction energy of molecules localized in the nearest sites of the ith and jth kinds. Since the potential energy of attraction is negative, we have put the negative signs in the sums in eq 7. This was done solely for the sake of convenience, since manipulation with the absolute values of the energies i > 0, wij > 0 is more straightforward. The interaction energy i with the adsorbent is taken from the zero point vibration energy i ) i* - v,i, where i* is the absolute value of the minimum of surface potential. This allows magnitudes v,i to be eliminated from the subsequent equations. With a constant N the quantization of the potential energy Φ in the system is connected with a possible redistribution of admolecules between the adsorption sites. The energy Φ varies in discrete manner and forms the set of levels {Φn}. However, summation in eq 6 can be omitted, since we shall have to use only one term in QN (one value of Φ) corresponding to the maximum term in Ξ. Taking account of eq 7, one may write the partition function QN in the form

Bi!/Ni!(Bi - Ni)! ∏ i)1

Now, the partition function (8) is expressed as I

QN )

[Bi!/Ni!(Bi - Ni)!]{qi exp{[i + ∏ i)1 I

(1/2)

wijNjBij(1 + ∆ij)/(BiBj)]/kT}}N ∑ j)1

i

By using the Stirling formula ln x! ≈ x ln x - x from the system (5), we have

∂ ln QN/∂Ni + µ/kT ) ln[(1 - θi)/θi] + ln ci + I

wijθjBij(1 + ∆ij)/(BikT) + µ/kT ) 0 ∑ j)1

Nijwij/kT)( ∏[qi × ∑ N i)1

QN ) g*({Ni},{Nij}) exp(

where ci ) qi exp(i/kT). The chemical potential29-31 of an ideal monatomic gas is expressed by the well-known formula

Taking into account eq 1, one can finally write the overall isotherm as

ij

exp(i/kT)]Ni) (8) where qi ) exp(v,i/kT)q˜ i. The degeneracy factor g* is equal to the number of ways in which a possible set {Ni} can generate the manifold of interacting pairs {Nij} while conserving the constancy of N and Φ. Let some molecule be localized in a site of the ith kind. The potential of its interaction with a nearest-neighbor site of the jth kind can be zero, if that site is empty, or equal to wij, if the site is occupied by another molecule. According to the molecular-field approximation the probability of detecting a molecule in the jth site is proportional, on average, to the occupancy of the jth sites. Therefore one may write the average potential of interactions of the chosen molecule with the molecule situated in the jth site as wijNj/Bj. The total potential of interactions of all molecules localized in the sites of the ith kind with their nearest surrounding is equal to ∑j)1IwijNjBij/Bj. The probability of occupancy of the ith site, in turn, is equal to Ni/Bi; hence, the energy Φ in such an approximation can be expressed as I

I

(11)

µ ) kT ln[h3/(2πm)3/2(kT)5/2] + kT ln P

I

Φ)-

(10)

I

Nii - (1/2) ∑ Ni/Bi ∑ wijNjBij(1 + ∆ij)/Bj ∑ i)1 i)1 j)1

(9)

The factor 1/2 in eq 9 is needed to avoid the double summation of the same energy contributions. The Kronecker symbol appears for the correct summation of the contributions of the same sites, since the Bii pairs contribute the magnitude of Biiwii(Ni/Bi)2 to the energy Φ. While determining the probability of the sites’ occupancy,

I

Θ(P) )

Fiθi(P,i) ∑ i)1

(12) I

θi(P,i) ) P/{P + A0,i exp[-(i + Zi

wijfijθj)/kT]} ∑ j)1

(13)

1eieI where A0,i ) (kT/qi)(2πmkT/h2)3/2 is the Langmurian constant for the sites of the ith kind. Equations 12 and 13 predict (just like, e.g., the van der Waals equation) that with large values of the elements {wij} (or at low temperatures) the function Θ(P) can exhibit, so-called, loops29 (see Figure 2). In the region of the loop, the function P(Θ) ) P has three roots and one root outside the loop. In conformity with each of the roots, the system of eqs 13 gives three manifold {θi}. The appearance of the loops means an occurrence of phase transitions of the first order. However, real isotherms in the subcritical region should have a shape of steps and not loops. Condition (5) defines the most probable thermodynamic state of the system in the grand ensemble. For loop-shaped isotherms this condition leads to three solutions. From these solutions one should find the only one corresponding to the global maximum of Ξ. The maximal value of Ξ (or, equivalently, the minimal value of Ω) calculated via eq (31) Ono, S.; Kondo, S. Molecular Theory of Surface Tension in Liquids; Springer-Verlag: Berlin, 1960.

Surface Topography

Langmuir, Vol. 18, No. 10, 2002 3967

chemical nature and, consequently, close adsorption energies, we can rewrite eq 15 in the form I

Θ(P) )

Figure 2. The treatment of the loop-shaped isotherm that exhibits two-phase transitions of the first kind. Inside the loops the isotherm is plotted point by point by estimation of a maximum thermodynamic probability via the grand potential.

4 serves as a measure of the maximal thermodynamic probability. By using eqs 10 and 11, one can find the following formula:

Ω/(kTB) ) -(ln Ξ)/B ) I

(1/2)

∑ i)1

I

θiFiZi

∑ j)1

I

wijfijθj/kT -

Fi ln[1/(1 - θi)] ∑ i)1

(14)

For plotting an isotherm in the subcritical region (see Figure 2), one should find three sets {θi} for each of the roots. Thereupon, expression (14) should be calculated for each set. One of three roots can be immediately rejected, since it corresponds to the minimum of the thermodynamic probability. Values of Ω/(kTB) for the two other roots must be compared. The minimal value of Ω/(kTB) corresponds to the sought for point on the isotherm. The isotherm being calculated in such a manner becomes stepwise and not loop-shaped. In general, eqs 12, 13, and 14 allow investigation of critical conditions for 2D transitions only numerically. However, one can conclude in advance that for systems with two kinds of sites (I ) 2) the double 2D transitions (see Figure 2) can take place independently of topography of a surface, if adsorption energies for different sites considerably differ. Indeed, if sites are occupied separately, the model with general topography turns into the patchwise model with effective parameters. 3. Reduction to the Integral Equation If lateral interactions are negligible, all adsorption sites on a surface are independent and, according to the model of localized adsorption, an overall isotherm can be expressed in the form B

Θ(P) ) B-1

θi(P,i) ∑ i)1

(15)

where θi denotes a probability of occupancy of a site with energy i. In contrast to eq 12, summation in eq 15 is carried out over all sites of the lattice (I ) B). Within the framework of the hypothesis of localized adsorption in the case of independence of sites, eq 15 is strict. However, the number of adsorption sites B is very large (say, B is equal to 6.02 × 1023 per kilogram of an adsorbent13). Thus, the analysis of heterogeneity by eq 15 is impossible. The next step of the theory can consist in a fragmentation of sites by groups. Indeed, if Bi sites have similar

∑ i)1

I

(Bi/B)θi(P,i) )

Fiθi(P,i) ∑ i)1

(16)

where θi now is the average occupancy of the ith kind of sites and Fi is the factor of degeneracy of energies of adsorption sites. For surfaces with crystalline structure, I can be sufficiently small, say, I = 3-10. Thereby, eq 16 can be used for the heterogeneity analysis. For amorphous surfaces eq 16 should be rewritten in the integral form

Θ(P) )

∫

max

min

F() θ(P,) d

(17)

where F() is the energy distribution function. The function F() can be evaluated by inverting eq 17. The known methods of regularization8-13 are based on smoothing of F(). However, any physical proofs of F() smoothness are absent, therefore, in contrast to the sets {Fi}, {i}, physical interpretation of F() is problematic. Equations 16 and 17 can be called sitewise models of heterogeneity.17 If lateral interactions are considerable, one should use cluster models23 taking into account interaction of some admolecule in the ith site with its neighbors. The model considered in the previous section is the simplest variant of the cluster approach. In conformity with eq 12 the cluster model can be written as follows: I

Θ(P) )

Fiθi(P,i,A0,i,Zi,{wij},{fij}) ∑ i)1

(18)

Equations 12 and 13 obtained in the previous section strictly correspond to series (18), and in this sense they are correct. Note that we omitted the variables A0,i, Zi, {wij}, and {fij} in the left-hand side of eq 13 only for simplicity, whereas, in fact: θi ) θi(P,i,A0,i,Zi,{wij},{fij}). According to the physical considerations, some sites belong to the ith manifold only if apart from the equal adsorption energies they have identical first coordination spheres consisting of other sites. This implies that the matrix {fij} in each line must contain only Zi nonzero elements. The analysis of heterogeneity by eq 18 is possible for crystalline surfaces23 with known structure. For example, if an adsorbent presents monocrystals with a surface consisting of different cleavage planes, one can preset {fij} from crystallographic data (see Figure 1). The values {wij} can be evaluated from heats of adsorption or from another thermodynamic information. Thus, the series (18) needs to be inverted only with respect to {Fi} and {i}. Such analysis means determination of fractions of different cleavage planes and energies of sites located on them. Two obstacles arise when using eq 18. First, if cleavage planes have microscopic sizes, one should take into account the boundary effects at edges of the planes. To overcome this difficulty Steele17 has offered to use the notion of a supersite. Actually it is the wellknown patchwise model but with correction for the boundary effects. In our opinion, Tovbin23 has developed the more promising concept. He has considered sites belonging to edges (e.g., between terraces and steps of growth17) as defined kinds of adsorption sites. If inversion of eq 18 would allow a detection of fractions not only of cleavage planes but also of their boundaries, one could estimate sizes of cleavage planes. Unfortunately, it seems

3968

Langmuir, Vol. 18, No. 10, 2002

Villie´ ras et al.

that accuracy of the existing theories is not sufficient for such computations. Second, the elements {fij} are connected with the fractions {Fi}. When solving the inverse problem, it is necessary to use the relationships

θi ∑ i)1

Z+1

∑

E i ) i + w

Z+1

1

θj + w1λ0θ + w2λΘ

(25) (26)

j*i∈AZ+1

fij ) fij({Fi})

(1 e i e I, 1 e j e I)

(19)

For cleavage planes with known crystallography (see Figure 1) an explicit definition of eqs 19 does not give rise to special difficulties but, before solving the inverse problem, a sequence of adsorption sites belonging to different planes in ascending order of adsorption energies should be predetermined. In other words, one cannot avoid preliminary estimations of adsorption energies, since otherwise after solving the inverse problem it will be impossible to attribute sites with the energy i to certain crystallographic plane. Both kinds of difficulties can be overcome for simple crystalline surfaces. Adaptation of the lattice model for amorphous surfaces and for surfaces with fuzzy crystallography is much more difficult. In this case we are forced to assume that sites with energy i are surrounded with different coordination spheres, and adsorption on sites of the ith kind implies the following averaging L

θi(P,i) ) L-1

θ)

∑θic(P,i,Ai,0c,Zic,{wij}c,{fij}c) c)1

(20)

where the summation using the index “c” is carried out over L possible coordination spheres of the ith kind of sites. Models of a surface allowing averaging (20) are absent. Therefore further approximations are necessary. For example, one can assume that the numbers Zic, wijc, and A0,ic are constant all over the lattice, i.e., Zic ) Z, wijc ) w, and A0,ic ) A0. Then the summation in eq 20 changes to

∑

λ0θ )

θ

(27)

j∉AZ+1,rijer0

λΘ )

∑

Θ

(28)

rij>r0

where r0 is some radius of correlation. The authors considered a cluster AZ+1 consisting of a central site surrounded with the first coordination sphere of other sites. The integers λ0 and λ denote the numbers of sites remaining beyond the cluster AZ+1; w1 and w2 are averaged energies of long-distance interactions. The averaging implied in eqs 27 and 28 can give rise to questions. However, the main difficulty is associated with eq 23. Accuracy of the calculation of lateral effects subsides from the center of the cluster to its periphery. That is why in this modeling the accuracy can be considerably increased, if the central site is considered separately. Without loss of a generality, we can rewrite eq 23 in the form

Θ(P) )

∫

max

min

∫

max

‚‚‚

min

F()fZ(,{˜ 1,‚‚‚,˜ Z})θc(P,,{˜ 1,‚‚‚,˜ Z}) d d˜ 1 ‚‚‚ d˜ Z (29)

We conserved the superscript “c” to emphasize that θc means occupancy of central sites in clusters and not the averaged occupancy. Let us retain in eq 26 only the terms expressing nearest interactions, i.e. Z

L

θi(P,i) )

fZθic(P,i,{1c,‚‚‚,Zc}) ∑ c)1

where fZ(1c,‚‚‚,Zc) is the probability to find a group of sites (cluster), where the ith (central) site with the energy i is surrounded with sites having the energies {1c,‚‚‚,Zc}. If L f ∞, eq 21 can be represented in the integral form

θi(P,i) )

∫

max min

∫

max

‚‚‚

E h )  + wZ

fZ(i,{˜ 1,‚‚‚,˜ Z})θic(P,i,{˜ 1,‚‚‚,˜ Z}) d˜ 1 ‚‚‚ d˜ Z (22)

Θ(P) ) max min

∫

max

∫

max

min

f(,ξ)θ(P,ξ) dξ

By transposing the averaged θc(P,) from under the multiple integral in eq 29, we have

Θ(P) )

∫

max

min

F()θ(P,) d

where

θ(P,) ) P/{P + A0 exp[-( +

‚‚‚

min

fZ+1(1,‚‚‚,Z+1)θ(1,‚‚‚,Z+1,P) d1 ‚‚‚ dZ+1 (23)

θi ) P/[P + A0 exp(-Ei/kT)]

θj ∑ j)1

Equation 29 is essentially simplified, if using an averaged value of Ei for each central site. Being averaged over the lattice, the average number of sites with the energy from ξ to ξ + dξ in the first coordination sphere of sites with energy from  to  + d is equal to Zf(,ξ) dξ, where f(,ξ) is the pair correlation function. Thus, the averaged value of Ei for sites with the energy from  to  + d is equal to

min

where fZ is the multivariate correlation function. We used the tildes in eq 22 to denote the sites at the coordination sphere. Riccardo et al.19 used the notion of multivariate correlation function to describe an overall isotherm. If our notation is used in some places, it has the following form:

∫

E i ) i + w

(21)

(i ) 1, ‚‚‚, Z + 1) (24)

wZ

∫

max

min

f(,ξ)θ(P,ξ) dξ)/kT]} (30)

Equations 23 and 25 are correct but they eliminate a specific character of different adsorption sites. Moreover, it is hard to imagine how the multivariate correlation

Surface Topography

Langmuir, Vol. 18, No. 10, 2002 3969

function could be measured or calculated from an experiment. Even the pair correlation function apparently cannot be measured. However, the latter may be preset on the basis of the existing theories of a surface. For example, if a surface consists of macroscopic patches and boundary effects are negligible, only similar adsorption sites are located in each coordination sphere, i.e., f(,ξ) ) δ( - ξ). Then eq 30 comes to the isotherm of the patchwise model13

θ(P,) ) P/{P + A0 exp[-( +wZθ)/kT]}

(31)

If adsorption sites are randomly “mixed” over a surface, we have f(,ξ) ) F(ξ). In this case eq 30 turns into the isotherm of the random model25

θ(P,) ) P/{P + A0 exp[-( + wZΘ)/kT]}

(32)

Figure 3. Typical adsorption isotherm on the logarithmic scale of pressure. The fractional contributions ni(P) (1 e i e I) of different sites (see Figure 4) are summed up in such a way that the overall isotherm becomes the smooth function without inflection points.

These examples show that for characterizing a surface we need two functions F() and f(,ξ). That is why it is not sufficiently clear which class of surfaces was spanned in the papers of Zgrablich, Riccardo et al.,19,20 where only one function, namely f(,ξ), was used for the simulation of heterogeneity. In a general case the mean field approximation gives23 the following equations:

Θ(P) )

∫

max

min

F()θ(P,) d

(33)

θ(P,) ) P/{P + A0() exp[-( + Z()

∫

max min

w(,ξ)f(,ξ)θ(P,ξ) dξ)/kT]} (34)

The system of eqs 33 and 34 describe adsorption on a crystalline surface, if I

F() )

Fiδ( - i) ∑ i)1

f(,ξ) ) fijδ( + ξ - i - ξj) where δ(x) is the Dirac delta function. Then eqs 33 and 34 turn to eqs 12 and 13. It is easy to show that patchwise model (31) and random model (32) for discrete distributions are obtained, if, in addition to the conditions Zic ) Z, wijc ) w, and A0,ic ) A0, the matrix {fij} has the form:

To calculate the overall isotherm (12), it is necessary to solve the system of transcendental equations (13). In the next section we shall discuss not only the method of the isotherm computation but also a more elaborate procedure of evaluating the {Fi}, {i}, and {fij} manifolds from a given overall isotherm by using relationships (19). 4. Validating Fractional Energies of Admolecules The volumetric technique allows one to measure an amount of a substance n(P) adsorbed per gram of an adsorbent as a function of the equilibrium pressure P (or logarithm of pressure) in a bulk phase (Figure 3). The n value is usually expressed either in the form of an amount of moles or in the form of the corresponding volume of a gaseous adsorbate per gram of an adsorbent at the standard temperature and pressure (STP).

Figure 4. The stages of occupancy of a monolayer with increasing pressure in a bulk phase.

The experiments are usually carried out in the subcritical region T < T0, where T0 is the critical temperature of an adsorbate condensation into a liquid. As the pressure P increases to the pressure P0 of a saturated vapor, the contribution from multilayer adsorption rises. Unfortunately, the analytical theory of multilayer adsorption on heterogeneous adsorbents is still at a stage of development. One can expect, however, that for a number of surfaces, in particular, for those of oxides,7 the occupancy of a monolayer (the first layer of adsorption) is completed, when the amount of adsorbate in the second layer and subsequent ones is still negligible. Thereby, the asymptote of an isotherm at P f P0 is determined solely by accumulation of an adsorbate in a multilayer film. An elementary correction for the effect of multilayer adsorption is possible by the Brunauer-Emmet-Teller (BET) equation: n(P) ) nB(P)(1 - P/P0), where nB(P) is an observed (experimental) isotherm of multilayer adsorption. Despite “roughness” of a surface, an arrangement of admolecules in a monolayer hardly differs greatly from a close-packed one. However, one should keep in mind that the number of adsorption sites B can differ considerably from the adsorption capacity na (see Figure 4), inasmuch as B is defined by a surface topography, whereas the value of na is defined by a density of two-dimensional lattice of an adsorbate. Only in the case of ideally complementary lattices may one expect B ) na. Monolayer adsorption occurs by way of sequential coverage of sites having a different affinity with admolecules (parts a, b, and c of Figure 4). Occupancy of the strongest sites obeys models of localized adsorption. On weak sites, mobile adsorption (Figure 4c) can take place. From the standpoint of the minimum internal energy of an adsorbed film itself, at a complete coverage of a monolayer the formation of a two-dimensional lattice of admolecules (Figure 4d) is preferable. If the tendency for closest packing prevails, admolecules must leave some

3970

Langmuir, Vol. 18, No. 10, 2002

Villie´ ras et al.

strong sites empty. Indeed, such a phenomenon was detected for adsorption of noble gases on carbon black.16 However, another phenomenon, where the adsorbate lattice is defective and the strong sites remain occupied (Figure 4c), in our opinion, takes place more frequently because in most experiments with heterogeneous surfaces steplike isotherms of multilayer adsorption are not observed.32 This fact confirms defectiveness of a monolayer due to “roughness” of a surface.16 In any case the effects at the stage of residual monolayer coverage do not hinder detection of strong sites, as these are manifested within the initial region of an isotherm. Weak sites (sites of the third kind in Figure 4c) can be formally considered as effective sites of one type.13 Some semiempirical local isotherm for these sites can describe averaged (apparent) properties of a surface at the stage of termination of monolayer formation. If original information includes only an initial fragment of an isotherm, one may estimate only the number of strong sites on the surface. Because of lack of data for occupancy of weak sites, it is impossible to carry out an extrapolation (see Figure 3) of an isotherm toward P f ∞. Attempts at an evaluation of the adsorption capacity na from isotherms restricted to the region of a monolayer adsorption were considered in a number of publications.10,33,34 In our opinion,35 this would only be possible in some special cases, e.g., if a degenerated energy distribution were obtained in the calculations.13 In this case an increase in the number of expansion terms for an overall isotherm does not influence the accuracy of approximation of the isotherm.13 However, in general, the extrapolation of an isotherm is illegitimate if the calculations do not include the normalization condition

Suppose an isotherm in the experiment be measured at the M points n(Pm) (1 e m e M). To construct the isotherm, one should solve the system of eqs 37 for each point Pm. While evaluating the unknown quantities in eqs 36 and 37 by a successive approximation technique, for each point Pm system (37) should be solved repeatedly; therefore it is desirable to use the fastest method for this purpose. For example, one can readily derive the equations for Newtonian iterations. Let θi(0) ) 0 (1 e i e I) be the initial approximation. The increments yi ) θi(s+1) - θi(s) satisfy the system I

yi - θi(s)(1 - θi(s))

I

P/{P + exp[(zi -

∑

Pf∞ i)1

(38)

where s is an iteration counter. After calculating {yi}, we assign: θi(s+1) ) yi + θi(s) (1 e i e I). For a superposition of Langmuirian isotherms, the method ensures exact convergence by the only iteration. Let us consider the problem of evaluation of the expansion coefficients in the equations for an overall isotherm. Suppose first that the matrix {fij} is known. Then our aim is the calculation of the sets {Fi} and {zi}. Their best values should ensure a minimum of the following function of 2 × I variables M

∆({Fi},{zi}) ) M-1

I

Fiθi(Pm,zi)}2 ∑ {1 - [na/n(Pm)] ∑ m)1 i)1

(39)

ni(P) ) lim n(P) ) na Pf∞

or I

Fi ) 1 ∑ i)1

ωijfijθj(s))]} - θi(s) ∑ j)1

1eieI

I

lim

ωijfijyj ) ∑ j)1

(35)

where ni(P) values are the fractional contributions of different sites to the overall isotherm. If experimental information comprises only one isotherm, a simultaneous calculation of the Langmuir constants {A0,i} and the absolute adsorption energies {i} is impossible in principle.13 That is why it is rational to introduce the reduced energies

zi ) -ηi/kT ) -i/kT + ln A0,i Besides, to simplify the notation it is convenient to represent the lateral energies in the reduced form

where ∆1/2 is a mean relative deviation (MRD) between a theoretical isotherm and an experimental one. If the value of na is indeterminate, only the magnitudes {naFi} can be calculated as a result of minimizing function (39). Thereupon, na results from the summation na ) ∑i)1InaFi. This variant of minimization actually means an extrapolation of an isotherm toward P f ∞. However, as mentioned above, such extrapolation can lead to erroneous na values. Therefore it is preferable to evaluate na by an independent method, e.g., by using the BET theory (see Figure 3). Then one may take into account normalization condition (35) by a penalty function method.35 Assume that Mp points in the isotherm are situated in the region of the plateau (see Figure 3). If the P values at those points are very large, one may accept for this region: n(Pm)/na ≈ 1, θi(Pm,zi) ≈ 1 (M + 1 e m e M + Mp, 1 e i e I). By adding function (39) to the weighted square of the error corresponding to these points, we obtain

∆R1({Fi},{zi}) ) ∆({Fi},{zi}) + M+Mp

ωij ) Ziwij/kT

R1Mp-1

Now, the isotherm given by eqs 12 and 13 can be rewritten as

∑ m)M+1

I

{1 - [na/n(Pm)]

Fiθi(Pm,zi)}2 ≈ ∑ i)1 I

∆({Fi},{zi}) + R1(1 -

Fi)2 ∑ i)1

I

n(P) ) na

Fiθi(P,zi) ∑ i)1

(36)

where R1 is the weighting coefficient (or the parameter of

(37)

(32) Gregg, S. J.; Sing, K. S. W. Adsorption, Surface Area and Porosity; Academic Press: London, 1982. (33) Ross, S.; Morrison, I. D. Surf. Sci. 1975, 52, 103. (34) House,W. A. J. Colloid Interface Sci. 1978, 67, 166. (35) Mamleev V. Sh.; Bekturov, E. A. Langmuir 1996, 12, 441.

I

θi(P,zi) ) P/{P + exp[(zi 1eieI

ωijfijθj)]} ∑ j)1

Surface Topography

Langmuir, Vol. 18, No. 10, 2002 3971

penalty). With an increase of R1, the statistical contribution of the points located on the plateau grows. In particular, if R1 f ∞, the normalization condition (35) is fulfilled with absolute accuracy. Depending on the available information concerning a surface, some additional constraints can be included in the objective function. For example, if the k*th and j*th types of adsorption sites belong to certain crystal face, it is natural to require that the magnitudes Fk* and Fj* would satisfy the preset proportion

Fk* ) νFj*

(40)

where ν is a known coefficient taken from crystallographic information. Constraint (40), as well as normalization condition (35), can be implemented by the penalty function method, namely

∆R1,R2({Fi},{zi}) ) ∆R1({Fi},{zi}) + R2(Fk* - νFj*)2

express the values related to the Jth iteration without the superscript “J” for brevity. At all iterations the overall isotherm is linearized as I

Θ(J+1) ≈ Θ(J) +

(1 e i e I)

(41)

36

The barrier function method seems to be the most simple and effective one in order to allow for constraint (41). In applying this method, one should insert an additional term in the equation for the objective function, that is

Βi ) (∂Θ(J+1)/∂βi)|βif0,γif0,

[ln(Fi) + ∑ i)1 2

D] (42) where Rs is a small positive parameter (Rs ∼ 10-9-10-6), D is a positive constant (D ∼ 10). One can see that the minimum of function (42) is ensured with the positive {Fi}. At the limit Rs f 0, the minimums of the functions ∆Rs,R1,R2 and ∆R1,R2 coincide; hence, the problem of the conditional minimization of ∆R1,R2 is equivalent to that of the unconstrained minimization of ∆Rs,R1,R2. The substitution

Fi ) exp(xi)

(1 e i e I)

together with the barrier additive in eq 42 immediately restricts the region of the search to the positive {Fi} components. So, one should solve the following optimization problem

∆Rs,R1,R2({xi},{zi}) w min

(43)

Γi ) (∂Θ(J+1)/∂γi)|βif0,γif0

The explicit expression for Βi is of the simple form

Βi ) xi exp(xi)θi

(1 e i e I)

Computation of {Γi} is a more difficult procedure. Expression (36) for the overall isotherm leads to the equation I

Γi ) (∂Θ(J+1)/∂γi)|βif0,γif0 )

exp(xj) dij ∑ j)1

(46)

where dij ) (∂θj(J+1)/∂γi)|γif0. By differentiating both parts of each of eqs 37, one can find the system of equations for computation of each ith line of the matrix {dij}. The system is of the form I

Akjdij ) bk ∑ j)1

I

∆Rs,R1,R2({Fi},{zi}) ) ∆R1,R2({Fi},{zi}) + Rs

γiΓi ∑ i)1

where

where R2 is the large positive penalty parameter. Finally, while solving the inverse problem it is necessary to provide for the natural condition of non-negativity

Fi g 0

∑ i)1

I

βiΒi +

(1 e k e I)

(47)

where

Akj ) -θk(1 - θk)ωkjfkj(k*j), bk ) 0(k*i),

Ajj ) 1 - θj(1 - θj)ωjjfjj

bi ) -ziθi(1 - θi)

After solving the system of linear eqs 47, only one component Γi can be calculated by formula (46). Thus, to find all the components Γi (1 e i e I), system (47) needs to be solved I times. By using the linearization, one may express explicitly the approximate value of ∆Rs,R1,R2 as the following function of 2 × I variables {βi} and {γi}: M

∆Rs,R1,R2({βi},{γi}) ≈ ∆*({βi},{γi}) ) M-1 I

[na/n(Pm)] [Θ(J)(Pm) +

∑ {1 -

m)1

I

βiΒi + ∑γiΓi]}2 + R1[1 ∑ i)1 i)1

I

In actualization of the minimization (43), one can encounter serious complications. For example, all attempts to use the direct Newtonian method can end with failure. However, a few simple ideas help to achieve success rather quickly. First of all, a scaling of variables needs to be envisaged in the algorithm. It can be made by replacing the absolute increments by the relative ones, namely

xi(J+1) ) (1 + βi)xi(J) ) (1 + βi)xi

(1 e i e I)

(44)

zi(J+1) ) (1 + γi)zi(J) ) (1 + γi)zi

(1 e i e I)

(45)

where J is an iteration counter; hereinafter, we shall

(1 + βixi) exp(xi)]2 + R2[(1 + βk*xk*) exp(xk*k) ∑ i)1 I

ν(1 + βj*xj*) exp(xj*)]2 + Rs

[(1 + βi)xi + D]2 ∑ i)1

(48)

The quantity ∆*({βi},{γi}) is the quadratic function of 2 × I variables {βi} and {γi}, therefore it has a unique minimum. However, due to errors of linearization in eq 48, a decrease of the magnitude ∆* at the iterations does not guarantee a decrease of ∆Rs,R1,R2. Therefore, it is necessary to perform an additional regularization, which (36) Mamleev V. Sh.; Bekturov, E. A. Langmuir 1996, 12, 3630.

3972

Langmuir, Vol. 18, No. 10, 2002

Villie´ ras et al.

is achieved by introducing the additional summands

∆Rβ,Rγ*({βi},{γi}) ) ∆*({βi},{γi}) + I

Rβ

∑ i)1

I

βi2 + Rγ

γi2 ∑ i)1

(49)

where Rβ and Rγ are the positive regularization parameters. If Rβ f ∞ and Rγ f ∞, the increments {βi} and {γi} tend to zeros; therefore a choice of suitable values of Rβ and Rγ ensures a high accuracy of the linearization of the expressions included in ∆Rs,R1,R2. As a consequence, a monotonic convergence of the function ∆Rs,R1,R2 to its minimum is reached by iteration. At all iterations the increments are calculated by solving the system of the linear equations resulting from the extremum conditions

∆Rβ,Rγ*({βi},{γi})/∂βi ) 0,

∂∆Rβ,Rγ*({βi},{γi})/∂γi ) 0

1eieI Upon solving this system, the magnitudes for the next iteration are computed by eqs 44 and 45. We note that by choosing very large values of multipliers Rβ and Rγ in eq 49, one can organize the minimization of the function ∆Rs,R1,R2 either relative to {xi} (small Rβ and large Rγ, e.g., Rβ ) 0.001, Rγ ) 1010) or relative to {zi} (large Rβ and small Rγ, e.g., Rβ ) 1010, Rγ ) 0.001). Thus, the special cases of minimization are carried out within the framework of the common algorithm. The Rs parameter of the barrier function can be taken equal to zero in all cases, where the absolute minimum of ∆R1,R2 is reached with positive values of the components {Fi}. It holds true, as a rule, with small I (I ∼ 1-5). In particular, the calculations described in the next sections were carried out with Rs ) 0. The parameter Rs needs to be assigned as some small constant, e.g., Rs ) 10-6 only when searching for degenerated distributions. The question remaining open is the problem of evaluating the set {fij}. According to the physical sense the elements of this matrix are connected to quotas of various adsorption sites on a surface, that is

fij ) fij({Fi}) ) fij({xi})

(1 e i e I, 1 e j e I)

(50)

Moreover, the functional relationships (50) determine the model used for a surface and are of great importance in the analysis of energy heterogeneity. It is obvious also, if the elements of the matrix {fij} are not varied in iterations, the model loses self-consistency. With an unknown matrix {fij} the optimization problem (43) can be formulated as follows:

∆Rs,R1,R2({xi},{zi},{fij({xi})}) w min

(51)

Although, as before, we have the minimization problem of the function of 2 × I variables, the algorithm of its solution is substantially more complicated. One rational idea consists of the use of simple iterations. Let {xji} be an arbitrary vector. By using this vector one can determine the matrix {fij({xji})}. Thereupon, applying the above algorithm, it is possible to find the components {xi*} and {zi*} satisfying the minimum:

∆Rs,R1,R2({xi*},{zi*},{fij({xji})}) ) min

(52)

Having made the assignment xji ) xi* (1 e i e I), one can repeat the minimization (52).

With small values of the matrix elements {ωij}, the overall isotherm is close to a superposition of Langmuirian isotherms. In such a case one may expect the convergence of the iterations that lead to a self-consistent solution, because constraint (50) weakly affects the results of the optimization. However, with large values of the elements {ωij} the convergence of the simple iterations is hardly possible. We note that in a general case it is wrong to minimize the function ∆Rs,R1,R2 ({xi},{zi},{fij({xji})}) with respect to 3 × I variables and consider that the sets {xi} and {xji} are mutually independent, since the minimum of the given function can be reached in the absence of the selfconsistency. Only in the case of absolutely exact model the global minimum of ∆Rs,R1,R2 is reached with identical {xi} and {xji}. Although self-consistency leads to an increase in the minimum of the function ∆Rs,R1,R2 ({xi},{zi},{fij({xji})}), it should be provided in the algorithm, since otherwise the model used loses physical sense. Thus, the only reliable variant for solving the problem consists of the use of direct minimization (51) relative to 2 × I variables. The algorithm developed is realized as follows. For each set of the variables {xji} it searches for an optimum set {zi*} of the reduced energies {zi} for the function ∆Rs,R1,R2, that is

φ({xji}) ) ∆Rs,R1,R2({xji},{zi*},{fij({xji})}) ) min

(53)

Minimum (53) is calculated by the above algorithm with Rβ ) 1010, Rγ ) 0.001. The value of φ({xji}) depends only upon the set of the variables {xji}. The minimization of φ({xji}) is realized by the Fletcher-Powell method.37 Note that adsorption sites can be numbered arbitrarily. The number of the variants in the indexing is equal to I!. Hence, the objective function possesses I! local minimums. It is obvious that the optimization should be realized only in a neighborhood of one of those minimums. In other words, the indexing of sites should not vary in the course of the iterations. In the algorithm this is ensured as follows. The sites are numbered in ascending order of the reduced energies, namely

z1 < z2 < ... < zI

(54)

If sequence (54) is violated during the iterations, the optimization process of the reduced energies ceases (Rγ ) 1010), and the algorithm optimizes only the set {xji}. At the next step of the iterations the algorithm attempts to optimize the reduced energies, and sequence (54) is checked again. Practical computations show that the algorithm converges to an exact minimum with preserved sequence (54), even if a very rough initial approximation is chosen. In the algorithm considered the magnitudes {ωij} are supposed to be preset. If parameters {ωij} are also variables, the problem of minimization relative to 2 × I variables {xi} and {zi} turns into the problem of minimization with respect to I(2+I) variables: {xi}, {zi}, and {ωij}. In a general case, such a problem is a multiextremal one. One may hope for its solution only if the matrix {fij} does not vary during the iterations. In other words, the evaluation of the set {ωij} is possible only within the framework of the patchwise model, where the matrixes {ωij} and {fij} are of the diagonal form. In this case the Fletcher-Powell method searches for the set {ωii}, while the sets {xi}, {zi} are corrected by the above algorithm. As (37) Fletcher, R.; Powell, M. J. D. Comput. J. 1963, 6, 163.

Surface Topography

Langmuir, Vol. 18, No. 10, 2002 3973

shown in the next section, the optimization relative to 3 × I variables leads to a correct solution, although when evaluating {xi}, {zi}, and {ωii} from real (experimental) data the solution is unlikely to have a high accuracy. The greater the number of target values is, the less is the accuracy of their evaluation. 5. Testing the Algorithm For the sake of better visualization of approximation errors it is convenient to represent in the figures not the isotherms but their derivatives with respect to the logarithm of pressure. The numerical differentiation does not entail any complications and can be carried out with high accuracy. Furthermore, it is desirable to inspect the errors caused by different terms in the decomposition of an overall isotherm. Let us assume that I functions ni*(P), which describe partial occupancies of different sites, really exist in nature. If a model were exact and the functions ni*(P) were known, the parameters of the model could be calculated by the minimization M

Figure 5. The overall derivative calculated for the model of patchwise heterogeneity with ωii ) 3.9 (1 e i e 4).

The results displayed in Figure 5 correspond to the derivative generated for the model of a patchwise surface with the following isotherm:

I

{∑[1 - ni*(Pj)/ni(Pj)]2} ) min ∑ j)1 i)1

4

(55)

under the condition

n(P) )

0.25θi(P) ∑ i)1

θi(P) ) P/{P + exp[-3 + 2(i - 1) - 3.9θi]}

I

ni*(Pj) ) n(Pj) ∑ i)1

(1 e j e M)

(56)

where the set {ni(Pj)} is calculated by using a model. Conditions (55) and (56) allow one to calculate the set {ni*(Pj)} best corresponding to {ni(Pj)}. By introducing the Lagrange multipliers {σj} we have M

∂{

I

{∑[1 - ni*(Pj)/ni(Pj)]2} + ∑ j)1 i)1 M

I

ni*(Pj)}/∂ni*(Pj) ) 0 ∑ ∑ j)1 i)1 σj

(57)

From eqs 56 and 57 one can readily obtain I

∑ nk(Pj) k)1

ni*(Pj) ) (1 - σj/2)ni(Pj) ) n(Pj)ni(Pj)/

(58)

The derivatives for different sites corresponding to the partial contributions calculated via eq 58 are shown as points in the figures below. These points visualize the regions, where the approximation loses accuracy. Such information is useful for an improvement of a model. Both accuracy of solution of the inverse problem and a rate of convergence of the algorithm depend on the nature of input information. A test of the algorithm with small values of the elements {ωij}, when an overall isotherm is close to a superposition of Langmurian isotherms, is not of interest, since such case was investigated earlier.13 The most stringent test relates to the situation, where values of the elements {ωij} are large, e.g., if these border on the critical magnitude 4.0. Note that in those models, where the matrix {fij} is densely completed, the values of the elements {ωij} in the supercritical region can exceed 4.0. Only in the case of a diagonal matrix {fij} is the value of ωii ) 4.0 critical.

As shown in Figure 5 the derivative represents clearly defined peaks. It is not hard to establish the distribution parameters in the case considered. These are calculated immediately by taking the derivative. However, such a case, at least for adsorption of noble gases at the temperature of liquid nitrogen, is not observed experimentally. A derivative usually represents a smooth function exhibiting one or two maximums. This fact allows two hypotheses to be made: either values of the energy of lateral interactions are small or the surfaces investigated contain topographical fragments combining several kinds of adsorption sites. If using the known5 estimation: ω ) 2T3D/T, where T3D is the critical temperature of condensation, we have ω ≈ 3.9 for argon (T3D ) 150.87 K) at the temperature of liquid nitrogen (T = 78 K). Thus, the second hypothesis seems to be more adequate. Suppose that a surface contains two adsorption sites with essentially different adsorption energies 2 . 1. Therefore these sites are occupied almost independently, namely

1 + f11ω11θ1 + f12ω12θ2 ≈ 1 + f11ω11θ1 + f12ω12 ) 1(eff) + ω11(eff)θ1 2 + f21ω21θ1 + f22ω22θ2 ≈ 2 + f22ω22θ2 ) 2(eff) + ω22(eff)θ2 Elementary considerations show that the surface with a complicated topography can be manifested as a patchwise surface with apparent (effective) values of the constants 1(eff), ω11(eff), 2(eff), ω22(eff). The general tendency is a decrease of the apparent values of the energies of lateral interactions due to the appearance of the multipliers fii < 1. The results shown in Figure 6 correspond to the solution of the inverse problem with the energy distribution being identical to that used for generating data in Figure 5, but

3974

Langmuir, Vol. 18, No. 10, 2002

Villie´ ras et al.

Figure 6. The overall derivative and the partial contributions corresponding to different sites calculated for the model of random heterogeneity (ωij ) 3.9, 1 e i e 4, 1 e j e 4) in the case where the energy distribution has the same parameters as the distribution corresponding to Figure 5: (a) initial approximation; (b) final result of optimization.

for the model of random heterogeneity. In the last case the overall isotherm was taken as

Figure 7. The overall derivative and the partial contributions of different sites calculated for the model of patchwise heterogeneity in the case of an energy distribution having the same parameters as the distribution corresponding to Figure 5 but with ωii ) 1.0 (1 e i e 4): (a) initial approximation; (b) result of the optimization after 50 iterations by the Fletcher-Powell algorithm.

In Figure 7 the results for the patchwise model are shown for the following isotherm: 4

n(P) ) 4

n(P) )

0.25θi(P) ∑ i)1

θi(P) ) P/{P + exp[-3 + 2(i - 1) - 1.0θi]} 4

θi(P) ) P/{P + exp[-3 + 2(i - 1) - 3.9

0.25θi(P) ∑ i)1

0.25θj]} ∑ j)1

In accordance with our conclusions, the bell-shaped curves of derivatives for different adsorption sites are dilated. The sets {Fi}, {zi}, and {fij({Fi})} are calculated as a result of the solution of the inverse problem. The rough initial approximation was taken (Figure 6a), nevertheless, the final values are found with excellent accuracy. In other words, the iterative algorithm ensures the convergence of the objective function to the global minimum. The test represented in Figure 6 is very rigid. For example, one may expect that in models containing a sparse matrix {fij}, accuracy and a rate of the convergence of the iterations will be increased. It is of interest to consider a situation with the patchwise and random model giving similar results. The above considerations allow one to conclude that in such a case the energies of lateral interactions in the patchwise model should be diminished approximately by a factor of 4.

Now, in the inverse problem the sets {Fi}, {zi}, and {ωii} are calculated. The initial set {ωii} is taken arbitrarily (Figure 7a). Numerical experiments show that a rate of convergence of the iterations is sharply decreased (3 × I variables) in comparison to the previous example (2 × I variables). However, the algorithm again copes with finding of the global minimum of the objective function. The overall derivative exhibits weak maximums, although the derivatives for the different sites resemble the corresponding curves for the random model. Therefore, to obtain results similar to those represented in Figure 6, one should not only decrease the values of the energies of lateral interactions but also move nearer with each other the values of the adsorption energies. In Figure 8 and Figure 9 the results of the calculations are shown for the distribution, where the interval between neighboring energies is reduced twice in comparison to the previous example. Now, the isotherm is of the form 4

n(P) )

0.25θi(P) ∑ i)1

θi(P) ) P/{P + exp[-1.5 + (i - 1) - 1.0θi]}

Surface Topography

Langmuir, Vol. 18, No. 10, 2002 3975

Figure 8. The overall derivative calculated for the model of patchwise heterogeneity in the case of an energy distribution having the same parameters as the distribution corresponding to Figure 9. Four terms (I ) 4) are used when simulating the isotherm, but solving of the inverse problem is carried out by using only one term (I )1).

Figure 10. The simulated overall derivative that is qualitatively similar to typical derivatives observed under Ar adsorption on clay minerals. The derivative consists of three summands I ) 3, but I ) 2 was used while solving the inverse problem: (a) initial approximation; (b) result of the complete optimization.

Even if I ) 1, the bell (see Figure 8) is approximated with good accuracy. One may conclude intuitively that with I > 1 the minimization problem relative to three manifolds {Fi}, {zi}, and {ωii} is multiextremal. Only with a very accurate initial approximation can the algorithm find the correct solution. While handling experimental information with errors included, an accuracy of evaluation of distribution parameters from isotherms with symmetric derivatives can be very low. For example, after 50 iterations with the initial approximation shown in Figure 9a, the error (∆1/2 × 100%) subsides from 0.236% to 0.419 × 10-2%. Although the latter is much less than an error of any experiment known at present, it is impossible to calculate distribution (Figure 9b) with good accuracy. Finally, we consider the results (see Figure 10 and Figure 11) of handling the artificial isotherm that was generated so that it would correspond to typical isotherms of argon adsorption on surfaces of clay minerals. This isotherm was represented as Figure 9. The overall derivative and the partial contributions of different sites calculated for the model of patchwise heterogeneity in the case of an energy distribution for which the difference between adsorption energies is lower by a factor of 2 than for the distribution corresponding to Figure 7: (a) initial approximation; (b) result of the optimization after 50 iterations by the Fletcher-Powell algorithm.

The overall derivative really represents an ideal bellshaped curve. When the optimization problem is solved in this case, there is the dilemma: the minimization should be carried out either with respect to the parameters of only one of the adsorption sites or with respect to the parameters of all adsorption sites simultaneously.

n(P) ) 0.1θ1(P) + 0.2θ2(P) + 0.7θ3(P) θ1(P) ) P/[P + exp(-4)] θ2(P) ) P/[P + exp(-1.5 - 0.5θ2)] θ3(P) ) P/[P + exp(1.5 - 1.5θ3)] The principal problem for the treatment of experimental information consists of the choice of a suitable number of terms in the isotherm expansion. In particular, in our example the value of I ) 2 does not ensure desirable accuracy (see Figure 10). We note that in the global optimization with I ) 2 the value of ω11 becomes negative;

3976

Langmuir, Vol. 18, No. 10, 2002

Villie´ ras et al.

Table 1. The Results of Solution of the Inverse Problem for the Artificial Isotherm Similar to Many Experimental Isotherms of Argon Adsorption on Clay Minerals exact calculated

F1

z1

ω11

F2

z2

ω22

F3

z3

ω33

0.1 0.999

-4. -4.000

0. -0.001

0.2 0.199

-1.5 -1.51

0.5 0.503

0.7 0.701

1.5 1.495

1.5 1.490

prefer a solution with larger I and positive ωij in comparison to a solution with smaller I and negative ωij. To take into account a priori information concerning positiveness of {ωij}, the same expedient is used in the next section while handling experimental data. One can conclude that the developed38 earlier interactive procedure called DIS (derivative isotherm summation), which has often been used38-40 for approximation of derivatives similar to that shown in Figure 11, can lead to correct results. The main condition for correctness of a solution is not mathematical optimization but some a priori information that is incorporated into the physical model of a surface. 6. Adsorption of Argon on Muscovite at 77.29 K

Figure 11. The overall derivative and the partial contributions corresponding to different sites calculated with I ) 3 for the example corresponding to Figure 10: (a) initial approximation; (b) result of the optimization after 50 iterations by the FletcherPowell algorithm. The exact result is ω11 ) 0.0, ω22 ) 0.5, ω33 ) 1.5.

therefore in the computations displayed in Figure 10b we assigned ω11 ) 0, and the optimization was carried out relative to ω22 only. In our artificial example the exact solution is known, and one can readily understand that strategy of the optimization is to add one more term in the expansion. Indeed, the evaluations made with I ) 3 lead to the correct solution (see Figure 11 and Table 1). However, the question may be raised as to whether such a strategy is suited to real systems. Theoretically, negative elements ωij can correspond to a system where a distance between adsorption sites is less than collision diameter of admolecules. If attractive forces of admolecules with adsorption sites prevail, the admolecules localized in the ith and jth adsorption sites can repulse each other. However, this situation seems to be a rare case in nature. Rather, the admolecules can be likened with hard spheres; i.e., one of the admolecules will extrude another one from this pair of adsorption sites (see Figure 4c). Note that the considered lattice model disregards both a change of distances between adsorption sites and competitive adsorption of admolecules. It fixes the imaginary nodes in accordance with an equilibrium structure of filled monolayer where both a minimum of potential interaction energy with the surface and a minimum of interaction energy between admolecules are simultaneously reached (see Figure 4c). This approximation seems reasonable, and so, at handling a real isotherm, one should

When interpreting experimental data, it is of importance to weigh carefully a priori information about the system under study. One should keep in mind that calculated parameters of energy distributions depend critically upon values of the matrix elements {ωij}. The general tendency13 is an increase in the number of calculated peaks in a distribution with increase in the elements {ωij}. Thus, while approximating an isotherm, it is necessary to solve the dilemma: either the elements {ωij} are large (hence, the distribution is continuous or consists of many closely spaced peaks) or the elements {ωij} are small (hence, the best approximation to the isotherm may be reached with a discrete distribution containing a small number of sharp peaks). In the latter case an increase of I in the algorithm does not influence accuracy, for the calculations result in a degenerated distribution.13 If values of the elements {ωij} are large, but the number I is selected in calculation as being not too large, the approximation loses an accuracy, and the problem becomes multiextremal. Since the surface of muscovite, like the surfaces of other clay minerals, is crystalline, a discrete distribution seems to be the most probable. This fact is accepted a priori in the DIS procedure. In calculations according to the patchwise model the DIS finds a minimum number of delta peaks in a distribution for which a derivative of an isotherm is approximated with prescribed accuracy. The elements of the diagonal matrix {ωii} are corrected during computations. With a minimal acceptable value of I the algorithm finds minimal acceptable values of the elements {ωii}. Such approach is justified by virtue of the hypothesis that, at least at a small coverage of a monolayer, lateral effects are small. The algorithm VADFRE is fulfilled automatically, although in contrast to DIS it requires a good initial approximation; therefore in practical applications it is rational to use both the algorithms sequentially. Let us consider application of VADFRE within the framework of the patchwise model (the diagonal matrixes (38) Villie´ras, F.; Cases, J. M.; Franc¸ ois, M.; Michot, L. J.; Thomas, F. Langmuir 1992, 8, 1789. (39) Villie´ras, F.; Michot, L. J.; Bardot, F.; Cases, J. M.; Franc¸ ois, M.; Rudzin˜ski, W. Langmuir 1997, 13, 1104. (40) Villie´ras, F.; Michot, L. J.; Cases, J. M.; Berend, I.; Bardot, F.; Franc¸ ois, M.; Ge´rard, G.; Yvon, J. In Equilibria and Dynamics of Gas Adsorption on Heterogeneous Solid Surfaces; Rudzin˜ski, W., Steele, W. A., Zgrablich, G., Eds.; Studies in Surface Science and Catalysis 104, Elsevier: Amsterdam, 1997.

Surface Topography

Figure 12. The normalized overall derivative and the partial contributions corresponding to different sites (a) calculated with I ) 4 from the isotherm of Ar adsorption on muscovite for the case of the patchwise model. The enlarged fragment (b) shows an accuracy of the decomposition in the region of low pressure.

{fij} and {ωij}) by using the isotherm of argon adsorption on muscovite as original information. We note that in order to increase accuracy, several initial points on the isotherm were extrapolated by Henry’s law.13 As shown,13 such extrapolation does not significantly influence the energy distributions calculated. The elements ωii ) 1.5 (1 e i e I) were used as the initial approximation in all calculations. DIS shows that the value of ωii = 1.5 ensures good approximation to the isotherm derivative in the region of high pressure. In Figure 12 the results of the calculation of the derivative are shown for I ) 4. The approximation accuracy of the isotherm is rather high (∆1/2 × 100% ) 0.451%), although appreciable disagreement between the calculated and experimental derivative is observed in the region of high pressure. However, the principal reason for rejection of the results shown in Figure 12 is that the value of ω22 is negative; as stated above, negative values of energy of lateral interactions contradict the lattice model. The value of ω11 is unexpectedly large; therefore it is of interest to elucidate what will occur, if ω11 ) 0 is fixed in the calculations (see Figure 13). Now, the approximation accuracy in the region of high pressure (compare Figure 12a with Figure 13a) is increased; however, in the region of low pressure the accuracy diminishes (compare Figure 12b with Figure 13b). The average error (∆1/2 × 100% ) 1.65%) is also appreciably increased. The next trial consists of increasing I from I ) 4 up to I ) 5 and fixing ω11 ) 0 (see Figure 14). The additional term in the expansion does not substantially increase the approximation accuracy (∆1/2 × 100% ) 1.59%). Only while optimizing with five values of {ωii} is it possible to reach high relative precision both at low and at high pressure

Langmuir, Vol. 18, No. 10, 2002 3977

Figure 13. The same as in Figure 12 with ω11 ) 0.

Figure 14. The same as in Figures 12 and 13 with I ) 5 and ω11 ) 0.

(see Figure 15). In this case all the elements {ωii} are non-negative; therefore these do not contradict the model used. Generally, the accuracy of evaluation of the elements {ωii} is not large. In particular, with ω22 ) 0 instead of ω22 ) 0.752, we obtain similar values of the sets {Fi}, {zi} (see Table 2). However, the value of ω11 ) 0, as mentioned above, leads to a considerable loss of accuracy.

3978

Langmuir, Vol. 18, No. 10, 2002

Villie´ ras et al.

Figure 16. The model of accumulation of admolecules around interparticle contacts at a low coverage of an adsorbent.

Figure 15. The same as in Figures 12, 13, and 14 when solving the inverse problem with I ) 5. The average relative error of the isotherm calculation 0.114% even less than estimated experimental error ∼0.2%. Table 2. The Distribution Parameters Calculated for the System Ar-Muscovite with the Optimized Value of ω22 and When Fixing ω22 ) 0 i 1

2

3

Fi 0.022 0.048 0.146 zi -4.76 -3.46 -1.19 Fi 0.022 0.073 0.132 zi -4.66 -3.29 -0.97

4 0.603 0.96 0.600 0.99

5

∑

ω22

∆1/2 × 100%

0.181 1.000 0.752 0.114 2.60 0.173 1.000 0.000 0.222 2.67

The initial fragment of an isotherm (or an isotherm derivative) corresponding to high adsorption energies is of extreme importance for checking a theory, inasmuch as within this fragment the assumption of localized adsorption holds true with good accuracy. Unexpectedly it is found that occupancy of the strongest sites does not obey the Langmuir isotherm. This fact contradicts the common viewpoint according to which lateral interactions can be ignored up to rather a high coverage of a monolayer Θ = 0.6. However, the experiment considered testifies that the strongest sites are grouped together. The authors of the fundamental work28,41 dedicated to Ar adsorption on TiO2 write (p 322, ref 28): “For ideal crystals all sites would have an identical energy or at least there would be a limited number of types of sites. However, for crystals as small as those used in the experiments (0.01 µ) the relative occurrence of corners, edges, irregularities, etc., may be expected to be rather high. Roughly 10% of the surface atoms are at edges or corners, even if the crystals are regular.” Corners and edges are linear formations. In addition to such defects one may imagine spiral dislocations, growth steps, fractures, etc. Adsorption sites within these formations are grouped together in two-dimensional chains. (41) Drain, L. E.; Morrison, J. A. Trans. Faraday Soc. 1952, 48, 316, 840.

Along these chains admolecules can be grouped with the coordination number 2. However, it is hard to explain why the adsorption sites on these chains are the strongest ones, and so an accumulation of the adsorbate around needle-shaped contacts (see Figure 16) between adsorbent particles seems to be the more preferable explanation. In the vicinity of such contacts, admolecules must undergo the strongest attractive force with an adsorbent, inasmuch as the potentials of two contacting particles at the point of the contact are additive. The space between particles near a contact actually represents a narrow slit pore that should be filled by an adsorbate first and foremost, i.e., at lowest pressures. Thus, an adsorbent, which originally behaves as nonporous, after grinding up to a subsieve powder, effectively becomes porous. Around interparticle contacts the admolecules (see Figure 16) can form rings that formally can be considered as a one-dimensional lattice with coordination number 2. However, one cannot exclude a formation of nested rings of admolecules, and so simplified lattice models may incorrectly describe adsorption at a low pressure. Fortunately, the fraction of such “abnormal” adsorption sites is not too large (=2% of common specific surface). Only very precise and laborious experiments such as those used in this work can detect these sites. At the same time one should bear in mind that their occurrence being negligible in adsorption isotherms can predetermine other thermodynamic dependencies at a low coverage, such as isosteric heats of adsorption. The next sites in descending order of adsorption energy  (in ascending order of the reduced energy z) can be related to point vacancy defects that may be expected to be really situated far enough from each other (ω22 ) 0) in space. The final conclusions, however, can be made only after a systematic study of a number of systems. At the same time the results obtained, even for one system, give rise to serious suspicions concerning the correctness of the use of some simplified models of a surface, such as Bernal’s model6,7 assuming an amorphous state of surfaces. For a correct explanation of adsorption phenomena on crystals, it is necessary to take into account mechanisms of formation of cleavages, growth steps, linear and point defects, etc. Random-number generators can hardly simulate statistics for such structures. Apparently, to a large extent these statistics are defined both by the nature of the crystals and by the conditions of their preparation. 7. Conclusions The algorithm developed computes stably not only fractions of different adsorption sites and the corresponding adsorption energies but also matrix elements reflecting correlation in a mutual disposition of these sites. To study an energy distribution, some relationships between these

Surface Topography

matrix elements and the fractions of adsorption sites should be preset in the model. This, in turn, requires an adequate hypothesis concerning crystallography of a specified surface. The most unexpected conclusion in this work is that the local Langmuir isotherm used in the overwhelming majority of publications contradicts the experimental data. In other words, a model that disregards lateral interactions does not ensure approximation of experimental data with required accuracy. Under small pressures, corresponding to large adsorption energy, purely localized adsorption must occur. Therefore the initial fragment of an isotherm is of extreme importance for checking the legitimacy of the assumptions accepted in the lattice model. Unfortunately, even at a very small pressure (∼10-3 Torr) lateral effects complicate the process. Their detection at small pressures can be explained by the grouping of the strongest sites into linear formations such as edges, growth steps, spiral dislocations, etc. However, an accumulation of the adsorbate around needle-shaped interparticle contacts seems to be the most probable. The fraction of these formations is relatively small; however, a precise experimental technique allows their content to be estimated (=2% of common specific surface). By use of the algorithm developed, the energies of lateral interactions are evaluated within the framework of the patchwise model, although the accuracy of such evaluation is rather low. The patchwise model unrealistically reflects crystallography of surfaces, but with empirically fitted energies of lateral interactions it can predict a content of different sites. Such possibility is explained by the fact

Langmuir, Vol. 18, No. 10, 2002 3979

that the derivatives of local isotherms with respect to the logarithm of pressure conserve the similar bell-shaped form for different variants of the matrix of conditional probabilities. It is necessary to keep in mind that the lattice model itself is rather imperfect. The energies of lateral interactions included in the model are semiempirical constants rather than strict physical magnitudes. However, phenomenologically the lattice model is correct, therefore, “calibrating” the model by standard systems, one can find the correspondence between values of real physical constants (e.g., in the Lennard-Jones potential) and those used in the lattice model. In such a case the model acquires predictive power and becomes a useful tool for study of surface structure. Comparison of calculated energy distributions with an identical adsorbate for different surfaces or with different adsorbates (noble gases, nitrogen, and methane) for an identical surface is of great interest for continuation of the work. Acknowledgment. This work was implemented according to the project INTAS-2000, Ref. No. 00-505. Authors are grateful to Professor Yuri Tovbin (Karpov Institute, Moscow, Russia) for consultations within the framework of the joint INTAS project. Vadim Mamleev is grateful to Ministry of National Education, Research, and Technology of France for the grant making his visit to Professor Cases’ laboratory possible. LA0116387