- Email: [email protected]
Studies of the adsorption of N2 on the graphite basal plane by computer simulation
Vol.
Carbon 2.5, No. 1, pp. 1-17, 1987 Printed in Great Britain.
0008-6223187 $3.00 + .Ml 0 1987 Pergamon Journals Ltd.
STUDIES OF THE ADSORPTION OF N2 ON THE GRAPHITE BASAL PLANE BY COMPUTER SIMULATION Department
WILLIAM A. STEELE of Chemistry, The Pennsylvania State University, University Park, PA 16802, U.S.A.
ALEXEI V. VERNOV Department of Physical and Colloid Chemistry, People’s Friendship University, 117302 Moscow, U.S.S.R. and
Department
DOMINIC J. TILDESLEY of Chemistry, The University, Southampton,
England
(Received 1 May 1986) Abstract-A molecular model for the physical adsorption of nitrogen on graphitized carbon black (gcb) has been utilized extensively in computer simulations of this system. After a brief review of the details of the model and of the simulation technique, a number of results produced in the course of these studies are summarized. These include the thermodynamic and structural properties of N, adsorbed on gcb at 73.6 K and at coverages ranging up to -2 l/2 monolayers. Changes in molecular orientation relative to the surface are found to have a significant effect on the lateral interaction in a layer; variations in packing density in the first layer are discussed; and the rate of molecular interchange between first and second layers is characterized. Likely directions for future simulation studies of physisorbed films are proposed. Key Words--Adsorption,
nitrogen, basal plane, computer simulation.
1. INTRODUCTION
For many decades, the characterization of surface area, pore structure and related properties of solid adsorbates has begun with measurement of nitrogen isotherms at temperatures near the normal boiling point of this inert nonpolar adsorbate, often followed by a BET surface area determination or, in the case of porous materials, a calculation of pore volume from either isotherm hysteresis or Dubinin-Radushkevich analysis of the data or both. In spite of the success of these procedures, there is much to learn about the molecular level properties of physisorbed nitrogen, even on homogeneous nonporous solids. Recently, a number of sophisticated experimental techniques have been applied to this problem. These include diffraction of electrons, neutrons and X-rays [l] to characterize the various solid-like ordered layers that form at low temperatures when nitrogen is adsorbed on the exposed basal planes of graphitic carbon. In addition, phase equilibria and the thermodynamic changes associated with two-dimensional (or quasi-two-dimensional) transitions have been obtained from high-precision isotherms and heat capacities [2]. These studies have revealed a rich phase diagram for the N,/gcb (graphitized carbon black) system with much of the detail appearing at T < 50
K. Nevertheless, some interesting and unexpected behavior has been observed at T = 78 K-for example, the presence of an ordered, solid-like film at coverages close to monolayer completion has been observed by neutron diffraction and confirmed by Xray diffraction. An alternative view of the molecular behavior in physisorbed films can be obtained by applying the modern algorithms of computer simulation to these systems [3]. This paper is devoted to a review of some of the information gained in this way. We will concentrate upon the coverage dependence of the properties of physisorbed nitrogen at 73.6 K, particularly at coverages in the dense monolayer to bilayer region. However, before discussing the results of these studies, the technique will be briefly together with a description of the assumptions needed to give a realistic description of the system.
2. SIMULATION METHODOLOGY
The technique used is a well known and wrdely applied method in which it is assumed that molecules obey the laws of classical Newtonian dynamics so that, for example, the center-of-mass position r, and velocity v, of particle i with mass m at a time I t Al
8
W.
A. STEELEet al.
can be obtained from r,(t + At) - r,(t) = v(t) At + correction terms,
(2.1)
v,(t + At) - v,(t) = ; F,(t) At + correction
terms,
(2.2)
where F, is the force on molecule i due to its interactions with other molecules and, in the case of physisorption, with the solid. Molecular trajectories for the system of interest are generated by assuming a realistic intermolecular potential (and molecule-solid interaction) and calculating the force by differentiating these functions. In the case of nonspherical molecules, equations for the angular coordinates and velocity w, must also be solved. For example, w,(t + At) - o,(t) = f N,(t) At + correction
terms,
(2.3)
where I and N, are the moment of inertia and torque, respectively. The calculation of the change in molecular orientation requires the introduction of a convenient set of coordinates. We here use quaternions-their definition and equations of motion are discussed in detail elsewhere [4]. In the computer, the slow step in generating the trajectories is the sum over all pairs of molecules (or sites within the molecules, for the model potentials discussed below) to obtain the force and torque on each molecule in the system. Thus in order to study the system properties over as long a time as possible, one wishes to maximize the timestep size At, consistent with an accurate set of trajectories. This is done by using one of several suggested prescriptions for the correction terms in eqns (2.1)-(2.3)-the algorithm used in the studies discussed here was a predictor-corrector method due to Gear [5]. The simulation is run for a sample of 100-250 molecules which are initially placed in some convenient set of positions, such as an ordered layer next to the surface, and with a random distribution of velocities chosen such that m(v’) = 3kT, I(d)
= 2kT,
(2.4)
where T is the desired temperature and the brackets denote an average over all molecules (and, when appropriate, over time steps as well). In the initial stages of the simulation, the velocities are resealed to keep T constant. For a timestep that is taken to be -1Om’5 sec. in the systems studied here, the systems appear to reach statistical equilibrium after 30006000 time steps. At this point, velocity resealing is usually terminated and the molecular system is allowed to evolve in time as a microcanonical ensemble
(constant particle number, volume and total energy). It is now well known that statistical averages obtained from such simulations conform nicely to the postulates of statistical mechanics, yielding Maxwellian velocity distributions and the expected values for the various fluctuation theorems. One other crucial assumption which must be made for this procedure to be successful arises in connection with the boundary conditions for a sample of (say) 100 molecules placed in a layer near the surface of a solid. We wish this sample to mimic the real adsorbed film containing -lOI molecules while maintaining the “adsorption volume” constant. This is conventionally achieved by surrounding the original sample with images of itself and then allowing the “walls” to be permeable. Thus a trajectory that causes a molecule to pass out through a wall will cause one of its images to pass in through the opposite wall, thus conserving both number and volume. In the adsorption case, these periodic boundary conditions are needed only in the two dimensions parallel to the surface. Other than requiring the box size match the periodicity of the surface so that molecules crossing a boundary suffers no discontinuity in molecule-solid interaction, the implementation is straightforward. As a result, each molecule in the computer box will interact with a mixture of molecules and, if it is near a boundary, images. Evidently, errors can be produced if a molecule interacts with its own image-a problem for small boxes or for long-range interaction laws. Of more importance is the fact that the longrange translational order that develops in solidified films can be materially affected by the periodicity of the images. Thus detailed studies of the thermodynamics of the transitions involving such phases require great care. The boundary condition for motion perpendicular to the solid surface is handled by simply reversing the z-component of those molecules at some large distance from the surface (-10 atomic diameters) if they have been desorbed because of an energetic fluctuation. This seldom happens under the conditions of the simulations performed to date, but a related occurrence which can cause significant problems is the transfer of a molecule from one adsorbed layer to another, particularly when the transfer involves first layer molecules. Since the molecule-solid potential is quite large for first layer molecules, a large change in potential energy will result. Conservation of total energy then requires a significant change in the total kinetic energy. For the restricted number of molecules under study, the consequence is an unwanted large jump in temperature. In the simulations at 73 K discussed here, this problem was avoided by utilizing one [6] of several recent suggestions for a constant temperature algorithm-in essence, velocity scaling is continued for the entire time span of the simulation according to a prescription which approximates [7] contact with a heat reservoir. In specifying the molecular interactions for the NJ GCB system, one makes maximum use of the em-
Adsorption of N2on the graphite basal plane pirical results obtained in simulations of bulk N, and in zero-coverage extrapolations of experimental heat and isotherm data. In particular, it is known that a reasonable compromise between computational simplicity and accuracy for the bulk properties is given by the site-site model for N,-N, interactions [8] in which it is assumed that each molecule contains two sites centered on the atoms and separated by the bond length L = 1.10 A. These interact with the sites on other molecules by Lennard-Jones (LJ) energies plus a Coulombic term that accounts for the electrostatic interaction associated with the N, quadrupole moments. The electrostatic calculation is completed by placing two charges at the center of the molecule which are of opposite sign to the “atomic” charge 4. The molecule quadrupole moment Q is qL?/2, and q can be adjusted to give the desired “effective” Q for a molecule in a dense adsorbed phase. The total N2-NZ interaction U,,(ij) for molecules i and j can now be written as
‘)
where we have switched to reduce units denoted by asterisks, defined for energy by dividing by Ed\, and for distances, by dividing by a = 2.46 A, the length of the lattice-vector for the graphite basal plane. The periodic term in the energy is given to good accuracy
by 2rrAh u;, (7. 2) = * 0, _
where A = a,la, u,* = area of a unit cell of the graphite lattice = J312: g* is the reduced length of a reciprocal lattice vector = 41~/,/3. K, and Kc are modified Bessel functions of the second kind and j”(7) is a periodic function of position T = x,y relative to the surface lattice which is defined to be f(7) = -2[cos
2n(P
+ y*/,:3)
+ cos 27F(X* - y*/\,‘3) f cos 4ay”/,
It is assumed that appropriate units for 4, are used; r,, is the separation between a site in molecule i and one in j, and enn,u,, are the empirical well-depth and size parameters for N-N. Values used in this work are egRlk = 36.4 K and uXg = 3.32 A. In the spirit of the site-site approach, one can devise [9] a molecular version of the atom-atom summation approximation to model N, molecule-carbon black potential. This is known to give a reasonably accurate account of rare gas/GCB interactions [lo], and its extension to the molecular U&i) can be written as
where the sites involved in r,, are nitrogen “atoms” in molecule i interacting with carbon “atoms” in the solid. Well-depth and size parameters that give good agreement with the low coverage experimental data for this system are EJk = 31.9 K and up, = 3.36 A. Equation (2.6) is not convenient for use in computer simulation since it involves a sum over a large number of carbon atoms and two N sites for each molecule. Thus one tranforms [lo] it to
Ah
+ 12.42(z* + 0.994)” _ 2z*’ + 9.66z* + 13.33 8.28 (z* + 1.38)’ +
u;,
(7. z),
(2.7)
3).
(2.9)
The T-dependent terms shown in eqns (2.8) and (2.9) are only the first in an infinite series which correctly represents any property possessing the periodicity of the graphite basal plane. This plane is made up of hexagonal adsorption sites defined by the carbon atoms; the hexagons are spaced by 1 unit along the x direction and by 4312 units in the y direction. Only l/3 of these sites are occupied by NZmolecules in the commensurate film, which thus is described as a (,‘3 x J3) R30” overlayer. Nearest neighbor hexagons cannot both be occupied in such films because of the large size of the Nz (>3.3 A) compared to the site spacing (2.46 A). It is interesting to note that the area per molecule is 15.7 AZ in the commensurate film, a value close to the 16.2 A: taken for the N, area in the BET calculation. It is evident that these model potential functions embody several desirable features. such as the nonsphericity of the NZ molecules. Consequently. the minimum energy configuration for the commensurate layer is the orientationally ordered arrangement shown in Fig. 1. One can readily imagine that this order, both in-plane and out-of-plane, will be lost as temperature is increased; furthermore, changes in surface density will also produce alterations in the angular ordering for the minimum energy structure as well as affecting the temperature dependence of the disordering process. These expectations have been confirmed by simulation, as will be discussed below. Another important feature of the molecule-solid potential is that it includes a reasonable estimate of the variation in energy that occurs as a molecule translates parallel to the surface. Figure 2 shows this variation for a molecule placed parallel to the surface at a distance corresponding to the minimum energy over the center of a hexagon. This variation can give rise to “localized” or “site-wise” adsorption, but no such restrictive approximation need be used here.
W. A. STEELEet at.
**
**
‘W. *.
**
*/ .
\
*.
.
.*
**
.*
**
**
**
**
/.*.**~..~**/ ,*.
.
.
.
.
hi’. :_** .‘,‘* .*,*. .*;. .‘,‘. .
B
* * . .+..
**
.*-**
. .
. .
.*
.*
.a
.
.
*:I:. v-- -_1 * .;,;. .-_+ +-_. .I{. .y_. *I<
Fig. 1. A top view of the minimum knergy configuration for model nitrogen molecules on the basal plane of graphite. Dots indicate carbon atoms. The N2 molecular axes are all parallel to the surface plane, with in-plane orientations indicated by the solid lines. The herringbone pattern comprises two sublattices with a glideline indicated by AB. The orientation angles 4 relative to the glideline are 245”.
The calculated energies of l~alization are rather small, and would hardly be expected to produce sitewise adsorption without significant enhancement from the N2-N2 potentials. Indeed, the successful representation of the site-site summation for the molecule-solid potential by eqns (2.4)-(2.6), which contain a very limited number of periodic terms, is only possible when the periodicity is small. The question of whether the magnitude of the periodicity obtained from this model is correct remains debatable. Physically, it is small because the carbon atoms are densely packed in the graphite basal plane so that a N site interacts simultaneously with a relatively large number of C sites, regardless of its position over the lattice. Although this aspect of the calculation is unIikely to be greatly altered in alternative models, the questionable feature actually arises from the assumption of a symmetric N-C site interaction. In view of the great assymmetry in the carbon-carbon bonding in graphite, it has been suggested [12] that the C sites are not strictly spherical in their polarizabilities and interactions, and as a result, the periodic terms may be as much as twice the size as those obtained from spherical C sites. Consequently, some simulations have been run using potentials in which these terms have been altered by inserting a multiplicative factor in eqn (2.9), and significant effects
on the two-dimensional melting temperature have been observed. In addition to increasing the periodicity as suggested by the theory, it is of interest to study adsorption on the planar surface produced when the periodicity is simply omitted. Evidently, such changes are straightforward to implement in a simulation, but impossible to produce in the real world. Another questionable aspect of the potential function used in the simulation arises from the choice of the bulk N-N site well-depth lgc Theoretical calculations [12,13] of the dispersion interaction of a pair of particles in the neighborhood of a semi-infinite dielectric solid indicates that the potential well should be roughly 80% of its bulk depth, a finding which is supported by a recent analysis [14] of experimental data for N2 on GCB. To date, no simulations have been performed using a modified N,-N, potential, with the exception of changes in the electrostatic terms made in response to the suggestion [8] that one should use a quadrupole moment which is smaller than the known experimental number [ 151by roughly 20%. Typically, simulations were run for 6~-12,~ time steps after equilibration. Coordinates and velocities were saved at intervals of -20 time steps to form a data set which could then be used to evaluate various averages of interest. Not all properties can be successfully studied over real-time intervals of lo-20 psec-for example, the growth of long-range translational order which occurs when a disordered phase freezes is too slow for observation. In general, orientational changes, including the growth of an ordered phase from a disordered one, were easily simulated. Translational motions involving interchange of molecules and similar large-scale displacements usually occurred only in the neighborhood of defects (i.e. vacancies); the formation of such defects could be a severely rate-limiting step in dense layers at relatively low temperatures. Figure 3 shows a phase diagram for nitrogen adsorbed on GCB at low temperature, with the state points studied by molecular dynamics computer simulation indicated by the filled circles. We will not attempt to discuss all the results obtained, but will here concentrate on the coverage-dependent series obtained at 73.6 IL 3. RESULTS
We begin with a brief consideration of the in-phase orientations of the coplanar N2 molecules in a commensurate monolayer. The herringbone ordering shown in Fig. 1 is lost at a low temperature of 27 K (experimental) or 33 K (from simulations which utilize the experimental Nz quadrupole moment [15]). Incomplete (and unpublished) simulations indicate that a decrease of 20% in Q causes a sharp drop in this simulated CO -+ CD transition temperature to a value near 23 IS. At first, this sensitivity is surprising, especially since the electrostatic terms are a minor part of the total NI_N, energy. In order to better
Adsorption of N, on the graphite basal plane
_
_
..---
I ,,’
-1195 _
_
_____-_-*
e_
,
.y
I’
’ ,’ ,’
SC’ _,200 ~-.-__.__.--~C .a0
.20
.40
.6E
.8E
I .E0
._.._..* -1266 ~, 40
,,I .20
, .e
,-, .6E
80
, 1.00
DISTANCE/a
Fig. 2. Variation in the minimum N,lgcb potential energy (divided by Boltzmann constant k) are shown for translations of a N2 molecule in the directions indicated. Energies were minimized with respect to orientation and distance from the surface. The various curves are for eq. (2.6) (A) and for modified versions of this expression in which the attractive term is multiplied by 1 - y,(3/2 cos’p - 1) and the repulsive term, by 1 - ~~(615cos2p-1), where p is the angle between the molecular axis and the surface normal. Curve B corresponds to ?A = 0.4, ‘R = -0.54 and D, to ?A = 0.4, ?R = -0.54.
understand the problem, it is helpful to calculate the dependence of this energy upon the in-plane angle 4. Several possibilities can be explored: for example, energy can be calculated as one molecule is rotated with the surrounding molecules held fixed in some ~n~guration [12], or all molecules can be rotated simultaneously. Figure 4 shows the lateral interaction energy of a molecule in a coplanar commensurate herringbone lattice plotted as a function of the angle between the molecular axes and the glide line (the two sublattices correspond to + and - angles of the same magnitude). In this case all molecules are rotated, but those in different sublattices rotate in opposite directions. The LJ site-site part of the interaction is much larger than the electrostatic, ranging from - 7 to - 10 reduced units relative to a range of + 1 to - 1 for the electrostatic However, it is important to realize that the temperature and energy for orientational disordering transition depend upon the change in the NrN2 energy with 4, rather than the absolute value. From the plots in Fig. 4, it is evident that the changes are of comparable magnitude for the two contributions. Qualitatively, a pronounced sensitivity of the disordering transition upon the magnitude of Q2, the constant of proportionality for the electrostatic energy, can now be understood. If the temperature of the commensurate solid N2 layer is raised to -49 K, the phase diagram of Fig. 1 shows that melting to a fluid occurs. Interestingly,
no two-dimensional liquid phase has yet been observed for Nz on gcb. It appears that the commensurate solid is sufficiently stabilized by the substrate to give a melting point higher than the expected (or “incipient”) gas-liquid critical temperature [ 161. One of the most interesting features of the Nz phase diagram is the “tongue” of commensurate solid that extends past 49 K up to -85 K. At these high temperatures, this is undoubtedly a defect solid that fills the surface-indeed, one that exerts a considerable two-dimensional spreading pressure. The precise nature of this solid phase is still a matter of speculation, but simulations performed over a range of coverage at 73.4 K do provide some info~ation about it, as well as less dense surface layers, and about the details of second layer formation in this system. We now discuss some of this work. The relevant equilibrium quantities obtained in the simulation studies include (a) distributions of in-plane angle 4 and the value of an associated order parameter OPl = (cos 2(+ +,,)), where & denotes the orientations in the minimum energy herringbone lattice. All simulations at 73.6 K indicate that in-plane order is quite negligible, i.e. that 4 angles are randomly distributed in these layers. (b) distributions of the out-of-plane angle p that denotes the angle between a molecular axis and the
W. A. STEELEet al.
Biloyer
fluid
I
0
0
I I I IO 20 3&O 50 60 70 80 90 Temperature Fig. 3. A temperature-coverage phase diagram for the NJ gcb system. Points subject to computer simulation are indicated by filled circles. The phase boundaries are based onexperimentaldata, with speculative boundariesindicated by dashed lines. Notation is: CO - commensurate, in-plane orientational ordering; CD - commensurate, in-plane disordered; UIO - uniaxial compressed, incommensurate, orientationally ordered; TIO - triangular, incommensurate, orientationally ordered; TID - triangular incommensurate orientationally disordered.
face such as that used in the simulation, one ordinarily ascribes this variation to a change in the average lateral interaction, with the mean number of neighbors to a molecule in the film increasing linearly with density. However, we can break up the total potential energy into its component parts in the simulation and show that this description is not quite correct. The assumption being made is that the average molecule-surface energy is independent of coverage, and Fig. 6 shows that this is in error. There are two contributions to the slow (linear) variation in this energy: the largest is due to the coveragedependence of the mean molecule-surface orientation angle, which carries with it a significant change in the average energy; and the second is that the molecules in the simulation tend to be noticeably more localized over the hexagonal site centers as the coverage approaches commensurate, which then gives rise to a small energy of localization that is opposite in sign to the energy associated with the orientational change. If one takes these factors into account, the linear variations give an average lateral interaction U,,, in a layer with the commensurate that is equal to 2.21 kJ/mole, compared to the value of 1.58 kJ/mole obtained if one ignores the change in molecule-surface energy. This number is very roughly half the value of 5.0 kJ/mole for the average potential energy of a molecule in bulk liquid N2 at 74 K, as one might expect from the decrease in the number of neighbors in going from three to two dimensions. We can pfogfess a bit further by using this result to estimate Upair, the average lateral energy per pair of N, molecules in the monolayer. This quantity can be written as
normal to the surface. In this case, a preferred orientation is found parallel to the surface; however, these distributions are broad (because of thermal agitation) and coverage dependent.
(3.2)
The computer-generated snapshot of a top view of the N2 layer shown in Fig. 5 reveals that a significant fraction of the molecules are standing more or less on end. It is not surprising to find that molecular crowding experienced near monolayer completion is partially relieved by an increase of this on-end fraction, relative to its value in a dilute film. This change has the consequence that the average molecule-solid energy changes noticeably as coverage increases. This point is illustrated in Fig. 6, where simulated values for the average potential energy per N2 molecule are plotted as a function of the surface coverage N:, measured in multiples of the commensurate layer coverage. The first point to notice is that the simulation yields 9.90 kJ/mole for the isosteric heat in the limit of zero coverage (equal to RT-potential energy (at 0 = 0)). This differs from the experimental result, corrected for surface heterogeneity, by less than 0.1 kJ/mole. As coverage increases, a linear variation in energy is obtained. For adsorption on a uniform sur-
where N/A is the surface density and ~*,a denote the center-center separation (in reduced units, defined here by dividing by a,,) and the molecular orientations, respectively. n is the angular integration weight = (47~)*,for two linear molecules. g(r,Q) is the correlation function that gives the excess (over random) probability of finding a pair of molecules in the monolayer with given ~,a, and U,, is the N,-N, interaction energy defined in eqn (2.50). Since N/A = 0.0636 molecules/A2 in the commensurate layer, we find apair = 3.1 kJ/mole from the simulation. Evidently, an evaluation of this quantity using eqn (3.1) directly will necessitate a careful study of the effects of restrictions in the orientational variables upon the monolayer pair correlation function. At present, neither theory or simulation is available to help us perform such a calculation. A steep drop in the total energy curve of Fig. 6 occurs at coverages greater than commensurate-
Adsorption of Nz on the graphite basal plane
I
+I
0
I 30
I
I
I I 60 90 Glideline Angle
13
I
I
I I20 (degl
I 150
-7
180
Fig. 4. Variation in lateral NTNl potential energy in a commensurate herringbone lattice. Glide-line orientation angles + for all molecules are varied simultaneously, and energy changes are calculated separately for the Lermard-Jones site-site part (dashed curve) and for the quadrupolar electrostatic part (sofid curve).
this is due primarily to the onset of adsorption in the second layer with a much diminished molecule-solid interaction. An interesting feature of the energy variation in this region of coverage is that the constant or gently rising curve that one might expect for simple growth in the second layer density is clearly absent in this simulation. The detailed molecular picture generated in this approach allows us to understand why this is so. It is simply that layer-by-layer growth
does not occur in this system at this temperature. This process undoubtedly occurs for other adsorbates, especially at relatively low temperature. However, the Nz is a supercritical fluid, at least in two dimensions, and the simulations show that in the region N;” = 1-2, coverages and molecular axis orientation are significantly changing in layer 1 and, to some extent, in layer 3 as well as in layer 2. Plots of the coverage-dependent distributions of molecular
Fig. 5. Snapshots of top views of the computer-generated instantaneous con~gurations of N2 molecules on gcb at 73.6 K. The right picture is for a layer at the commensurate density and the left one, for a partially complete monolayer. In-plane and out-of-plane orientations relative to the surface are indicated by dumbbell and double circle shapes, respectively. Carbon atoms are shown as dots.
W.
14 1
I
I
A.
STEELE
et al.
I
102 101 NI 100 I
Y
zr’ 0
5000
10000
Timesteps
Fig. 8. A typical plot of the number of molecules in layer 1 as a function of time at 73.6 K. A commensurate layer corresponds to 96 N2 molecules; total coverage amounted to 124. The length of one timestep was chosen such that 10,000 timesteps = 2.0 picosec. The equilibrium value of N, in this case is -102 molecules. Fig. 6. Computed curves for the average potential energy per molecule as a function of coverage, in fractions of the commensurate monolayer. The total energy is the sum of the NJgcb and the N2/N2 interactions. The sharp break at N,* = 1 is due primarily to the onset of second layer adsorption, as indicated in Figure 7.
orientation are shown elsewhere; Fig. 7 shows simulations of the changes in layer densities with total coverage. Another noteworthy point deals with the question of equilibration between layers. Recollect that a bilayer simulation is begun with an arbitrary partition of molecules between layers which is almost certainly incorrect. To reach equilibrium, it is necessary that a reasonable number of interlayer transfers occur during the simulation. As shown in the plots of the number of molecules in layer [l] as a function of time (see Fig. 8), the number of molecules per layer fluc-
tuates by one (or occasionally two) molecules per 100 on the time scale of roughly one per picosecond over a run of perhaps lo-psec duration. This appears to be the first determination of the monolayer lifetime as well as the magnitude of the fluctuation in layer density. Of course, the values obtained are not really quantitative, but this can be remedied either by running longer simulations or going to higher temperature or both.
4. DISCUSSION
In this brief review, we have attempted to show the power of the computer simulation technique as well as discussing the assumptions that underly it. At least some of the questionable aspects can easily be improved upon, at the expense of more computer
Fig. 7. Coverage in the i’th layer at 73.6 K as a function of total coverage, calculated as fractions of a commensurate monolayer. The break in the smooth variation for the monolayer (i = 1) is due to a preference for commensurate density. Note that adsorption into second and third layers occurs simultaneously at N,* > 2, rather than as a simple layer-by-layer filling.
Adsorption of N, on the graphite basal plane time. For example parameters of the pair interaction or the molecule-solid potential can be altered to more realistic values, if such changes are indicated by new evidence. Alternatively, one can study adsorbates that are polar or flexible; certain simple types of surface heterogeneity (such as steps in the surface) can be introduced; adsorption of mixtures is readily handled; and extension of this work to higher coverages can also be done. Time-dependent properties such as molecular diffusive behavior in these films have been obtained from the molecular dynamics simulations, but are not discussed here. One could even relax the assumption of a rigid solid by including the carbon atoms in the dynamics, but this would probably require an inordinate amount of computer time. Another time-consuming improvement would be to handle the many-body terms in the intermolecular interactions explicitly rather than in the handwaving, but economical, way of using “effective” pair-wise models. The attentive reader will note that one important topic has been avoided in this discussion of results. This is the simulation of the adsorption isotherm. The problem here is that the isotherm is essentially a measure of the density dependence of the adsorbate chemical potential, and this is notoriously difficult to evaluate by computer simulation in either the Monte Carlo or the molecular dynamics approach. Specifically, one has
where a,g denote adsorbed and gas phase, respectively, &is the standard-state chemical potential and p is the isotherm pressure at equilibrium with an adsorbed film containing n, molecules per unit area. There are several possibilit&s for evaluating u, by simulation: First, note that E,, the partial molal energy of the adsorbed film, is readily simulated but s,, the partial molal entropy, is not. In principle, the components of the pressure tensor perpendicular and parallel to the surface can be simulated. Although the perpendicular term piZ must be equal to the gas phase pressure, it is generally obtained as a very small difference between two very large terms arising from the attractive and the repulsive parts of the z-component of the force on a molecule in the film. Consequently, the statistical errors in this calculation are often unacceptably large. On the other hand, one could use the spreading pressure to calculate the isotherm from the Gibbs adsorption equation, which we write as
ln(pK,in,)
with
= r&
c&
- 1) 2,
(4.2)
15
and KH is the Henry’s law constant defined by lim (Pins) *.-,I
= KH.
(4.4)
The components of the pressure tensor P,,, P,, and P,, vary from point to point in the fluid near the surface, and thus an average over the entire volume V, of the film is needed in eqn (4.3). Although there is a well-known algorithm to evaluate these pressures in a simulation based on the virial theorem of Clausius, a number of complications arise when attempting to treat a molecular gas in the external field due to the molecule-solid potential, especially if this field has a significant periodic variation in the plane of the surface. (Work on this problem is currently in progress.) Alternatives to this route to the isotherm include various versions of the potential-distribution theory, which are essentially calculations of the work or energy of adding a particle to the fluid phase [ 181. We explicitly mention only one of these which has been employed by Knight and Monson [19] in an adsorption calculation. Although their work is based on the Monte Carlo algorithm, this distinction is unimportant here since both algorithms are capable of generating datasets containing the molecular coordinates for configurations that characterize the statistical ensemble for the adsorbed phase. The potential-distribution calculation relies on the generalized isotherm equation.
$$T)=-(%$ = -In
5 , ( %-Ii
(4.5)
(45)
where &a is the configurational integral for n, molecules near the surface. Knight and Monson then write
z nyZ “PI
h
(exp(-uWWno_,,
(45)
where u(r) is the interaction energy of a molecule at point r with the neighboring adsorbate and with the solid, and p(r) is the local density at that point. The average is taken over the configurations for n, - 1 particles. For this reason, the addedmolecule is often described as a ghost particle, since it does not affect the trajectories or configurations of the n, - 1 others. The computation thus consists in evaluation of the exponential factor in eqn (4.5) for a large number of ghost particles inserted in a large number of molecular configurations. Of course, the result will often be that the ghost experiences a large positive interaction, especially in a dense phase where the probability of overlap with a real molecule is quite large. If the exponential is always close to zero, it is difficult to get good statistics in this region of the isotherm. A final possibility for evaluation of the chemical potential is to recollect the fundamental theorem of
W. A.
16
STEELE
statistical mechanics which states that the fluctuation in the number of molecules N in a volume of fluid is related to chemical potential by 772
-
772
&=-
F
(4.5)
The quantity 6 can be evaluated in a simulation as long as the sample is large enough to eliminate the effects of the periodic boundary conditions which, for example, totally suppress the ~uctuations in a volume coincident with the computer box size, If one can obtain accurate values of S by statistical averaging using small volumes, one can integrate (+./ an,) with respect to n, from the Henry’s law limit up to monolayer or even multilayer coverage. Of course, this necessitates a series of simulations closely spaced
et
al.
over the relevant range of coverage as well as working with a one-phase system, since the fluctuation theorem is no longer useful at densities corresponding to two-dimensional phase coexistence. To date, simulations of adsorption isotherms are almost nonexistent, other than Knight and Monson and some early Monte Carlo studies [20] of rare gases on gcb using the grand canonical ensemble. Still, the range of molecular such properties that can be generated by computer simulation is impressive and enlightening, and it is highly likely that studies of other molecules as adsorbates and other solids as adsorbents will yield much more useful info~ation in the future. Acknowledgments-Partial support for this work was provided by the Division of Materials Research of the National Science Foundation,
REFERENCES
J. K. Kjems, L. Passell, H. Taub, J. G. Dash and A. D. Novaco, Whys.Rev. B 13,1446 (1976); J. Eckert, W. D. Ellenson, J. B. Hastings and L. Passe& P&r. Rev. Lett. 43,1329 (1979); R. D. Diehl, M. F. Toney and S. C. Fain, Jr., Phys. Rev. Letf. 48, 177 (1982); R. D. Diehl and S. C. Fain, Jr., Surf Sci. 125, 116 (1983); S. C. Fain, Jr., M. D. ChinnandR. D. Diehl, Phvs. Rev. B 21. 4170 119801: M. Nielsen. K. Kiaer and J. Bohr, J. &iectron’Spec: Relat. Phenom. 30,‘lll (1983); K. Morishige, C. Mowforth and R. K. Thomas, surf. .Sci. 151, 289 (1985); R. D. Diehl and S. Fain, Phys. Rev. B 26, 4785 (1982); R. D. Diehl and S. C. Fain, Jr., J. Chem. Phys. 77,5065 (1982); H. You andS. C. Fain, Jr.,D&x~s. FaradaySot. 81,159(1985). D. M. Butler. G. B. Huff. R. W. Toth and G. A. Stewart, Phys: Rev. Lett. 34, 1718 (1975); A. D. Migone, H. K. Kim, M. H. W. Chan, J. Talbot, D. J. Tildesley and W. A. Steele, Phyr. Rev. Left. 51, 192 (19831: 0. M. Zhana. H. K. Kim and M. H. W. Chan. Phys.‘Rev. B 33,41~~1986); Q. M. Zhang, H. K. Kim and M. H. W. Chart, Phys. Rev. R 32,182O (1985); M. H. W. Ghan, A. D. Migone, K. D. Miner and Z. R. Li, Phys. Rev. B 30,268l (1984); J. Piper, J. A. Morrison, C. Peters and Y. Ozaki, J. Chem. Sot. Faraday Trans. 179,2863 (1983); T. T. Chung and J. G. Dash, Surf. Sci. 66,559 (1977). J. Talbot, D. J. Tildeslev and W. A. Steele, Mol. Phys. 51, 1331 (1984); J. Taldot, D. J. Tildesley and W. A. Steele. Surface Sci.. 169.71 (1986). A. V. Vernov and W. A.‘Steele, Langmuir 2, 219 (i986); A. V. Vernov and W. A. Steele, Surf. Sci., 171,83 (1986); C. Peters and M. Klein, Mol. Phys. 54, 985 (1985); Y. P. Joslin and D. J. Tildesiey, Mol. Phys. 55,999 (1985); J. Talhot, D. J. Tildeslev and W. A. Steele, Discuss. Fara~v Sot. 81,91(1985j. See, for example, W. B. Streett and D. 3. Tildesley, in Computer Mode~l~ngof Matter (Edited by P. Lykos). A.C.S. Symposium Series, Vol. 86; S. Murad and K.
5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16.
17.
18. 19. 20.
Gubbins, in Computer ModelZingof Matter (Edited by P. Lvkos). A.C.S. Svmuosium Series. Vol. 86: D. J. Evans, kol. Phys. $, 317 (1977). Washington,‘D.C. C. W. Gear, Numerical Init~aIValue Problems in Ord~nary Di~erent~u~Equation. Prentice-Hall, Englewood Cliffs, NJ (1971). D. J. Evans and G. P. Morriss, Comput. Phys. Rep. 1, 297 (1984); Chem. Phys. 87,451 (1984). _ S. Nose. J. Chem. Phvs. 81. 511 (19841. C. S. ikurthv, K. Smger, M. L. Klein and I. R. MacDonald, &iol. Phys,41,1387 (1980). W. A. Steele. J. Phvs. (Paris) C4. 61 (1977). W. A. Steele; .I. P/&s. &em: 82,817 (1978). G. Vidali and M. Cole, Phys. Rev. B 29,6736 (1984). L. W. Bruch, J. Chem. Phys. 79,3148 (1983). 0. Sinanoelu and K. S. Pitzer. 1. Chem. Phvs. 32.1279 (1960); A.“D. MacLachlan, &or. Phys. 7,381 (Ib64). M. J. Bojan, Ph.D. thesis, Penn State University, 1986. A. D. Buckingham, R. L. Disch and D. A. Dunmur, J. Am. Chem. Sot. !Mt,3104 (1968); N. J. Bridge and A. D. Buckingham, J. Chem. Phys. 40,2733 (1964). A. D. Migon< M. H. W. Ghan,.K. J. Niskanen and R. B. Griffiths, J. Phys. (Paris) C 16, L1115 (1983); D. M. Butler, J. A. Litzinger, G. A. Stewart and R. B. Griffiths, Phys. Rev. Len. 42, 1289 (1979). W. A. Steele. The Interaction of Gases with Solid Surfaces, Sec. ‘4.8. Pergamon, New York (1974), sec. 4.8. B. Widom, J. Chem. Phys. 86,869 (1982); K. Mon and R. B. Griffiths, Phys. Rev. A 31,956 (1985). J. F. Knight and P. A. Monson, J. Chem. Phys. 84, 1909 (1986). D. Nicholson, N. G. Parsonage and L. A. Rowley, Mol. Phys. 44,629 (1981); L. A. Rowley, D. Nicholson and N. G. Parsonage, Gol. Phys. 31, k5,389 (1976); J. Commit. Phvs. 26.66 (1978): J. E. Lane and T. H. Spurling, Aus; J. khem. 291.2103 (1976); 31, 933 (1978).
DISCUSSION Question by P. Pendleton, American Cyanamid In these systems, one would expect to observe an induction effect through the graphite substrate as introduced by McLachan. Firstly, for sub-picosec intervals, will future models include these sub-surface
Adsorption of N, on the graphite basal plane atom-mirror image interactions? Secondly, if these terms are included, introduced in a model with step heterogeneities?
how much complication
17 would be
Reply by Steele There is little doubt that McLachan-type substrate-mediated interactions are significant in these systems and indeed, the problem has been directly studied by Bruch (ref. 12). One attempts uses “effective” twobody functions which crudely accounts for these alterations for the simple reason that the full many-body potential would run so slowly in the available computers as to eliminate all meaningful results. We are among those who would very much like to be able to do a better job in this respect. As for the second question, the answer is unknown, but would make a worthwhile thesis project. Question by K. S. W. Sing, Brunel University Could you please comment on the factors which are likely to influence the form of the multilayer isotherm? In our experience the form of the nitrogen multilayer at 77 K (as analyzed by the FHH method) is almost independent of the nature of the substrate surface. Reply by Steele This is a difficult question to answer, especially since simulations for a single substrate do not necessarily help in understanding the differences or similarities with other substrates. The only suggestive feature that comes to mind is that both these and the older Ar/gcb multilayer simulations indicate that several layers are filling simultaneously at total coverages greater than 8 = 1.5 and at temperatures equal to or greater than the bulk boiling point. This smears out the energetics somewhat and might be construed as indicating an entropy-dominated process. Much more work needs to be done to fully answer the question. Question by H. A. Resing, Naval Research Lab Is the failure to observe orientational or translational computer time?
melting of the 2-D solid perhaps due to lack of allowed
Reply by Steele Orientational degrees of freedom come to equilibrium quite rapidly in these systems and no problems of equilibration time are encountered at any temperature of interest. The situation is more complex when translational melting (or motion) is studied. A strong dependence upon temperature and surface density is observed which appears to reflect the fact that molecular interchanges or site-to-site diffusion is more difficult in 2-D than in 3-D. This may indicate that local density fluctuations are smaller because the degree of freedom perpendicular to the surface is strongly hindered. We believe that translational equilibration times have become marginal on the time-scale of the computation for T > 65, and that they become satisfactorily short as temperature increases by another 15”. Question by A. Groszek, British Petroleum Co. What are the densities of adsorbed NZ before and after the melting transition? Reply by Steele The solid layer remains at the commensurate density at temperatures up to the melting point of -47 K for a patch of solid. The patch actually melts to a 2-D vapor with density determined by (number of molecules/ total area available). If the surface is completely filled with solid, melting appears to occur at constant density but at temperatures that are very dependent upon the initial density of the solid. Question by V. Deitz, Naval Research Lab Does the calculation include interactions with basal planes below the surface basal plane? How many basal planes contribute to the energy? Reply by Steele The interactions are summed over all basal planes. The summation converges rapidly even for the relatively long-range attractive terms. In particular, in the region of the minimum energy, the entire contribution of all planes excluding the surface is less than 15% of the total.
Carbon 2.5, No. 1, pp. 1-17, 1987 Printed in Great Britain.
0008-6223187 $3.00 + .Ml 0 1987 Pergamon Journals Ltd.
STUDIES OF THE ADSORPTION OF N2 ON THE GRAPHITE BASAL PLANE BY COMPUTER SIMULATION Department
WILLIAM A. STEELE of Chemistry, The Pennsylvania State University, University Park, PA 16802, U.S.A.
ALEXEI V. VERNOV Department of Physical and Colloid Chemistry, People’s Friendship University, 117302 Moscow, U.S.S.R. and
Department
DOMINIC J. TILDESLEY of Chemistry, The University, Southampton,
England
(Received 1 May 1986) Abstract-A molecular model for the physical adsorption of nitrogen on graphitized carbon black (gcb) has been utilized extensively in computer simulations of this system. After a brief review of the details of the model and of the simulation technique, a number of results produced in the course of these studies are summarized. These include the thermodynamic and structural properties of N, adsorbed on gcb at 73.6 K and at coverages ranging up to -2 l/2 monolayers. Changes in molecular orientation relative to the surface are found to have a significant effect on the lateral interaction in a layer; variations in packing density in the first layer are discussed; and the rate of molecular interchange between first and second layers is characterized. Likely directions for future simulation studies of physisorbed films are proposed. Key Words--Adsorption,
nitrogen, basal plane, computer simulation.
1. INTRODUCTION
For many decades, the characterization of surface area, pore structure and related properties of solid adsorbates has begun with measurement of nitrogen isotherms at temperatures near the normal boiling point of this inert nonpolar adsorbate, often followed by a BET surface area determination or, in the case of porous materials, a calculation of pore volume from either isotherm hysteresis or Dubinin-Radushkevich analysis of the data or both. In spite of the success of these procedures, there is much to learn about the molecular level properties of physisorbed nitrogen, even on homogeneous nonporous solids. Recently, a number of sophisticated experimental techniques have been applied to this problem. These include diffraction of electrons, neutrons and X-rays [l] to characterize the various solid-like ordered layers that form at low temperatures when nitrogen is adsorbed on the exposed basal planes of graphitic carbon. In addition, phase equilibria and the thermodynamic changes associated with two-dimensional (or quasi-two-dimensional) transitions have been obtained from high-precision isotherms and heat capacities [2]. These studies have revealed a rich phase diagram for the N,/gcb (graphitized carbon black) system with much of the detail appearing at T < 50
K. Nevertheless, some interesting and unexpected behavior has been observed at T = 78 K-for example, the presence of an ordered, solid-like film at coverages close to monolayer completion has been observed by neutron diffraction and confirmed by Xray diffraction. An alternative view of the molecular behavior in physisorbed films can be obtained by applying the modern algorithms of computer simulation to these systems [3]. This paper is devoted to a review of some of the information gained in this way. We will concentrate upon the coverage dependence of the properties of physisorbed nitrogen at 73.6 K, particularly at coverages in the dense monolayer to bilayer region. However, before discussing the results of these studies, the technique will be briefly together with a description of the assumptions needed to give a realistic description of the system.
2. SIMULATION METHODOLOGY
The technique used is a well known and wrdely applied method in which it is assumed that molecules obey the laws of classical Newtonian dynamics so that, for example, the center-of-mass position r, and velocity v, of particle i with mass m at a time I t Al
8
W.
A. STEELEet al.
can be obtained from r,(t + At) - r,(t) = v(t) At + correction terms,
(2.1)
v,(t + At) - v,(t) = ; F,(t) At + correction
terms,
(2.2)
where F, is the force on molecule i due to its interactions with other molecules and, in the case of physisorption, with the solid. Molecular trajectories for the system of interest are generated by assuming a realistic intermolecular potential (and molecule-solid interaction) and calculating the force by differentiating these functions. In the case of nonspherical molecules, equations for the angular coordinates and velocity w, must also be solved. For example, w,(t + At) - o,(t) = f N,(t) At + correction
terms,
(2.3)
where I and N, are the moment of inertia and torque, respectively. The calculation of the change in molecular orientation requires the introduction of a convenient set of coordinates. We here use quaternions-their definition and equations of motion are discussed in detail elsewhere [4]. In the computer, the slow step in generating the trajectories is the sum over all pairs of molecules (or sites within the molecules, for the model potentials discussed below) to obtain the force and torque on each molecule in the system. Thus in order to study the system properties over as long a time as possible, one wishes to maximize the timestep size At, consistent with an accurate set of trajectories. This is done by using one of several suggested prescriptions for the correction terms in eqns (2.1)-(2.3)-the algorithm used in the studies discussed here was a predictor-corrector method due to Gear [5]. The simulation is run for a sample of 100-250 molecules which are initially placed in some convenient set of positions, such as an ordered layer next to the surface, and with a random distribution of velocities chosen such that m(v’) = 3kT, I(d)
= 2kT,
(2.4)
where T is the desired temperature and the brackets denote an average over all molecules (and, when appropriate, over time steps as well). In the initial stages of the simulation, the velocities are resealed to keep T constant. For a timestep that is taken to be -1Om’5 sec. in the systems studied here, the systems appear to reach statistical equilibrium after 30006000 time steps. At this point, velocity resealing is usually terminated and the molecular system is allowed to evolve in time as a microcanonical ensemble
(constant particle number, volume and total energy). It is now well known that statistical averages obtained from such simulations conform nicely to the postulates of statistical mechanics, yielding Maxwellian velocity distributions and the expected values for the various fluctuation theorems. One other crucial assumption which must be made for this procedure to be successful arises in connection with the boundary conditions for a sample of (say) 100 molecules placed in a layer near the surface of a solid. We wish this sample to mimic the real adsorbed film containing -lOI molecules while maintaining the “adsorption volume” constant. This is conventionally achieved by surrounding the original sample with images of itself and then allowing the “walls” to be permeable. Thus a trajectory that causes a molecule to pass out through a wall will cause one of its images to pass in through the opposite wall, thus conserving both number and volume. In the adsorption case, these periodic boundary conditions are needed only in the two dimensions parallel to the surface. Other than requiring the box size match the periodicity of the surface so that molecules crossing a boundary suffers no discontinuity in molecule-solid interaction, the implementation is straightforward. As a result, each molecule in the computer box will interact with a mixture of molecules and, if it is near a boundary, images. Evidently, errors can be produced if a molecule interacts with its own image-a problem for small boxes or for long-range interaction laws. Of more importance is the fact that the longrange translational order that develops in solidified films can be materially affected by the periodicity of the images. Thus detailed studies of the thermodynamics of the transitions involving such phases require great care. The boundary condition for motion perpendicular to the solid surface is handled by simply reversing the z-component of those molecules at some large distance from the surface (-10 atomic diameters) if they have been desorbed because of an energetic fluctuation. This seldom happens under the conditions of the simulations performed to date, but a related occurrence which can cause significant problems is the transfer of a molecule from one adsorbed layer to another, particularly when the transfer involves first layer molecules. Since the molecule-solid potential is quite large for first layer molecules, a large change in potential energy will result. Conservation of total energy then requires a significant change in the total kinetic energy. For the restricted number of molecules under study, the consequence is an unwanted large jump in temperature. In the simulations at 73 K discussed here, this problem was avoided by utilizing one [6] of several recent suggestions for a constant temperature algorithm-in essence, velocity scaling is continued for the entire time span of the simulation according to a prescription which approximates [7] contact with a heat reservoir. In specifying the molecular interactions for the NJ GCB system, one makes maximum use of the em-
Adsorption of N2on the graphite basal plane pirical results obtained in simulations of bulk N, and in zero-coverage extrapolations of experimental heat and isotherm data. In particular, it is known that a reasonable compromise between computational simplicity and accuracy for the bulk properties is given by the site-site model for N,-N, interactions [8] in which it is assumed that each molecule contains two sites centered on the atoms and separated by the bond length L = 1.10 A. These interact with the sites on other molecules by Lennard-Jones (LJ) energies plus a Coulombic term that accounts for the electrostatic interaction associated with the N, quadrupole moments. The electrostatic calculation is completed by placing two charges at the center of the molecule which are of opposite sign to the “atomic” charge 4. The molecule quadrupole moment Q is qL?/2, and q can be adjusted to give the desired “effective” Q for a molecule in a dense adsorbed phase. The total N2-NZ interaction U,,(ij) for molecules i and j can now be written as
‘)
where we have switched to reduce units denoted by asterisks, defined for energy by dividing by Ed\, and for distances, by dividing by a = 2.46 A, the length of the lattice-vector for the graphite basal plane. The periodic term in the energy is given to good accuracy
by 2rrAh u;, (7. 2) = * 0, _
where A = a,la, u,* = area of a unit cell of the graphite lattice = J312: g* is the reduced length of a reciprocal lattice vector = 41~/,/3. K, and Kc are modified Bessel functions of the second kind and j”(7) is a periodic function of position T = x,y relative to the surface lattice which is defined to be f(7) = -2[cos
2n(P
+ y*/,:3)
+ cos 27F(X* - y*/\,‘3) f cos 4ay”/,
It is assumed that appropriate units for 4, are used; r,, is the separation between a site in molecule i and one in j, and enn,u,, are the empirical well-depth and size parameters for N-N. Values used in this work are egRlk = 36.4 K and uXg = 3.32 A. In the spirit of the site-site approach, one can devise [9] a molecular version of the atom-atom summation approximation to model N, molecule-carbon black potential. This is known to give a reasonably accurate account of rare gas/GCB interactions [lo], and its extension to the molecular U&i) can be written as
where the sites involved in r,, are nitrogen “atoms” in molecule i interacting with carbon “atoms” in the solid. Well-depth and size parameters that give good agreement with the low coverage experimental data for this system are EJk = 31.9 K and up, = 3.36 A. Equation (2.6) is not convenient for use in computer simulation since it involves a sum over a large number of carbon atoms and two N sites for each molecule. Thus one tranforms [lo] it to
Ah
+ 12.42(z* + 0.994)” _ 2z*’ + 9.66z* + 13.33 8.28 (z* + 1.38)’ +
u;,
(7. z),
(2.7)
3).
(2.9)
The T-dependent terms shown in eqns (2.8) and (2.9) are only the first in an infinite series which correctly represents any property possessing the periodicity of the graphite basal plane. This plane is made up of hexagonal adsorption sites defined by the carbon atoms; the hexagons are spaced by 1 unit along the x direction and by 4312 units in the y direction. Only l/3 of these sites are occupied by NZmolecules in the commensurate film, which thus is described as a (,‘3 x J3) R30” overlayer. Nearest neighbor hexagons cannot both be occupied in such films because of the large size of the Nz (>3.3 A) compared to the site spacing (2.46 A). It is interesting to note that the area per molecule is 15.7 AZ in the commensurate film, a value close to the 16.2 A: taken for the N, area in the BET calculation. It is evident that these model potential functions embody several desirable features. such as the nonsphericity of the NZ molecules. Consequently. the minimum energy configuration for the commensurate layer is the orientationally ordered arrangement shown in Fig. 1. One can readily imagine that this order, both in-plane and out-of-plane, will be lost as temperature is increased; furthermore, changes in surface density will also produce alterations in the angular ordering for the minimum energy structure as well as affecting the temperature dependence of the disordering process. These expectations have been confirmed by simulation, as will be discussed below. Another important feature of the molecule-solid potential is that it includes a reasonable estimate of the variation in energy that occurs as a molecule translates parallel to the surface. Figure 2 shows this variation for a molecule placed parallel to the surface at a distance corresponding to the minimum energy over the center of a hexagon. This variation can give rise to “localized” or “site-wise” adsorption, but no such restrictive approximation need be used here.
W. A. STEELEet at.
**
**
‘W. *.
**
*/ .
\
*.
.
.*
**
.*
**
**
**
**
/.*.**~..~**/ ,*.
.
.
.
.
hi’. :_** .‘,‘* .*,*. .*;. .‘,‘. .
B
* * . .+..
**
.*-**
. .
. .
.*
.*
.a
.
.
*:I:. v-- -_1 * .;,;. .-_+ +-_. .I{. .y_. *I<
Fig. 1. A top view of the minimum knergy configuration for model nitrogen molecules on the basal plane of graphite. Dots indicate carbon atoms. The N2 molecular axes are all parallel to the surface plane, with in-plane orientations indicated by the solid lines. The herringbone pattern comprises two sublattices with a glideline indicated by AB. The orientation angles 4 relative to the glideline are 245”.
The calculated energies of l~alization are rather small, and would hardly be expected to produce sitewise adsorption without significant enhancement from the N2-N2 potentials. Indeed, the successful representation of the site-site summation for the molecule-solid potential by eqns (2.4)-(2.6), which contain a very limited number of periodic terms, is only possible when the periodicity is small. The question of whether the magnitude of the periodicity obtained from this model is correct remains debatable. Physically, it is small because the carbon atoms are densely packed in the graphite basal plane so that a N site interacts simultaneously with a relatively large number of C sites, regardless of its position over the lattice. Although this aspect of the calculation is unIikely to be greatly altered in alternative models, the questionable feature actually arises from the assumption of a symmetric N-C site interaction. In view of the great assymmetry in the carbon-carbon bonding in graphite, it has been suggested [12] that the C sites are not strictly spherical in their polarizabilities and interactions, and as a result, the periodic terms may be as much as twice the size as those obtained from spherical C sites. Consequently, some simulations have been run using potentials in which these terms have been altered by inserting a multiplicative factor in eqn (2.9), and significant effects
on the two-dimensional melting temperature have been observed. In addition to increasing the periodicity as suggested by the theory, it is of interest to study adsorption on the planar surface produced when the periodicity is simply omitted. Evidently, such changes are straightforward to implement in a simulation, but impossible to produce in the real world. Another questionable aspect of the potential function used in the simulation arises from the choice of the bulk N-N site well-depth lgc Theoretical calculations [12,13] of the dispersion interaction of a pair of particles in the neighborhood of a semi-infinite dielectric solid indicates that the potential well should be roughly 80% of its bulk depth, a finding which is supported by a recent analysis [14] of experimental data for N2 on GCB. To date, no simulations have been performed using a modified N,-N, potential, with the exception of changes in the electrostatic terms made in response to the suggestion [8] that one should use a quadrupole moment which is smaller than the known experimental number [ 151by roughly 20%. Typically, simulations were run for 6~-12,~ time steps after equilibration. Coordinates and velocities were saved at intervals of -20 time steps to form a data set which could then be used to evaluate various averages of interest. Not all properties can be successfully studied over real-time intervals of lo-20 psec-for example, the growth of long-range translational order which occurs when a disordered phase freezes is too slow for observation. In general, orientational changes, including the growth of an ordered phase from a disordered one, were easily simulated. Translational motions involving interchange of molecules and similar large-scale displacements usually occurred only in the neighborhood of defects (i.e. vacancies); the formation of such defects could be a severely rate-limiting step in dense layers at relatively low temperatures. Figure 3 shows a phase diagram for nitrogen adsorbed on GCB at low temperature, with the state points studied by molecular dynamics computer simulation indicated by the filled circles. We will not attempt to discuss all the results obtained, but will here concentrate on the coverage-dependent series obtained at 73.6 IL 3. RESULTS
We begin with a brief consideration of the in-phase orientations of the coplanar N2 molecules in a commensurate monolayer. The herringbone ordering shown in Fig. 1 is lost at a low temperature of 27 K (experimental) or 33 K (from simulations which utilize the experimental Nz quadrupole moment [15]). Incomplete (and unpublished) simulations indicate that a decrease of 20% in Q causes a sharp drop in this simulated CO -+ CD transition temperature to a value near 23 IS. At first, this sensitivity is surprising, especially since the electrostatic terms are a minor part of the total NI_N, energy. In order to better
Adsorption of N, on the graphite basal plane
_
_
..---
I ,,’
-1195 _
_
_____-_-*
e_
,
.y
I’
’ ,’ ,’
SC’ _,200 ~-.-__.__.--~C .a0
.20
.40
.6E
.8E
I .E0
._.._..* -1266 ~, 40
,,I .20
, .e
,-, .6E
80
, 1.00
DISTANCE/a
Fig. 2. Variation in the minimum N,lgcb potential energy (divided by Boltzmann constant k) are shown for translations of a N2 molecule in the directions indicated. Energies were minimized with respect to orientation and distance from the surface. The various curves are for eq. (2.6) (A) and for modified versions of this expression in which the attractive term is multiplied by 1 - y,(3/2 cos’p - 1) and the repulsive term, by 1 - ~~(615cos2p-1), where p is the angle between the molecular axis and the surface normal. Curve B corresponds to ?A = 0.4, ‘R = -0.54 and D, to ?A = 0.4, ?R = -0.54.
understand the problem, it is helpful to calculate the dependence of this energy upon the in-plane angle 4. Several possibilities can be explored: for example, energy can be calculated as one molecule is rotated with the surrounding molecules held fixed in some ~n~guration [12], or all molecules can be rotated simultaneously. Figure 4 shows the lateral interaction energy of a molecule in a coplanar commensurate herringbone lattice plotted as a function of the angle between the molecular axes and the glide line (the two sublattices correspond to + and - angles of the same magnitude). In this case all molecules are rotated, but those in different sublattices rotate in opposite directions. The LJ site-site part of the interaction is much larger than the electrostatic, ranging from - 7 to - 10 reduced units relative to a range of + 1 to - 1 for the electrostatic However, it is important to realize that the temperature and energy for orientational disordering transition depend upon the change in the NrN2 energy with 4, rather than the absolute value. From the plots in Fig. 4, it is evident that the changes are of comparable magnitude for the two contributions. Qualitatively, a pronounced sensitivity of the disordering transition upon the magnitude of Q2, the constant of proportionality for the electrostatic energy, can now be understood. If the temperature of the commensurate solid N2 layer is raised to -49 K, the phase diagram of Fig. 1 shows that melting to a fluid occurs. Interestingly,
no two-dimensional liquid phase has yet been observed for Nz on gcb. It appears that the commensurate solid is sufficiently stabilized by the substrate to give a melting point higher than the expected (or “incipient”) gas-liquid critical temperature [ 161. One of the most interesting features of the Nz phase diagram is the “tongue” of commensurate solid that extends past 49 K up to -85 K. At these high temperatures, this is undoubtedly a defect solid that fills the surface-indeed, one that exerts a considerable two-dimensional spreading pressure. The precise nature of this solid phase is still a matter of speculation, but simulations performed over a range of coverage at 73.4 K do provide some info~ation about it, as well as less dense surface layers, and about the details of second layer formation in this system. We now discuss some of this work. The relevant equilibrium quantities obtained in the simulation studies include (a) distributions of in-plane angle 4 and the value of an associated order parameter OPl = (cos 2(+ +,,)), where & denotes the orientations in the minimum energy herringbone lattice. All simulations at 73.6 K indicate that in-plane order is quite negligible, i.e. that 4 angles are randomly distributed in these layers. (b) distributions of the out-of-plane angle p that denotes the angle between a molecular axis and the
W. A. STEELEet al.
Biloyer
fluid
I
0
0
I I I IO 20 3&O 50 60 70 80 90 Temperature Fig. 3. A temperature-coverage phase diagram for the NJ gcb system. Points subject to computer simulation are indicated by filled circles. The phase boundaries are based onexperimentaldata, with speculative boundariesindicated by dashed lines. Notation is: CO - commensurate, in-plane orientational ordering; CD - commensurate, in-plane disordered; UIO - uniaxial compressed, incommensurate, orientationally ordered; TIO - triangular, incommensurate, orientationally ordered; TID - triangular incommensurate orientationally disordered.
face such as that used in the simulation, one ordinarily ascribes this variation to a change in the average lateral interaction, with the mean number of neighbors to a molecule in the film increasing linearly with density. However, we can break up the total potential energy into its component parts in the simulation and show that this description is not quite correct. The assumption being made is that the average molecule-surface energy is independent of coverage, and Fig. 6 shows that this is in error. There are two contributions to the slow (linear) variation in this energy: the largest is due to the coveragedependence of the mean molecule-surface orientation angle, which carries with it a significant change in the average energy; and the second is that the molecules in the simulation tend to be noticeably more localized over the hexagonal site centers as the coverage approaches commensurate, which then gives rise to a small energy of localization that is opposite in sign to the energy associated with the orientational change. If one takes these factors into account, the linear variations give an average lateral interaction U,,, in a layer with the commensurate that is equal to 2.21 kJ/mole, compared to the value of 1.58 kJ/mole obtained if one ignores the change in molecule-surface energy. This number is very roughly half the value of 5.0 kJ/mole for the average potential energy of a molecule in bulk liquid N2 at 74 K, as one might expect from the decrease in the number of neighbors in going from three to two dimensions. We can pfogfess a bit further by using this result to estimate Upair, the average lateral energy per pair of N, molecules in the monolayer. This quantity can be written as
normal to the surface. In this case, a preferred orientation is found parallel to the surface; however, these distributions are broad (because of thermal agitation) and coverage dependent.
(3.2)
The computer-generated snapshot of a top view of the N2 layer shown in Fig. 5 reveals that a significant fraction of the molecules are standing more or less on end. It is not surprising to find that molecular crowding experienced near monolayer completion is partially relieved by an increase of this on-end fraction, relative to its value in a dilute film. This change has the consequence that the average molecule-solid energy changes noticeably as coverage increases. This point is illustrated in Fig. 6, where simulated values for the average potential energy per N2 molecule are plotted as a function of the surface coverage N:, measured in multiples of the commensurate layer coverage. The first point to notice is that the simulation yields 9.90 kJ/mole for the isosteric heat in the limit of zero coverage (equal to RT-potential energy (at 0 = 0)). This differs from the experimental result, corrected for surface heterogeneity, by less than 0.1 kJ/mole. As coverage increases, a linear variation in energy is obtained. For adsorption on a uniform sur-
where N/A is the surface density and ~*,a denote the center-center separation (in reduced units, defined here by dividing by a,,) and the molecular orientations, respectively. n is the angular integration weight = (47~)*,for two linear molecules. g(r,Q) is the correlation function that gives the excess (over random) probability of finding a pair of molecules in the monolayer with given ~,a, and U,, is the N,-N, interaction energy defined in eqn (2.50). Since N/A = 0.0636 molecules/A2 in the commensurate layer, we find apair = 3.1 kJ/mole from the simulation. Evidently, an evaluation of this quantity using eqn (3.1) directly will necessitate a careful study of the effects of restrictions in the orientational variables upon the monolayer pair correlation function. At present, neither theory or simulation is available to help us perform such a calculation. A steep drop in the total energy curve of Fig. 6 occurs at coverages greater than commensurate-
Adsorption of Nz on the graphite basal plane
I
+I
0
I 30
I
I
I I 60 90 Glideline Angle
13
I
I
I I20 (degl
I 150
-7
180
Fig. 4. Variation in lateral NTNl potential energy in a commensurate herringbone lattice. Glide-line orientation angles + for all molecules are varied simultaneously, and energy changes are calculated separately for the Lermard-Jones site-site part (dashed curve) and for the quadrupolar electrostatic part (sofid curve).
this is due primarily to the onset of adsorption in the second layer with a much diminished molecule-solid interaction. An interesting feature of the energy variation in this region of coverage is that the constant or gently rising curve that one might expect for simple growth in the second layer density is clearly absent in this simulation. The detailed molecular picture generated in this approach allows us to understand why this is so. It is simply that layer-by-layer growth
does not occur in this system at this temperature. This process undoubtedly occurs for other adsorbates, especially at relatively low temperature. However, the Nz is a supercritical fluid, at least in two dimensions, and the simulations show that in the region N;” = 1-2, coverages and molecular axis orientation are significantly changing in layer 1 and, to some extent, in layer 3 as well as in layer 2. Plots of the coverage-dependent distributions of molecular
Fig. 5. Snapshots of top views of the computer-generated instantaneous con~gurations of N2 molecules on gcb at 73.6 K. The right picture is for a layer at the commensurate density and the left one, for a partially complete monolayer. In-plane and out-of-plane orientations relative to the surface are indicated by dumbbell and double circle shapes, respectively. Carbon atoms are shown as dots.
W.
14 1
I
I
A.
STEELE
et al.
I
102 101 NI 100 I
Y
zr’ 0
5000
10000
Timesteps
Fig. 8. A typical plot of the number of molecules in layer 1 as a function of time at 73.6 K. A commensurate layer corresponds to 96 N2 molecules; total coverage amounted to 124. The length of one timestep was chosen such that 10,000 timesteps = 2.0 picosec. The equilibrium value of N, in this case is -102 molecules. Fig. 6. Computed curves for the average potential energy per molecule as a function of coverage, in fractions of the commensurate monolayer. The total energy is the sum of the NJgcb and the N2/N2 interactions. The sharp break at N,* = 1 is due primarily to the onset of second layer adsorption, as indicated in Figure 7.
orientation are shown elsewhere; Fig. 7 shows simulations of the changes in layer densities with total coverage. Another noteworthy point deals with the question of equilibration between layers. Recollect that a bilayer simulation is begun with an arbitrary partition of molecules between layers which is almost certainly incorrect. To reach equilibrium, it is necessary that a reasonable number of interlayer transfers occur during the simulation. As shown in the plots of the number of molecules in layer [l] as a function of time (see Fig. 8), the number of molecules per layer fluc-
tuates by one (or occasionally two) molecules per 100 on the time scale of roughly one per picosecond over a run of perhaps lo-psec duration. This appears to be the first determination of the monolayer lifetime as well as the magnitude of the fluctuation in layer density. Of course, the values obtained are not really quantitative, but this can be remedied either by running longer simulations or going to higher temperature or both.
4. DISCUSSION
In this brief review, we have attempted to show the power of the computer simulation technique as well as discussing the assumptions that underly it. At least some of the questionable aspects can easily be improved upon, at the expense of more computer
Fig. 7. Coverage in the i’th layer at 73.6 K as a function of total coverage, calculated as fractions of a commensurate monolayer. The break in the smooth variation for the monolayer (i = 1) is due to a preference for commensurate density. Note that adsorption into second and third layers occurs simultaneously at N,* > 2, rather than as a simple layer-by-layer filling.
Adsorption of N, on the graphite basal plane time. For example parameters of the pair interaction or the molecule-solid potential can be altered to more realistic values, if such changes are indicated by new evidence. Alternatively, one can study adsorbates that are polar or flexible; certain simple types of surface heterogeneity (such as steps in the surface) can be introduced; adsorption of mixtures is readily handled; and extension of this work to higher coverages can also be done. Time-dependent properties such as molecular diffusive behavior in these films have been obtained from the molecular dynamics simulations, but are not discussed here. One could even relax the assumption of a rigid solid by including the carbon atoms in the dynamics, but this would probably require an inordinate amount of computer time. Another time-consuming improvement would be to handle the many-body terms in the intermolecular interactions explicitly rather than in the handwaving, but economical, way of using “effective” pair-wise models. The attentive reader will note that one important topic has been avoided in this discussion of results. This is the simulation of the adsorption isotherm. The problem here is that the isotherm is essentially a measure of the density dependence of the adsorbate chemical potential, and this is notoriously difficult to evaluate by computer simulation in either the Monte Carlo or the molecular dynamics approach. Specifically, one has
where a,g denote adsorbed and gas phase, respectively, &is the standard-state chemical potential and p is the isotherm pressure at equilibrium with an adsorbed film containing n, molecules per unit area. There are several possibilit&s for evaluating u, by simulation: First, note that E,, the partial molal energy of the adsorbed film, is readily simulated but s,, the partial molal entropy, is not. In principle, the components of the pressure tensor perpendicular and parallel to the surface can be simulated. Although the perpendicular term piZ must be equal to the gas phase pressure, it is generally obtained as a very small difference between two very large terms arising from the attractive and the repulsive parts of the z-component of the force on a molecule in the film. Consequently, the statistical errors in this calculation are often unacceptably large. On the other hand, one could use the spreading pressure to calculate the isotherm from the Gibbs adsorption equation, which we write as
ln(pK,in,)
with
= r&
c&
- 1) 2,
(4.2)
15
and KH is the Henry’s law constant defined by lim (Pins) *.-,I
= KH.
(4.4)
The components of the pressure tensor P,,, P,, and P,, vary from point to point in the fluid near the surface, and thus an average over the entire volume V, of the film is needed in eqn (4.3). Although there is a well-known algorithm to evaluate these pressures in a simulation based on the virial theorem of Clausius, a number of complications arise when attempting to treat a molecular gas in the external field due to the molecule-solid potential, especially if this field has a significant periodic variation in the plane of the surface. (Work on this problem is currently in progress.) Alternatives to this route to the isotherm include various versions of the potential-distribution theory, which are essentially calculations of the work or energy of adding a particle to the fluid phase [ 181. We explicitly mention only one of these which has been employed by Knight and Monson [19] in an adsorption calculation. Although their work is based on the Monte Carlo algorithm, this distinction is unimportant here since both algorithms are capable of generating datasets containing the molecular coordinates for configurations that characterize the statistical ensemble for the adsorbed phase. The potential-distribution calculation relies on the generalized isotherm equation.
$$T)=-(%$ = -In
5 , ( %-Ii
(4.5)
(45)
where &a is the configurational integral for n, molecules near the surface. Knight and Monson then write
z nyZ “PI
h
(exp(-uWWno_,,
(45)
where u(r) is the interaction energy of a molecule at point r with the neighboring adsorbate and with the solid, and p(r) is the local density at that point. The average is taken over the configurations for n, - 1 particles. For this reason, the addedmolecule is often described as a ghost particle, since it does not affect the trajectories or configurations of the n, - 1 others. The computation thus consists in evaluation of the exponential factor in eqn (4.5) for a large number of ghost particles inserted in a large number of molecular configurations. Of course, the result will often be that the ghost experiences a large positive interaction, especially in a dense phase where the probability of overlap with a real molecule is quite large. If the exponential is always close to zero, it is difficult to get good statistics in this region of the isotherm. A final possibility for evaluation of the chemical potential is to recollect the fundamental theorem of
W. A.
16
STEELE
statistical mechanics which states that the fluctuation in the number of molecules N in a volume of fluid is related to chemical potential by 772
-
772
&=-
F
(4.5)
The quantity 6 can be evaluated in a simulation as long as the sample is large enough to eliminate the effects of the periodic boundary conditions which, for example, totally suppress the ~uctuations in a volume coincident with the computer box size, If one can obtain accurate values of S by statistical averaging using small volumes, one can integrate (+./ an,) with respect to n, from the Henry’s law limit up to monolayer or even multilayer coverage. Of course, this necessitates a series of simulations closely spaced
et
al.
over the relevant range of coverage as well as working with a one-phase system, since the fluctuation theorem is no longer useful at densities corresponding to two-dimensional phase coexistence. To date, simulations of adsorption isotherms are almost nonexistent, other than Knight and Monson and some early Monte Carlo studies [20] of rare gases on gcb using the grand canonical ensemble. Still, the range of molecular such properties that can be generated by computer simulation is impressive and enlightening, and it is highly likely that studies of other molecules as adsorbates and other solids as adsorbents will yield much more useful info~ation in the future. Acknowledgments-Partial support for this work was provided by the Division of Materials Research of the National Science Foundation,
REFERENCES
J. K. Kjems, L. Passell, H. Taub, J. G. Dash and A. D. Novaco, Whys.Rev. B 13,1446 (1976); J. Eckert, W. D. Ellenson, J. B. Hastings and L. Passe& P&r. Rev. Lett. 43,1329 (1979); R. D. Diehl, M. F. Toney and S. C. Fain, Jr., Phys. Rev. Letf. 48, 177 (1982); R. D. Diehl and S. C. Fain, Jr., Surf Sci. 125, 116 (1983); S. C. Fain, Jr., M. D. ChinnandR. D. Diehl, Phvs. Rev. B 21. 4170 119801: M. Nielsen. K. Kiaer and J. Bohr, J. &iectron’Spec: Relat. Phenom. 30,‘lll (1983); K. Morishige, C. Mowforth and R. K. Thomas, surf. .Sci. 151, 289 (1985); R. D. Diehl and S. Fain, Phys. Rev. B 26, 4785 (1982); R. D. Diehl and S. C. Fain, Jr., J. Chem. Phys. 77,5065 (1982); H. You andS. C. Fain, Jr.,D&x~s. FaradaySot. 81,159(1985). D. M. Butler. G. B. Huff. R. W. Toth and G. A. Stewart, Phys: Rev. Lett. 34, 1718 (1975); A. D. Migone, H. K. Kim, M. H. W. Chan, J. Talbot, D. J. Tildesley and W. A. Steele, Phyr. Rev. Left. 51, 192 (19831: 0. M. Zhana. H. K. Kim and M. H. W. Chan. Phys.‘Rev. B 33,41~~1986); Q. M. Zhang, H. K. Kim and M. H. W. Chart, Phys. Rev. R 32,182O (1985); M. H. W. Ghan, A. D. Migone, K. D. Miner and Z. R. Li, Phys. Rev. B 30,268l (1984); J. Piper, J. A. Morrison, C. Peters and Y. Ozaki, J. Chem. Sot. Faraday Trans. 179,2863 (1983); T. T. Chung and J. G. Dash, Surf. Sci. 66,559 (1977). J. Talbot, D. J. Tildeslev and W. A. Steele, Mol. Phys. 51, 1331 (1984); J. Taldot, D. J. Tildesley and W. A. Steele. Surface Sci.. 169.71 (1986). A. V. Vernov and W. A.‘Steele, Langmuir 2, 219 (i986); A. V. Vernov and W. A. Steele, Surf. Sci., 171,83 (1986); C. Peters and M. Klein, Mol. Phys. 54, 985 (1985); Y. P. Joslin and D. J. Tildesiey, Mol. Phys. 55,999 (1985); J. Talhot, D. J. Tildeslev and W. A. Steele, Discuss. Fara~v Sot. 81,91(1985j. See, for example, W. B. Streett and D. 3. Tildesley, in Computer Mode~l~ngof Matter (Edited by P. Lykos). A.C.S. Symposium Series, Vol. 86; S. Murad and K.
5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16.
17.
18. 19. 20.
Gubbins, in Computer ModelZingof Matter (Edited by P. Lvkos). A.C.S. Svmuosium Series. Vol. 86: D. J. Evans, kol. Phys. $, 317 (1977). Washington,‘D.C. C. W. Gear, Numerical Init~aIValue Problems in Ord~nary Di~erent~u~Equation. Prentice-Hall, Englewood Cliffs, NJ (1971). D. J. Evans and G. P. Morriss, Comput. Phys. Rep. 1, 297 (1984); Chem. Phys. 87,451 (1984). _ S. Nose. J. Chem. Phvs. 81. 511 (19841. C. S. ikurthv, K. Smger, M. L. Klein and I. R. MacDonald, &iol. Phys,41,1387 (1980). W. A. Steele. J. Phvs. (Paris) C4. 61 (1977). W. A. Steele; .I. P/&s. &em: 82,817 (1978). G. Vidali and M. Cole, Phys. Rev. B 29,6736 (1984). L. W. Bruch, J. Chem. Phys. 79,3148 (1983). 0. Sinanoelu and K. S. Pitzer. 1. Chem. Phvs. 32.1279 (1960); A.“D. MacLachlan, &or. Phys. 7,381 (Ib64). M. J. Bojan, Ph.D. thesis, Penn State University, 1986. A. D. Buckingham, R. L. Disch and D. A. Dunmur, J. Am. Chem. Sot. !Mt,3104 (1968); N. J. Bridge and A. D. Buckingham, J. Chem. Phys. 40,2733 (1964). A. D. Migon< M. H. W. Ghan,.K. J. Niskanen and R. B. Griffiths, J. Phys. (Paris) C 16, L1115 (1983); D. M. Butler, J. A. Litzinger, G. A. Stewart and R. B. Griffiths, Phys. Rev. Len. 42, 1289 (1979). W. A. Steele. The Interaction of Gases with Solid Surfaces, Sec. ‘4.8. Pergamon, New York (1974), sec. 4.8. B. Widom, J. Chem. Phys. 86,869 (1982); K. Mon and R. B. Griffiths, Phys. Rev. A 31,956 (1985). J. F. Knight and P. A. Monson, J. Chem. Phys. 84, 1909 (1986). D. Nicholson, N. G. Parsonage and L. A. Rowley, Mol. Phys. 44,629 (1981); L. A. Rowley, D. Nicholson and N. G. Parsonage, Gol. Phys. 31, k5,389 (1976); J. Commit. Phvs. 26.66 (1978): J. E. Lane and T. H. Spurling, Aus; J. khem. 291.2103 (1976); 31, 933 (1978).
DISCUSSION Question by P. Pendleton, American Cyanamid In these systems, one would expect to observe an induction effect through the graphite substrate as introduced by McLachan. Firstly, for sub-picosec intervals, will future models include these sub-surface
Adsorption of N, on the graphite basal plane atom-mirror image interactions? Secondly, if these terms are included, introduced in a model with step heterogeneities?
how much complication
17 would be
Reply by Steele There is little doubt that McLachan-type substrate-mediated interactions are significant in these systems and indeed, the problem has been directly studied by Bruch (ref. 12). One attempts uses “effective” twobody functions which crudely accounts for these alterations for the simple reason that the full many-body potential would run so slowly in the available computers as to eliminate all meaningful results. We are among those who would very much like to be able to do a better job in this respect. As for the second question, the answer is unknown, but would make a worthwhile thesis project. Question by K. S. W. Sing, Brunel University Could you please comment on the factors which are likely to influence the form of the multilayer isotherm? In our experience the form of the nitrogen multilayer at 77 K (as analyzed by the FHH method) is almost independent of the nature of the substrate surface. Reply by Steele This is a difficult question to answer, especially since simulations for a single substrate do not necessarily help in understanding the differences or similarities with other substrates. The only suggestive feature that comes to mind is that both these and the older Ar/gcb multilayer simulations indicate that several layers are filling simultaneously at total coverages greater than 8 = 1.5 and at temperatures equal to or greater than the bulk boiling point. This smears out the energetics somewhat and might be construed as indicating an entropy-dominated process. Much more work needs to be done to fully answer the question. Question by H. A. Resing, Naval Research Lab Is the failure to observe orientational or translational computer time?
melting of the 2-D solid perhaps due to lack of allowed
Reply by Steele Orientational degrees of freedom come to equilibrium quite rapidly in these systems and no problems of equilibration time are encountered at any temperature of interest. The situation is more complex when translational melting (or motion) is studied. A strong dependence upon temperature and surface density is observed which appears to reflect the fact that molecular interchanges or site-to-site diffusion is more difficult in 2-D than in 3-D. This may indicate that local density fluctuations are smaller because the degree of freedom perpendicular to the surface is strongly hindered. We believe that translational equilibration times have become marginal on the time-scale of the computation for T > 65, and that they become satisfactorily short as temperature increases by another 15”. Question by A. Groszek, British Petroleum Co. What are the densities of adsorbed NZ before and after the melting transition? Reply by Steele The solid layer remains at the commensurate density at temperatures up to the melting point of -47 K for a patch of solid. The patch actually melts to a 2-D vapor with density determined by (number of molecules/ total area available). If the surface is completely filled with solid, melting appears to occur at constant density but at temperatures that are very dependent upon the initial density of the solid. Question by V. Deitz, Naval Research Lab Does the calculation include interactions with basal planes below the surface basal plane? How many basal planes contribute to the energy? Reply by Steele The interactions are summed over all basal planes. The summation converges rapidly even for the relatively long-range attractive terms. In particular, in the region of the minimum energy, the entire contribution of all planes excluding the surface is less than 15% of the total.