Article pubs.acs.org/Langmuir
Thermodynamic Analysis of Ordered and Disordered Monolayer of Argon Adsorption on Graphite E. A. Ustinov†,‡ and D. D. Do*,‡ †
Ioffe Physical Technical Institute, 26 Polytechnicheskaya, St. Petersburg 194021, Russia School of Chemical Engineering, University of Queensland, St. Lucia, QLD 4021, Australia
‡
ABSTRACT: We presented a detailed thermodynamic analysis of argon adsorption on a graphitized carbon black with a kinetic Monte Carlo scheme. In this study, we particularly paid attention to the formation of a hexagonal two-dimensional molecular layer on a graphite surface and discuss conditions of its stability and thermodynamic properties of the adsorbed phase as a function of loading. It is found that the simulation results are substantially affected by the dimensions of the simulation box when the monolayer forms a hexagonal ordered structure. This is due to the fact that the lattice constant is constrained by the dimensions of the surface. To circumvent this, we presented a thermodynamic technique, which allows for the variation of the box size as a function of loading, to determine the “intrinsic” lattice constant (rather than apparent average value because of the fixed dimensions of the simulation box) and the thermodynamic functions for the adsorbed phase: the Helmholtz free energy, the chemical potential, and the surface tension. The tangential and normal pressures as a function of the distance from the surface are also discussed.
1. INTRODUCTION Understanding the adsorption of gases in a confined space is always preceded by a thorough study of adsorption on an open surface that has the same chemistry as the walls of the confined space. This allows us to investigate the effects of confinement on adsorption isotherm and isosteric heat of adsorption, and it has been the practice that is adopted by researchers both experimentally and theoretically.1−7 On the theoretical front, molecular simulation has been increasingly used because it can be readily executed on desktop computers. A development of simulation techniques and methodology to understand better adsorption phenomena often forces researchers to reconsider seemingly simple systems to analyze those properties that were previously overlooked or did not attract enough attention. For example, the heat spike observed calorimetrically in argon and nitrogen adsorption on graphitized carbon black8,9 at the completion of the first layer is now attributed to the formation of the 2D hexagonal lattice; this was not unambiguously reproduced by molecular simulations for about 30 years until recent work based on a refined Monte Carlo (MC) method.10,11 With regard to the crystalline two-dimensional © 2012 American Chemical Society
molecular layer, it is now realized that the properties of the monolayer are highly affected by the size and shape of the simulation box, a serious restriction due to the finite size used in the simulation. If the graphite surface is modeled as a square in the XY directions and periodic conditions are applied at the x and y boundaries, then the perfect hexagonal lattice cannot be modeled without its asymmetry and some dislocations or gaps along the sides of the box. This difficulty can be resolved, for example, by taking a rhombic cell in the XY plane; however, it does not eliminate the correlation of the lattice constant and the box size unless the box is extremely large. Such a correlation appears because each side of the rhomb equals an integer number of the lattice constant, and, therefore, a relatively small change in the area of the XY plane does not change the number of lattice sites in the monolayer, but rather it forces a change in the distance between neighboring sites. This means that there is infinite number of adsorption isotherms at the same temperReceived: March 30, 2012 Revised: May 28, 2012 Published: May 31, 2012 9543
dx.doi.org/10.1021/la301328x | Langmuir 2012, 28, 9543−9553
Langmuir
Article
ature that can be simulated using an arbitrary box size. This then raises the question of how the lattice constant depends on the loading and temperature and, consequently, how the box size can be adjusted in the course of simulation to avoid the artifacts of dislocation or gaps between domains. This question can be resolved with a careful thermodynamic analysis, which is one of the main objectives of this Article. The order−disorder transition in the monolayer adsorbed on a homogeneous surface is a very complicated issue, which has attracted the attention of researchers for many decades. Lots of effort was directed to the analysis of melting in the adsorbed phase.12−20 In most cases, the order−disorder transition is coupled by a registry and orientation transition if the molecules are not spherical.19,20 At a low temperature and loading that is below monolayer coverage, the adsorbed phase consists of clusters that can coalesce, grow, rotate, solidify, or melt. At a high multilayer loading, the structure of each layer is affected by that of neighboring layers. In view of such interplay of different factors, we find it reasonable to simplify the task to elucidate the net effect of ordering transition resulted from intermolecular interaction in an external potential exerted by unstructured homogeneous and smooth surface. In this Article, we focus on the ordering transition in the monolayer when it uniformly covers the surface, and this transition is manifested as a spike on the differential heat of adsorption as a function of loading. For this reason, the melting of clusters, registry, and orientation ordering is out of the scope of the present study. It was shown recently21 that equilibrium systems can be effectively simulated using a kinetic Monte Carlo (kMC) method.22−25 There are several advantages over the conventional Metropolis scheme realizing the MC method.26 First, there are no discarded trial moves in the kMC, which makes this method very efficient in dealing with a dense phase, where the acceptance ratio in the Metropolis scheme is so small, which results in very small maximum displacement, typically of the order of 10% of the collision diameter. This leads to a strong correlation among configurations and therefore requires additional computational time. Second, the pressure and density in the rarefied regions of an inhomogeneous system can be determined accurately with kMC, and the same cannot be said for the conventional MC because of the poor sampling in those regions. In kMC, a molecule in trial move is chosen with the Rosenbluth scheme, which accounts for the interaction potentials of all molecules, in the manner that the probability of selecting a molecule is larger for the rarefied section,21 which allows for reliable determination of the bulk density and pressure directly without invoking the chemical potential. Another advantage of the kMC is that the determination of the chemical potential is built in the simulation scheme, without recourse to a separate scheme, such as the Widom scheme implemented in the conventional MC simulation. We show in this Article that the kMC provides a very robust analysis of the monolayer order−disorder transition in the system Ar− graphitized carbon black.
consists of interaction potential energy of molecule i with other molecules in the system and the external field:
2. THEORY The main basis of the kinetic Monte Carlo scheme is the residence time of a molecule i proportional to the Boltzmann factor:
where μ°(T) is kBT ln(Λ ) and Λ is the de Broglie thermal wavelength. Equation 8 was derived as an extension of the ideal gas equation to nonideal conditions21 and was shown to be exactly in agreement with the chemical potential determined with the Widom scheme in the framework of the Metropolis MC. Nevertheless, it should be noted that eq 8 leads to the true chemical potential in the limiting case of a sufficiently large system. The physical meaning of μN is the increment of the
τi = exp( −ui /kBT )
N
∑ φi ,j + Vi
ui =
j=1
(2)
j≠i
Here, N is the total number of molecules in the system, φi,j is the interaction potential of molecule i and molecule j, and Vi is the external potential exerted onto molecule i. The inverse of the residence time, νi = 1/τi, which is called the mobility rate in the literature, is proportional to the probability of molecule i to move from its position per unit time. Therefore, the probability of the event that one of the N molecules leaves its position per unit time is proportional to N
R=
∑ νi
(3)
i=1
The average lifetime of a configuration j is then Δτj = 1/Rj, with the term “average” to mean that it is the average over the large number of replicas of the same total rate R. The actual lifetime of the configuration j is given by Δτj = Δτjln(1/p)
(4)
where p (0 < p < 1) is a random number. Equation 4 corresponds to the Poisson distribution for the time step.24 To choose a molecule to be displaced, we consider the following partial sum: k
Rk =
∑ νi
(5)
i=1
The value Rk/R (k = 1,..., N) is between zero and unity and represents a probability of the event that the first molecule that will leave its position is a molecule having number less than k. Therefore, the molecule to be displaced (molecule k) is chosen to satisfy the following inequality: R k − 1 ≤ pR < R k
(6)
where p is the random number (0 < p < 1). The new position of the molecule is chosen randomly and uniformly over the volume V of the system. The time average of any thermodynamic function, A, can be evaluated as: A̅ =
1 τ
M
∑ AjΔτj (7)
j=1
Here, τ is the sum of Δτj and M is the number of time steps, or, which is the same, the number of configurations. 2.1. Chemical Potential. In our previous work,21 we derived the following expression for the chemical potential: μN = μ°(T ) + kBT ln(N /V ) − kBT ln(Nτ /M )
(8)
3
(1)
where ui is the potential energy of the molecule, T is the temperature, and kB is the Boltzmann constant. The potential ui 9544
dx.doi.org/10.1021/la301328x | Langmuir 2012, 28, 9543−9553
Langmuir
Article
Helmholtz free energy of the canonical system, that is, FN − FN−1, and it is defined for a closed system containing N molecules and must be distinguished from the chemical potential μ of the grand canonical ensemble of systems, each of which has N molecules on average. In the limiting case of a very large system, there is no numerical difference between μ and μN. According to the grand canonical distribution, the average number of molecules in the system corresponding to the chemical potential μ is N=
qst = −
(9)
PT(z) = ρ(z)kBT −
where j
Fj =
∑ μi i=1
(10)
N−1
A T (z ) =
ξT (z , zi , zj) =
2rij
drij
ξT(z , zi , zj) (16)
1 [H(z − zi) − H(z − zj)] |zij|
(17)
Here, xij, yij are projections of the separation distance rij between molecules i and j onto the corresponding axes. H(x) is the Heaviside function. With regard to the normal pressure, we distinguished two expressions. The first one, PN(z), accounts for all forces acting onto a molecule including the force exerted by the solid. In this case, the normal pressure does not change with the distance from the surface and equals the bulk pressure, which is negligibly small as compared to the tangential pressure in the adsorbed subcritical phase (in the case of adsorption on the flat surface). The second expression, P*N(z), accounts for only intermolecular interactions, which is more informative because it is closely related to the density distribution and reflects deviations in the behavior of the highly compressed gas in the adsorbed phase from that of the ideal gas. We used the following equation for the intrinsic normal pressure:
(11)
P*N (z) = ρ(z)kBT −
AN̅ (z) Sxy
(18)
where A̅ N is the time average of N−1
AN (z) =
(12)
N
zij dφij
∑ ∑ i=1 j=i+1
Having taken the derivatives on the bottom of eq 12 with respect to the chemical potential μ and accounting for the expressions 9 and 11 for N and A̅ , respectively, we obtain: ∑j j(Aj − A̅ )exp[(jμ − Fj)/(kBT )] dA̅ = ∑j j(j − N )exp[(jμ − Fj)/(kBT )] dN
∑ ∑
and
For a very large system involving say a million molecules, only one term would contribute to sums in the numerator and denominator of eq 11, and if the chemical potential is determined with the system of N molecules, then A̅ =AN. In the case of smaller systems, which are commonly used in simulations, several values around AN contribute to the average value A̅ , which substantially decreases its uncertainty due to fluctuations. Equation 11 allows for evaluation of the derivative of a thermodynamic function with respect to the number of molecules using the following equality:
dA̅ /dμ dA̅ = dN dN /dμ
(15)
xij2 + yij2 dφij
N
i=1 j=i+1
∑j Aj exp[(jμ − Fj)/(kBT )] ∑j exp[(jμ − Fj)/(kBT )]
A̅ T (z) Sxy
Here, ρ(z) is the density of the adsorbed phase at a distance z from the surface, and Sxy is the area of the simulation cell in the XY plane parallel to the surface. A̅ T is the time average of AT, which is given by28
Therefore, at a given number of molecules N, which is proportional to the amount adsorbed, the chemical potential is the root of eq 9. The set of μi for all integers i determined by the canonical simulation defines the adsorption isotherm. However, the accuracy of such an isotherm is usually insufficient because of the fluctuations resulting from a relatively small simulation cell. Equation 9 presents an effective way to circumvent this difficulty. In fact, eq 9 defines the smoothed adsorption isotherm in the form of dependence of the true chemical potential on the average number of molecules in the simulation cell. Once the chemical potential μ at an average number N of molecules in the system is determined, any thermodynamic function A corresponding to a grand canonical ensemble at the same number N of molecules can be determined as A̅ =
(14)
where u is the internal (configuration) energy of the adsorbed phase per mole; upg , p, and ρ are the potential energy, pressure, and density of the bulk phase, respectively. 2.2. Pressure in the Adsorbed Phase. The pressure in a nonuniform system is a tensor, and in the case of adsorption of gases on a flat surface it consists of the tangential and normal pressures. For the tangential pressure, we used the Irving−Kirkwood description:28
∑j j exp[(jμ − Fj)/(kBT )] ∑j exp[(jμ − Fj)/(kBT )]
p ∂Nu + ugp + ∂N ρ
rij drij
(19)
In the above equation, one must count only those pair of molecules i and j, for which the vector rij crosses the XY plane at a distance z. In addition to the pressures calculated from the virial route, we can determine a pressure via a thermodynamic route. As is known, the Gibbs free energy of a heterogeneous phase is
(13)
G = F + PV + γSxy
We used eq 13 to determine the derivative of the internal energy with respect to amount adsorbed to calculate the isosteric heat of adsorption:27
(20)
where P (=PN) is the bulk pressure, V is the volume of the system, and γ is the surface tension given by 9545
dx.doi.org/10.1021/la301328x | Langmuir 2012, 28, 9543−9553
Langmuir
Article
Figure 1. Argon adsorption isotherm on graphite at 77 K simulated with kMC in canonical ensemble. The simulation box lengths Lx, Lz are 18.0 and 10σff, respectively. Ly = Lx(√3/2). Symbols are for simulation; line is plotted with eq 9.
γ=
Lz
∫z=0 (PT(z) − PN )dz
columns along the X and Y directions. If the Lx side equals an even number of the lattice constant, the common boundary conditions can be applied. It is evident that the choice of the box size affects the simulation results because it defines the lattice constant. In our case, the lattice constant (the average distance between neighboring molecules in the monolayer) changes from 1.1125 to 1.1469σff. The true lattice constant depends not only on temperature, but also on loading, and its choice is dictated by a thermodynamic constraint. We will discuss it in section 4. The number of kMC steps was not less than 25 × 106 for sampling and not less than 5 × 106 for equilibration. The pair Ar−Ar interaction energy was modeled by the 12-6 Lennard-Jones potential with parameters εff/kB = 118.4 K and σff = 0.3393 nm, which were obtained from the matching of the kMC simulation results against the vapor−liquid interface data.21 The solid−fluid potential was modeled for the lamella structureless graphite structure by the following equation:
(21)
The integral is taken over the normal to the surface from zero to Lz, which is the height of the simulation box. Because SxyLz = V and P = PN, eq 20 can be rearranged as follows: Lz
∫z=0 PT(z)dz = G S− F xy
(22)
The left-hand side of eq 22 is evaluated with eqs 15−17. The right-hand side of the above equation can be determined explicitly accounting for that G = Nμ and eq 10 for the Helmholtz free energy. However, the Helmholtz free energy of the adsorbed phase can be approximated by the following integral: F (N ) =
∫0
N
μ di
(23)
where i is taken to be a continuous variable. The integral can be obtained explicitly using eq 9: ⎛ jμ − Fj ⎞ F = μN − kBT ln ∑ exp⎜ ⎟ ⎝ kBT ⎠ j=0
n ⎡ 10 ⎛ σsf ⎞4 ⎤ 2 ⎛ σsf ⎞ ⎟ ⎟ ⎥ Vsf (z) = 2πεsf ρs Δσsf2 ∑ ⎢ ⎜ −⎜ ⎝ z + kΔ ⎠ ⎝ z + kΔ ⎠ ⎦ 5 ⎣ k=0
N
(24)
(26)
Substitution of F from eq 24 to 22 and accounting for μN = G yields:
Here, ρs (=114 × 10 ) is the number density of the graphite, and Δ (=0.335 nm) is the interlayer spacing. The number of graphene sheets n was set to 10. The solid−fluid potential well depth, εsf, and the solid−fluid collision diameter, σsf, were calculated using the Lorentz−Berthelot mixing rule: 57.58 K and 0.3397 nm for εsf/kB and σsf, respectively.
1 Lz
∫0
Lz
PT(z)dz =
N ⎛ jμ − Fj ⎞ 1 kBT ln ∑ exp⎜ ⎟ V ⎝ kBT ⎠ j=0
27
(25)
The right-hand side of eq 25 is the so-called spreading pressure defined as the integral V−1∫ N0 N dμ. Equation 25 allows one to not only check the thermodynamic consistency of the system, but also to determine the lattice constant of the hexagonally packed molecular layer on the surface. We will consider this in detail in section 4.
4. RESULTS AND DISCUSSION Figure 1 shows a typical argon adsorption isotherm at 77 K calculated in canonical ensemble with kMC. The box length Lx was 18.0σff. For the sake of clarity, Figure 1b presents a segment of the isotherm in the submonolayer region where a transition from low to high density states is observed. It is seen that the isotherm has two distinct sections. The first one (lower branch) corresponds to a disordered structure of the adsorbed phase, and the other section corresponds to a crystalline twodimensional structure of the monolayer. The smoothed red line is plotted with eq 9, which is the grand canonical distribution, with the chemical potentials μi calculated in a canonical ensemble. Values of the Helmholtz free energy Fi were determined with eq 10. The transition zone connecting the two branches of the adsorption isotherm is also given by eq 9 as a segment of the smoothed isotherm.
3. SIMULATION DETAILS Simulations of argon adsorption on graphitized carbon black were performed with the kMC simulation in a rectangular box having dimensions Lz = 10σff and Ly = Lx(√3/2), where σff is the fluid−fluid collision diameter. The ratio Ly/Lx = (√3/2) was taken to allow for formation of the hexagonally packed molecular layer without asymmetry or distortion along the sides of the box, with one symmetry axis of the lattice being parallel to the Lx side. The length of this side was varied in the range from 17.8 to 18.35 collision diameters to ensure 16 rows and 9546
dx.doi.org/10.1021/la301328x | Langmuir 2012, 28, 9543−9553
Langmuir
Article
Figure 2. Density distribution of argon on graphite at 77 K corresponding to points A and B in Figure 1. The amount adsorbed is 12.337 μmol/m2 (a) and 13.365 μmol/m2 (b). The dimensionless density ρσff3 changes from 0 to 15 according to the rainbow.
The transition from the disordered to ordered monolayer structure is confirmed by the density contour plots in Figure 2, which presents the average density distribution in the first monolayer (z < 1.5σff) at the same chemical potential, but at different loadings. The density distributions presented in Figure 2 were obtained by averaging over 25 × 106 kMC steps. The number of molecules that constitute perfect hexagonal lattice at a given box size is 256. Point A of the isotherm in Figure 1b and the corresponding density distribution shown in Figure 2a involve 240 molecules in the simulation box, which is not enough to form a crystalline lattice of 256 molecules. Thus, the monolayer structure of point A is disordered, although there are some small hexagonally packed 2D clusters visible in Figure 2a. Point B in Figure 1b corresponds to 260 molecules in the simulation box, with 256 molecules to form the crystalline first layer (Figure 2b) and the rest in the second and higher layers. The order−disorder transition always appears in experimental adsorption isotherms of argon and nitrogen on a highly homogeneous graphitized carbon black surface at 77 K as a relatively small kink at a relative pressure of about 0.01.8,29−32 Evidence of this transition is also shown in the differential heat curves8 measured calorimetrically as a spike at a loading corresponding to the completion of the first layer. Numerically, it was reliably confirmed only recently.10,11 Our kMC also leads to a very distinct heat spike using eqs 14 and 13 where A is the internal energy. The behavior of the heat of adsorption with loading is shown in Figure 3. The increase of the heat of adsorption from zero coverage to 11 μmol/m2 is due to densification of the adsorbed layer, which results in an increase in the intermolecular interactions. The sharp decrease of the heat of adsorption after passing through the maximum indicates the onset of formation of the second molecular layer. This means that molecules start to adsorb mainly in the second layer where the solid−fluid interaction is markedly weaker, and, for this reason, the heat released is smaller than that during the first layer formation. However, a further densification of the first monolayer continues as the chemical potential is increased, which at a some point leads to its rearrangement to yield a spike in the heat curve versus loading. The heat spike occurs between 12.2 and 13.5 μmol/ m2, which is similar to that observed experimentally,8 and its maximum at 12.9 μmol/m2 corresponds to the inflection point of the red line in Figure 1b. This means that maximum heat
Figure 3. Isosteric heat of Ar adsorption on graphitized carbon black at 77 K. Kinetic Monte Carlo simulation was performed in the box having length Lx = 18σff. Ly = Lx(√3/2).
releases due to the coalescence of 2D clusters in the course of formation of the defect-free hexagonal lattice. It might seem that the problem of the heat spike and the kink on the adsorption isotherm is completely resolved. However, one should bear in mind that once the hexagonally packed monolayer is formed, the lattice constant is kept unchanged with further loading because of the constant box size used in the simulation. However, as we showed in our previous paper,27 the lattice constant decreases with an increase of loading. This feature significantly complicates the interpretation of experimental adsorption isotherms even for the simple argon− graphite system. Thus, it is known that the order−disorder transition viewed in experimental adsorption isotherms is continuous. This phenomenon was called a “degenerate” firstorder transition8,9 and was attributed to a transition between a hypercritical fluid and a solid phase. At the present state of the art, the continuity of the transition is not completely understood, although most researchers attribute this transition as the first-order type. One of the difficulties is that it is hard to determine accurately the position of the inflection point in Figure 1b. This position depends on the lower limit of the upper branch of the adsorption isotherm corresponding to an ordered monolayer and the upper limit of the lower branch corresponding to the completely disordered adsorbed phase. In both cases, we deal with the limit of mechanical stability of the monolayer, which is hard to determine correctly because of the stochastic nature of the system. The order−disorder transition occurs via formation of small clusters, which grow and coalesce 9547
dx.doi.org/10.1021/la301328x | Langmuir 2012, 28, 9543−9553
Langmuir
Article
the equilibrium order−disorder transition pressure gradually shifts to higher loadings, resulting in a continuous transition in a manner that can be regarded as degenerate first order. At this stage, the question is raised about how one can determine thermodynamically the true lattice constant at a given temperature and loading among infinite number of solutions. We consider this issue via analysis of the pressure tensor of the adsorbed phase. 4.1. Pressure Tensor in the Adsorbed Phase. First, we considered the argon adsorption isotherm at 90 K at which the crystalline monolayer does not occur. The simulated adsorption isotherm is shown in Figure 5a, and the corresponding tangential pressure is presented in Figure 5b. The 90 K adsorption isotherm does not show any kink, indicating that there is no crystallization in the monolayer. The pressure calculated via the thermodynamic route (the spreading pressure, solid line in Figure 5b) practically coincides with the tangential pressure calculated with the mechanical (virial) route. The normal pressure accounting for the interaction by the solid is equal to the bulk pressure at any distance from the surface and is negligible as compared to the tangential pressure. Figure 5b, therefore, suggests that the “thermodynamic pressure” is the same as the “mechanical” tangential pressure. However, this was not the case when we dealt with isotherms at lower temperatures at which the ordering transition occurs. Indeed, Figure 6a shows the compressibility PTV/(NkBT) versus loading evaluated thermodynamically and with the virial route at 77 K, where we see that the difference between the two is eminent. The pressures determined from both routes coincide up to the formation of the 2D hexagonally packed monolayer. As expected, in both cases, the compressibility factor tends to unity at zero loading, because the adsorbed phase behaves like an ideal gas. The increase of the amount adsorbed up to 6 μmol/ m2 results in decreasing the adsorbed phase compressibility due to interaction of molecules. Because on average molecules are at a relatively large distance from each other, the potential of such interactions is mainly negative. Further growth of the amount adsorbed leads to the increase of the tangential pressure because repulsive forces in the monolayer are becoming dominant. The tangential and normal pressures versus distance from the surface are shown in Figure 7 at loadings corresponding to completion of one, one and half, and two monolayers. At a side
with loading, and it is not clear whether the transition becomes a vertical isobaric step (true first-order) in the case of an extremely large simulation cell. Even if a perfect hexagonal lattice has formed over the whole graphite surface, the 2Dcompression of the monolayer can result in a continuous transition, rather than a vertical transition of first order. Figure 4 presents some selected argon adsorption isotherms simulated for boxes of different size (which was maintained
Figure 4. Argon adsorption isotherms on graphite at 77 K evaluated with the kMC. The box length Lx (in units of σff) from left to right: 18.3, 18.1, 18.0, 17.9.
constant during the course of simulation) at 77 K. The difference in the transitions and the segment of isotherms after the transition clearly highlights the complexity of the issue. The isotherms were plotted with eq 9 derived from the grand canonical distribution, but the set of μi (i = 1,..., N) was calculated with kMC in a canonical ensemble. As is clearly seen in Figure 4, all isotherms coincide in the low loading region where the monolayer is disordered, but they deviate from each other in the higher loading region where the crystalline monolayer exists. Most importantly, the chemical potential and, consequently, the bulk pressure at which the ordered and the disordered phase coexist depend on the lattice constant (that is proportional to the box length Lx). The decrease of the lattice constant increases the coexisting pressure and the chemical potential in the range of about one kBT. On the other hand, if the rigid boundary conditions are not imposed on the system, the lattice constant decreases with the increase of the amount adsorbed.27 It allows us to conclude that
Figure 5. Argon adsorption isotherms on graphite at 90 K (a) and corresponding pressures in the adsorbed phase versus loading (b). The box length Lx was 18.0σff. Points (×) are from simulation. The line in (a) is plotted with eq 9. In (b), the solid line is for the “thermodynamic pressure” (the spreading pressure) calculated with eq 24. The dashed line is for the tangential pressure calculated with the virial route (eqs 15, 16). 9548
dx.doi.org/10.1021/la301328x | Langmuir 2012, 28, 9543−9553
Langmuir
Article
Figure 6. Compressibility (a) and the surface tension (b) of argon adsorbed on graphite at 77 K. Other conditions are the same as in Figure 5. Points are referred to the Irving−Kirkwood (IK) equation for the tangential pressure. Black lines correspond to the spreading pressure. Red lines are dependences of the compressibility and the surface tension smoothed with the grand canonical distribution (GCD) using eq 11.
Figure 7. Tangential (a) and normal (b) pressures versus distance from the graphite surface in the system argon−graphite at 77 K. Numbers at lines are amounts adsorbed (μmol/m2). The length Lx = 18.0σff.
4.2. Thermodynamic Functions of the Adsorbed Phase. The difference of the “thermodynamic” (spreading) pressure and the “mechanical” pressure derived via the virial route can be attributed to an influence of the periodic boundary conditions on the behavior of the ordered molecular layer. The lattice constant is artificially kept constant during simulations in the box of constant volume, which does not correspond to the real system. Therefore, despite that the simulated system is in equilibrium and thermodynamically consistent, its thermodynamic functions including tangential pressure can differ from those of real systems at the same loading because of the difference in lattice constant, which is expected to decrease in real systems. To resolve this, we resort to a thermodynamic analysis, which requires a thermodynamic functional to be minimal, and in canonical ensemble this functional is the Helmholtz free energy. We achieved this by optimizing the lattice constant at a specified amount adsorbed. It can be done by sequentially adding one molecule to the system and changing the box size to keep the amount adsorbed per unit area constant. Such a procedure allows for a gradual increase of the lattice constant at an imaginary closed system having constant volume and number of molecules at a given temperature. The Helmholtz free energy per one molecule will also change, and its minimum would indicate the true equilibrium of the real system. Let us now consider a two-stages process of adding one molecule to the system and then increasing its volume at constant temperature to maintain the amount adsorbed. In the
length of the simulation box Lx of 18σff, the maximum number of molecules that can reside in the perfect crystalline monolayer is 256. To ensure the formation 2D ordered molecular layer, we used 257 molecules. For the two other cases, the number of molecules in the simulation cell was 384 and 512, respectively. Additionally, we plot in the figure the pressure distribution (the dashed line) at a loading slightly below the complete monolayer coverage (255 molecules in the simulation cell) to highlight the change in the distribution at the point of formation of the ordered hexagonally packed monolayer. As seen in Figure 7a, the 2D ordering is accompanied by an instant rise of the tangential pressure in the first layer. Figure 7b shows the normal pressure that accounts for only intermolecular interactions. In this case, we consider the pressure in the same way as the barometric pressure in the atmosphere, which is a function of the distance from the surface. Similar to the tangential pressure, the highest normal pressure corresponds to the maximum density, so the formation of the second molecular layer appears in the growth of the second pressure peak. Interestingly, there is a plateau in the normal pressure distribution in the region between the first and second layers, with the pressure in this region increasing with loading. It must be remembered that the total normal pressure including the contribution of the solid−fluid interaction is very small and equal to the bulk pressure, which is negligible as compared to the normal pressure when the solid−fluid interaction is excluded. 9549
dx.doi.org/10.1021/la301328x | Langmuir 2012, 28, 9543−9553
Langmuir
Article
first stage, the Helmholtz free energy increases by μ, and in the second stage it decreases by PTΔV, that is: ΔF = μ − PTΔV
The lattice constant versus the amount adsorbed determined with the points of intersection is presented in Figure 9. As expected, the lattice constant decreases with loading.
(27)
Here, ΔV is the volume change. The amount adsorbed is constant if adding a molecule is accompanied by an increase of the surface area of the simulation box by ΔSxy = Sxy/N. Because the box length Lz is constant, ΔV = V/N. Hence, eq 27 can be rewritten as follows:
PV μ ΔF = − T kBT kBT NkBT
(28)
The compressibility factor decreases with an increase of the chemical potential along constant values of amount adsorbed as presented in Figure 8. This means that the Helmholtz free energy change of the system is a monotonic function along the constant loading without minimum, which might seem incorrect.
Figure 9. Dependence of the lattice constant on the amount adsorbed for the system Ar−graphite at 77 K. The dashed line designates the distance of the potential minimum for the LJ potential of an isolated pair of argon molecules.
The lowest point in the curve shown in Figure 9 is calculated for the box length Lx of 17.8σff. Because there are 16 lattice sites along this dimension, the lattice constant is 1.1125σff. This is between the points where the 12-6 LJ potential passes through zero (z = σff) and the minimum (z = 21/6σff). Our tentative analysis seems to suggest that further growth of loading does not change the lattice constant of the first layer, which remains approximately equal to 1.1125σff. Some proof of that is provided in Figure 10 that shows the tangential pressure Figure 8. The change of the compressibility factor with the chemical potential at specified values of amount adsorbed for the system Ar− graphite at 77 K.
However, it must be remembered that keeping constant the amount adsorbed, we model a process referred to as a change of the lattice constant in an imaginary canonical system having the same volume at constant number of molecules. Therefore, we have to find the minimum of the molecular Helmholtz free energy, rather than the total Helmholtz free energy of the system that we use. Let f(N) be the Helmholtz free energy per one molecule at N molecules in our flexible system. Then ΔF = (N + 1)f (N + 1) − Nf (N ) ≅ N Δf + f (N )
(29)
Figure 10. Tangential compressibility factor of argon adsorbed on graphite at 77 K. The simulation box length Lx is 17.8σff. (1) The spreading pressure determined with eq 25. (2) The tangential pressure calculated via virial route with eqs 15−17. The two pressures coincide at the amount adsorbed above 26 μmol/m2.
Substitution of ΔF from eq 29 to 28 yields: N
Δf μ−f PV = − T kBT kBT NkBT
(30)
In the right-hand side of eq 30, the difference (μ − f) is PV/N, where P is the spreading pressure determined with eq 25. The minimum of the specific Helmholtz free energy requires Δf to be zero. This leads us to the conclusion that given the amount adsorbed and the temperature, the lattice constant takes a value, which ensures equality of the tangential pressure determined by the virial equation and that determined thermodynamically (spreading pressure). In other words, at any specified amount adsorbed, the true state of a large adsorption system involving an ordered monolayer corresponds to the point of intersection of “thermodynamic” and “mechanical” branches of the tangential pressure as shown with arrows in Figure 6.
determined via thermodynamic and virial routes with the box length Lx of 17.8σff. In this case, there is no crossing between the two pressures, but rather they merge for loadings greater than 26 μmol/m2. A possible explanation of the lower limit of the lattice constant is that each molecule in 2D hexagonal lattice interacts with all other molecules in this lattice. The distance between a given molecule and any molecule beyond the first coordination sphere is larger than the distance (=21/6σff) at which the LJ potential passes through the minimum. This means that the molecule experiences only attractive forces from all other molecules except its six nearest neighbors, which yields a compression of the lattice and hence a 9550
dx.doi.org/10.1021/la301328x | Langmuir 2012, 28, 9543−9553
Langmuir
Article
volume is changed. This means that we need to interpolate the set of the chemical potential and other thermodynamic functions to reduce them to a constant volume of the system, which contains an integer number of molecules. Finally, we obtain the full set of the thermodynamic functions that can be treated with the grand canonical distribution in a standard way. In Figure 12, we present two simulated adsorption isotherms. The first isotherm (red line and circles) is calculated in the box
decrease of the lattice constant. An estimation of the limiting value of the lattice constant can be obtained by minimization of the interaction potential of the perfect hexagonal lattice per one molecule. This potential versus the lattice constant is shown in Figure 11.
Figure 11. The potential of intermolecular interaction in a 2D hexagonal lattice versus the lattice constant. The dashed vertical line designates the minimum of the potential at a lattice constant of 1.1115σff. The dash−dotted line corresponds to the minimum of the pair LJ potential at a distance 1.1225σff.
Figure 12. Argon adsorption isotherm at 77 K simulated with kMC in the box of constant volume having length Lx = 18.35σff (1, 2) and in the box of variable size (3, 4) after the monolayer has ordered. Lines 2 and 4 are plotted with the grand canonical distribution using eq 9.
The minimum of the interaction potential per molecule in the hexagonal lattice corresponds to the lattice constant of 1.1115σff, which is very close to that in the simulation box of 17.8σff (=1.1125σff). Further contraction of the lattice would lead to a rapid increase of the contribution from repulsive forces, which explains the existence of the lower limit of the lattice constant. 4.3. Thermodynamically Consistent Adsorption Isotherm and Tangential Pressure. Once the dependence of the lattice constant on loading is determined, the adsorption isotherm can be calculated as follows. First, one should take the simulation box that provides the largest lattice constant and carry out simulations up to formation of the 2D hexagonally packed molecular layer. It is then necessary to gradually decrease the box size as the amount adsorbed grows according to the pattern shown in Figure 9. Thus, if the initial box length Lx is taken equal to 18.35σff, meaning that the number of sites in the monolayer is 256, it has to be gradually decreased up to 17.80σff at the amount adsorbed of 26 μmol/m2. The dependence of the lattice constant on the amount adsorbed shown in Figure 9 was approximated by a polynomial, so given the number of lattice site along the box dimension Lx and the number of molecules in the box, one can determine the amount adsorbed and the length Lx corresponding to the chosen number of molecules. The grand canonical ensemble can be represented as a set of systems of equal volume, but the lattice constant in each system depends on the number of molecules in it. The state of each system containing j molecules can be determined in a canonical ensemble. Thus, given j molecules in the system, one can determine the average configuration energy per one molecule, the compressibility factor, and the chemical potential μj. In the case of relatively small systems that we usually use in simulations, it can be achieved by adjusting the box size to obtain the lattice constant, which corresponds to the amount adsorbed at a given number of molecules in the box. However, it should be kept in mind that in this case the amount adsorbed is not proportional to the number of molecules in the system because each time one molecule is added the box
of constant volume (and, consequently, the lattice constant), having length Lx = 18.35σff. The second isotherm (green line and crosses) is calculated accounting for the change of the lattice constant with the amount adsorbed. The isotherms presented in Figure 12 deviate from each other once the hexagonal monolayer has appeared. The adsorption isotherm accounting for the dependence of the lattice constant on the amount adsorbed corresponds to the minimum of the Helmholtz free energy calculated in a canonical ensemble at any loading. As seen in Figure 12, the standard technique of simulation in the box of constant volume results in a lower slope of the isotherm at a bulk pressure exceeding the pressure of the disorder−order transition. The difference in the tangential pressure determined with the proposed method and that calculated in a simulation box of constant volume is very significant as shown in Figure 13, where we plotted the surface tension versus loading calculated with the two schemes. The surface tension is proportional to the average tangential pressure, but contrary to that, does not depend on the normal box length Lz, which is quite convenient. The surface tension determined with the thermodynamic (line 1) and mechanical (line 2) routes using the developed technique now agree very well, which confirms the thermodynamic consistency of our method. At the same time, the surface tension calculated with the standard procedure in a box of constant volume (lines 3 and 4) drastically deviates from the real value. A small increase of the box length Lx by 0.55σff increases the surface tension nearly twice! However, an even more important feature of the standard technique is that the tangential pressure decreases with the increase of the chemical potential starting from the point of ordering transition. This is a clear sign that such a system cannot be stable without imposing the periodic boundary conditions in a relatively small simulation box, which artificially keeps the lattice constant unchanged. Our method overcomes this difficulty, and, as seen in Figure 13, the surface tension (and, consequently, the 9551
dx.doi.org/10.1021/la301328x | Langmuir 2012, 28, 9543−9553
Langmuir
Article
Figure 14. Heat of adsorption in the system argon adsorbed on graphite at 77 K simulated with kMC accounting for the dependence of the lattice constant on loading (−) and in the box of constant volume having length Lx = 18.35σff (− − −) and 17.90σff (− · − · −). Red line is experimental data of Grillet et al.9
Figure 13. Surface tension of argon adsorbed on graphite at 77 K simulated with kMC accounting for the dependence of the lattice constant on loading (1, 2) and in the box of constant volume having length Lx = 18.35σff (3) and 17.80σff (4). Lines 2−4 are plotted with eq 11 using virial route for the tangential pressure (eqs 15−17). Line 1 is plotted with the grand canonical distribution using thermodynamic route (eq 25).
chosen arbitrarily. Otherwise, the results of simulation would be not only incorrect, but they would unpredictably change from simulation to simulation. It should be also kept in mind that an incorrect choice of the box size and shape could lead to a sudden change of a lattice symmetry angle relative to the box sides followed by an instant change in the lattice constant. It would then result in the change of the heat of adsorption, which might be erroneously interpreted as a simulation noise. Figure 14 also shows a comparison of the simulation result with experimental data of Grillet et al.9 The heat spike measured experimentally has the same magnitude as that obtained in our simulation. The latter is slightly shifted toward the lager value of the amount adsorbed, which can be due to an overestimation of the graphite surface area in the experimental measurement using the BET method.
tangential pressure) determined by a virial route is monotonically increasing as a function of loading. 4.4. Thermodynamically Consistent Heat of Adsorption. As is mentioned above, the spike observed in a heat curve just after the formation of ordered monolayer is due to the heat released during two-dimensional crystallization of that layer. For this reason, to reproduce the heat spike correctly in terms of thermodynamics, it is necessary to account for the dependence of the lattice constant on loading. To this end, we calculated the isosteric heat of adsorption with eq 14 in which the partial derivative of the internal energy of the system with respect to the number of molecules N was calculated with eq 13 meaning that A is Nu, where u is the configuration potential per one molecule. The set of uj for all j incremented by one starting from unity was calculated using our new method. First, we determined the sets of uj and μj in a simulation box of variable size accounting for the lattice constant change with the number of molecules in the system. Next, we recalculated these sets by an interpolation to provide integer number of molecules in the system of constant volume. Having done this, we determined the isosteric heat of adsorption with eq 13. The result of application of this method is shown in Figure 14. The heat of adsorption is plotted in Figure 14 along with that calculated in a box of constant volume. The comparison shows that in the latter case the heat curve substantially depends on the box size. In all cases, neglecting the lattice constant−loading dependence leads to underestimation of the isosteric heat of adsorption in the region to the right of the heat spike. This is because in reality the compressing of the 2D ordered layer leads to the decrease of the intermolecular potential energy, which is not accounted for in the case of constant box size. The magnitude of the heat spike is also affected by the box size. It takes the maximum value only in the case when the box size allows the lattice constant to be relatively large at the instant of ordering transition. In all other cases, the magnitude of the heat spike would be underestimated. Thus, in case the box length Lx is underestimated by 0.45σff as compared to the thermodynamically grounded value in the region of the ordering transition, the heat spike is less than the true value by about 1 kJ/mol and shifted toward higher loading by about 1 μmol/m2. This shows that the simulation box size and shape must be by no means
5. CONCLUSION We presented a novel approach based on the kinetic Monte Carlo method to describe two-dimensional ordering transition in molecular layer in the system argon−graphite at 77 K. We showed that at 77 K argon molecules are arranged in the hexagonal lattice whose lattice constant depends on loading. The thermodynamic consistency requires that the box size must be changed with loading in a manner that the lattice constant corresponds to the minimum molecular Helmholtz free energy. This requirement is that the tangential pressure of the adsorbed fluid determined with the mechanical (virial) route is the same as that determined thermodynamically via the difference between the Gibbs and Helmholtz free energies. In the latter case, the tangential pressure is the spreading pressure. This method, accounting for the dependence of the lattice constant on the amount adsorbed, allows us to describe very accurately the spike observed experimentally on the heat curve versus loading as the phenomenon associated with the ordering transition.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Notes
The authors declare no competing financial interest. 9552
dx.doi.org/10.1021/la301328x | Langmuir 2012, 28, 9543−9553
Langmuir
■
Article
ACKNOWLEDGMENTS We acknowledge fruitful discussion with D. Nicholson. This work is supported by the Russian Foundation for Basic Research (project no. 11-03-00129-a). Support from the Australian Research Council is also acknowledged.
■
REFERENCES
(1) Ohba, T.; Murata, K.; Kaneko, K.; Steele, W. A.; Kokai, F.; Takahashi, K.; Kasuja, D.; Yudasaka, M.; Lijima, S. Nano Lett. 2001, 1, 371. (2) Ohba, T.; Kaneko, K. J. Phys. Chem. B 2002, 106, 7171. (3) Kowalczyk, P.; Tanaka, H.; Kaneko, K.; Terzyk, A. P.; Do, D. D. Langmuir 2005, 21, 5639. (4) Kruk, M.; Jaroniec, M. J. Phys. Chem. B 2002, 106, 4732. (5) Kruk, M.; Jaroniec, M. Chem. Mater. 2003, 15, 2942. (6) Jiang, J.; Sandler, S. I. Langmuir 2004, 20, 10910. (7) Miura, K.; Yanazawa, H.; Nakai, K. Adsorption 2007, 13, 139. (8) Rouquerol, J.; Partyka, S.; Rouquerol, F. J. Chem. Soc., Faraday Trans. 1 1977, 73, 306. (9) Grillet, Y.; Rouquerol, F.; Rouquerol, J. J. Colloid Interface Sci. 1979, 70, 239. (10) Wongkoblap, A.; Do, D. D.; Nicholson, D. Phys. Chem. Chem. Phys. 2007, 10, 1106. (11) Fan, Ch.; Birkett, G.; Do, D. D. J. Colloid Interface Sci. 2010, 342, 485. (12) Migone, A. D.; Li, Z. R.; Chan, M. H. W. Phys. Rev. Lett. 1984, 53, 810. (13) Zhu, D. M.; Dash, J. G. Phys. Rev. Lett. 1986, 57, 2959. (14) Morrison, J. A. Pure Appl. Chem. 1987, 59, 7. (15) Pettersen, M. S.; Lysek, M. J.; Goodstein, D. L. Phys. Rev. B 1989, 40, 4938. (16) D’Amico, K. L.; Bohr, J.; Moncton, D. E.; Gibbs, D. Phys. Rev. B 1990, 41, 4368. (17) Day, P.; Lysek, M.; Madrid, M.; Goodstein, D. Phys. Rev. B 1993, 47, 10716. (18) Larese, J. Z.; Zhang, Q. M. Phys. Rev. B 1995, 51, 17023. (19) Kuchta, B.; Etters, R. D. Phys. Rev. B 1987, 36, 3400. (20) Golebiowska, M.; Firley, L.; Kuchta, B.; Fabianski, R. J. Chem. Phys. 2009, 130, 204703. (21) Ustinov, E. A.; Do, D. D. J. Colloid Interface Sci. 2012, 366, 216. (22) Bortz, A. B.; Kalos, M. H.; Lebowitz, J. L. J. Comput. Phys. 1975, 17, 10. (23) Gillespie, D. T. J. Phys. Chem. 1977, 81, 2340. (24) Fichthorn, K. A.; Weinberg, W. H. J. Chem. Phys. 1991, 95, 1090. (25) Battaile, C. C. Comput. Methods Appl. Mech. Eng. 2008, 197, 3386. (26) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (27) Ustinov, E. A.; Do, D. D. J. Chem. Phys. 2012, 136, 134702. (28) Irving, J. H.; Kirkwood, J. G. J. Phys. Chem. 1950, 18, 817. (29) Kruk, M.; Jaroniec, M.; Gadkaree, K. P. J. Colloid Interface Sci. 1997, 192, 250. (30) Kruk, M.; Li, Z.; Jaroniec, M. Langmuir 1999, 15, 1435. (31) Gardner, L.; Kruk, M.; Jaroniec, M. J. Phys. Chem. B 2001, 105, 12516. (32) Olivier, J. P. In Adsorption by Carbons; Bottani, E. J., Tascón, J. M. D., Eds.; Elsevier: Amsterdam, 2008; Chapter 7, p 147.
9553
dx.doi.org/10.1021/la301328x | Langmuir 2012, 28, 9543−9553
