Langmuir 1994,10, 1550-1555
1550
The Grand Canonical Ensemble Monte Carlo Simulation of Nitrogen on Graphite E. J. Bottanit and V. A. Bakaev*J Department of Chemistry, The Pennsylvania State University, 152 Davey Laboratory, University Park, Pennsylvania 16802 Received July 14, 1993. I n Final Form: January 26, 1994" The adsorption of Nz molecules on the basal plane of graphite was computer simulated by the canonical (CMC) and grand canonical (GCMC) ensemble Monte Carlo methods. The dependence of the average energy of adsorption on coverage obtained by the CMC method coincides with the results of independent molecular dynamics simulations. The isotherm of adsorption obtained by the GCMC method is close to the experimental data. The small deviations of the simulated isotherm from the experiment in the submonolayer region may be explained by the residual heterogeneity of the real surface of graphite. The deviation in the region of the second layer is probably due to the many-body interactions which were not properly taken into account in the computer simulations. The results for the rigid nonspherical model of Nz are compared to those for a spherical model. 1. Introduction
Computer simulations of the isotherms of adsorption are important for theoretical studies and practical application. Two main techniques used in adsorption studies are the grand canonical ensemble Monte Carlo (GCMC)' and the test particle insertion method.2 Both methods have been applied to spherical molecules, e.g. rare gases. The first GCMC simulation of adsorption of nonspherical molecules appears to be due to Spurling and Lane.3 They specified the orientation of a molecule by two angles and described how to make a creation move and a displacement move. The algorithm was used for simulation of ammonia adsorbed onto a graphite surface. However, they were unable to find parameters for a Lennard-Jones (LJ) potential which reproduded experimental adsorption isotherms. So they limited their simulation to one coverage (ca. 0.75) and studied the structural characteristics (distribution functions) of the adsorbed molecules. Recently there appeared several communications on the GCMC simulations of N24*5 and 0 2 6 on graphite as well as 0 2 in slitlike pores.' In refs 5-7, isotherms of adsorption are presented. However these were not compared against experimental data; instead, the isotherms and distribution functions obtained in the simulations are used to test approximate theories of adsorption of nonspherical moleculesa or to consider the peculiarities of adsorption and capillary condensation in slitlike pores.7 In these papers, there is no description of the GCMC algorithm which refers to the nonspherical character of Nz and 02 molecules. In ref 7, e.g., it is only mentioned that the classical GCMC method was used, as described in detail in ref 1. But in ref 1 (see p 177) the discussion t Permanent address: Instituto de Investigaciones Fisicoqulmicas Tedricas y Aplicadas (INIFTA), C.C. 16, SUC.4. 1900-La Plata, Argentina. 8 Permanent address: Institute of Phvsical Chemistrv. Russian Academy of Sciences, Leninsky Prospect 31, Moscow 117915,Russia. m . Abstract published in Advance ACS Abstracts, April 1, 1994. (1)Nicholson, D.;Parsonage, N. G . Computer Simulation and the Statistical Mechanics of Adsorption; Academic Press: London, 1982. (2) Cheng, A.; Steele, W. A. Mol. Simulation 1990,4, 349. Lane, J. E. Aust. J . Chem. 1980,33, 1967. (3)Spurling, T. H.; (4)Lajtar, L.;Patrykiejew, A.; Sokolowski, S. Acta Phys. Pol. 1987, -A-.7 7- ,. 56.1. - - -. (5) Lajtar, L.; Patrykiejew, A.; Sokolowski, S. J.Chem. SOC.,Faraday Trans. 2 1987,83,473. (6)Patrykiejew, A.; Sokolowski,S. Fluid Phase Equilib. 1989,48,243. (7)Sokolowski, S.Mol. Phys. 1992,75,999.
0743-7463/94/2410-1550$04.50/0
of methods of computer simulation has been concentrated on systems composed of monatomic molecules. The description of Monte Carlo methods for rigid nonspherical molecules is very brief and unspecific. One can find a more detailed discussion of the Monte Carlo algorithms for nonspherical molecules in ref 8. Recently, several papers on computer simulation of phase equilibria of nonspherical molecules and even polymers by the Gibbs ensemble method have appeared (see ref 9 and references therein). The problems and algorithms of computer simulation of phase equilibria are very similar but not identical to those of adsorption equilibria. In this connection, the purpose of this paper is to give a more detailed description of the algorithm of GCMC for rigid nonspherical molecules and to compare the results of computer simulation of adsorption isotherms with experimental data. Specifically we consider here adsorption of nitrogen on the basal plane of graphite. The system has been a subject of numerous computer simulations. The references may be found in the recent papers on the path-integral Monte Carlo simulation of this system.1° The structure of the adsorbate, its average energy, and heat capacity together with other functions of molecular configurations as well as translational and orientational phase transitions have been extensivelystudied by the Monte Carlo and molecular dynamics methods and compared with the relevant experiments. However, all these papers (except those cited above) did not consider the computer simulation of the isotherm of adsorption of nitrogen on the basal plane of graphite. That problem boils down actually to the direct or indirect (e.g., GCMC) calculation of the entropy of adsorption which cannot be calculated as an average value over configurations generated by the canonical ensemble Monte Carlo or molecular dynamics methods. The importance of the entropy of adsorption for the system under consideration consists in the fact that its variation with coverage determines the variation of the system's free energy (the shape of the isotherm of adsorption). The reason for that is that the basal plane of graphite is a homogeneous surface and the adsorption (8) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (9)Gubbins, K.E.Fluid Phase Equilib. 1993,83, 1. (10)Marx, D.;Sengupta, S.; Nielaba, P. J.Chem. Phys. 1993,99,6031.
0 1994 American Chemical Society
Simulation of Nitrogen on Graphite energy varies only slightly in the submonolayer region due to adsorbate-adsorbate interaction so that the main contribution to the change of the free energy (chemical potential) in that region comes from the entropy of adsorption. Thus in what follows, we calculate the isotherm of adsorption of nitrogen on the basal plane of graphite and (for the first time) compare it with experimental data. 2. Grand Canonical Monte Carlo for Rigid Nonspherical Molecules Consider an open system of rigid nonspherical molecules. Its configuration is determined by their number Nand by the set of their coordinates {si], i = 1, ...,N, where each qi consists of s coordinates: three Cartesian coordinates determining the position of the center of a molecule and two (for alinear molecule) or three coordinates determining its Orientation. The average value of any function of configuration is determined in the grand canonical ensemble by the equality " ZN
( F ) = e-'&-J
...
=ON!
Here z = exp(p/kT)/A where p is chemical potential and l l A designates the space and angular momenta partition function.'*" The latter is actually excluded from the computation since z is activity-a measurable quantity. The pressure of a gas is connected with activity through the equation
P
-=
kT
" cbjz3
where coefficients bj may be expressed through the virial coefficients Bj.12 Invert the power series13to obtain
The virial coefficients B2 and Bs may be either calculated directly or evaluated from Tables 1-A,E of ref 14 where the parameters of effective LJ potentials for a number of molecules and virial coefficients corresponding to the LJ potential are compiled. The interaction of nonspherical molecules is approximated by the spherically symmetric LJ potential with effective parameters obtained from the best fit of the experimental values of the second virial coefficients. The latter way is accepted in the current paper as the simplest since we simulated isotherms at low pressures where the second and third terms on the righthand side of eq 3 are only small corrections. For simulated isotherms of adsorption in Figure 3, e.g., the highest correction to z (at P = 700 Torr) due to B2 is -3.3% and due to B3 is -3 X At higher pressures where these terms are comparable in magnitude to the first one, it is better not to use the virial expansion (eqs 2 and 3) at all and describe a system by activity instead of pressure. (11) Gray, C. G.; Gubbins, K. E. Theory of Molecular Fluids. Vol. I : Fundamentals; Clarendon Press: Oxford, 1984; p 157. (12) Hill, T.L. Statistical Mechanics, principles and selected applications; McGraw-Hill: New York, 1956; pp 131, 140, 141. (13) Gradshtein, I. S.;Ryzhik, I. M. Table of Integrals, Series, and Products; Academic Press: New York, 1980; Formula 0.315. (14) Hirschfelder, J. 0.; Curtis, C. F.; Bird, R. B. Molecular Theory of Gases and Liquids; J. Wiley & Sons: New York, 1954.
Langmuir, Vol. 10, No. 5, 1994 1551 If F = 1 and UN = 0 in eq 1, the system is the ideal gas andz = P/kT. It shows that d e is actually dridQJn, where Q is equal to 87r2for nonlinear and 47r for linear molecules.11 Thus integration with respect to angular coordinates should be normalized so that the configurational integral in eq 1 be VN. Then eq 1 holds since E = exp(PV/kT). Multiply z and divide dri by V so that dqilVQ is the probability to choose an ith element of volume from the s-dimensional configuration space for one molecule at random. Suppose also that all these probabilities are equal to l l B , where B is very large but finite number. The value of B is determined by the finite precision of numbers in a computer. Now integrals in eq 1 may be approximated by sums
Here q,j determines the position and orientation of the jth molecule in the ith configuration (allthe configurations are supposed to be numbered somehow). This equation is the starting point for the algorithm first described in ref 15 (see also refs 1 and 8). Formally it is the same algorithm as that for spherical molecules so that one may follow, e.g., the description of algorithm in ref 8. The only difference is the dimensionality of qj which is more than 3 now, so that one has to make a random choice not only of the position of the center of a molecule but also of its orientation in the creation and displacement trials. The latter may be done as in ref 3 but we used another method which takes more random numbers but makes the algorithm of simulation of nonspherical molecules closer to that for spherical ones. We determined the orientation of a linear ith molecule by a unit vector ei directed along the axis of the molecule from the first atom to the second (atoms in a molecule are numbered). In the creation trial this vector should be chosen at random on the surface of a unit sphere (as in ref 3; see also the algorithm of the choice in ref 8, p 349). The coordinates of ei (we used all three of them though only two are independent) together with Cartesian coordinates of the center of a molecule constitute a vector qi which should be chosen at random in the creation trial of the ith particle in the Norman and Filinov algorithm.lS One may use any of three methods described in ref 8 for changing the orientation of a molecule in the canonical Monte Carlo algorithm. We, however, used a slightly different one which better conforms to our choice of the orientation of an inserted molecule. Describe a sphere of a radius 6 around the tip of a vector ei and choose a point inside this sphere at random. The new vector will not be of a unit length so we have to normalize it to the unit length. This algorithm is similar to that for choosing a random point on a surface of a unit sphere by the rejection method (see ref 8, p 349) and to that for choosing a random displacement of the center of a molecule. In particular, one may choose one parameter 6 which determines the range of the orientational as well as translational displacements (if the length of a molecule is taken as a unit of length). The difference to the latter is that the probabilities of the angular displacements are not uniformly distributed around the initial point. This, however, does not violate the microscopic reversibility which is (15) Norman, G. E.;Filinov, V. S. High Temp. (Engl. Transl.) 1969, 7,216.
Bottani and Bakaev
1552 Langmuir, Vol. 10, No. 5, 1994 11
A
A' 1'21
0.8
~
P s
$ 0.6
a
0.4
0.5
i?
4
8
7.5
0.2
'
0
I
0.5
1
I
1.5
2
2.5
-.. 1-
3
Coverage
Figure 1. Average values of potential energy per molecule for N2adsorbed on graphitevs coveragesat 73.6 K: solid line, results of constant temperature molecular dynamics computer simulation;17 symbols, our results (canonical ensemble MC; error bars, nonsphericalmodelof N2;diamonds,sphericallysymmetricmodel of N2; monolayer capacity is ca. 49 molecules). essential for the Monte Carlo algorithms8 because if we designate the probability of angular displacement of the kth molecule from eki to e$ as pij then pij = pji. 3. Results and Discussion 3.1. Rigid Nonspherical Model of N2. In the simulation of the adsorption of N2 on the basal plane of graphite, the interaction between N2 molecules was assumed to be exactly as in refs 16and 17. Two interaction sites were located at N atoms (0.11 nm apart) so that LJ interaction of two N2 molecules was approximated by four LJ 12-6 potentials with parameters e"/k = 36.4 K and UNN = 0.332 nm. The electrostatic interaction was approximated with the help of a positive charge 12.98 X 10-20 C at the center of symmetry of the molecule and two negative charges of half of that value at the nitrogen atoms. Some discussion of this adsorbate-adsorbate interaction potential will be given below. The N2-graphite energy (adsorption potential, AP) also was calculated almost as in ref 17. The laterally averaged part of AP was calculated as in ref 18
Here h is the distance of a nitrogen atom from the uppermost basal plane of graphite and h; = h + id, where d is the distance between planes. The parameters in eq 5 had the following values: t/k = 31.5 K; u = 0.336 nm; p a d = 4.313 and d = 0.335 nm and the summation was limited to ten first terms. The part of AP which depends on lateral coordinates was calculated as in ref 19, that is the Fourier components were multiplied by the factor 1.5 to increase the corrugation of the surface. We have found that inclusion of the first set of Fourier component terms in AP had a negligible effect on the results of our caleulations. Thus all the results in Figures 1and 2 have been obtained with a laterally independent AP (calculated by eq 5) and those in Figure 3 with the laterally dependent AP. (16) Murthy, C. S.;Singer, K.; Klein, M. L.; McDonald, I. R. Mol. Phys. 1980,41, 1387. (17) Vernov, A. V.;Steele, W. A. Langmuir 1986,2, 219. (18) Steele, W.A. The Interaction of Gases with Solid Surfaces; Pergamon Press: Oxford, 1974; p 53. (19) Kim, H.-Y.; Steele, W. A. Phys. Rev. B 1992, 45,6226.
I,-
I
0
2
0
I
I
4
6
0
10
Pressure (Torr) Figure 2. Isotherm of adsorption of N2 on graphite (submonolayer region): solid line,experimentaldata22at 77.8 K; diamonds, experimentaldata;23dashed line, experimentaldataz2at 90.1 K error bars, GCMC computer simulation for the nonspherical model of Nz; squares, GCMC computer simulation for the spherically symmetric model of Nz at 90.1 K. 40
4
5 0
100
200
300
I
400
500
600
700
Pressure (Torr) Figure 3. Isotherm of adsorption of Nz on graphite at 77.8 K
(multilayer region): solid line, experimental data;22 error bars, GCMC computer simulation for the nonspherical model of Nz; diamonds, the same as the error bars but the well depth of the N-N interaction potential is diminished by about 5%; squares, GCMC computer simulation for the spherically symmetricmodel of Nz. In Figure 1, the energies of adsorption obtained in canonical ensemble Monte Carlo (MC) computer simulations are compared with the results of the independent constant temperature molecular dynamics computer simulations from ref 17. Our results were based on programs described earlier which allow canonical as well as grand canonical MC computer simulations.20 In particular, error bars in Figures 1-3 correspond to the method of estimation of errors described in ref 20. The difference was that now the trial shift of the position of a particle consisted not only of a random choice of a point in a small cube with the length of edge 6 but also of a random choice of a point in a small sphere whose diameter we also took equal to 6. Those six random numbers determine the new position and orientation (as was described above) of a molecule in a trial test. In our simulations 6 was 0.09~" and the ratio of successful trials to all trials was about 0.4-0.6 for the shift trials and an order of magnitude less for the creation and destruction trials. Periodic boundary conditions were imposed in the lateral directions. Periods in x - and y-directions were q x and qy correspondingly. We tried (20) Bakaev, V. A.; Steele, W. A. Langmuir 1992,8, 148.
Langmuir, Vol. 10, No. 5, 1994 1553
Simulation of Nitrogen on Graphite several simulation boxes (qL.= SUNN, q y = qx or qy = ((31f2/ 2)qJ and one close to that used in ref 17) but results were essentially the same. In the vertical direction the boundary condition was that of a mirror reflection from the ceiling of the rectangular simulation box. The height of the box provided negligible (almost gas) density of adsorbate at the ceiling. In Figures 1 and 2, it was about 3am. In Figure 3, it was determined by the following way. First, the effective thickness of the adsorption layer was determined by dividing the number of adsorbed molecules by the density of bulk liquid and by the surface area. The height of the simulation box was established as the effective thickness of the adsorption layer plus 3am. The small height of the simulation box and the low density of gas phase in our simulations allow us to consider all molecules in the simulation box as adsorbed. The difference between the value of adsorption determined in this way and that determined as the Gibbs excess' is negligibly small in comparison to the errors of computer simulation. Finally, the nominal monolayer capacity was chosen as that of a commensurate monolayer (as in ref 17) which corresponds to the molecular area of N2 equal to 0.157 nm2. The general conclusion that may be drawn from Figure 1is that our simulations give essentially the same results as the totally independent method of molecular dynamics computer simulations which in turn are rather close to experimental data (cf. Figure 7 from ref 17). At the completion of the monolayer, the experimental curve of isosteric heat of adsorption vs coverage for N2 on graphite displays a sharp spike.21 Similar spikes were observed for adsorption of rare gases on graphite in real experiments and computer simulations (see ref 1, pp 241-265). These were explained either by a two-dimensional "crystallization"l or by phase transitions connected with variation of adsorption potential in the lateral directions,' e.g., transition from a two-dimensional fluid phase to the registered 4 3 X 4 3 two-dimensional ordered phase.21 In our case there was no variation of adsorption potential in the lateral directions. Nevertheless,we observed on snapshot pictures structure which displayed a kind of crystalline order. Besides, the point near maximum of our curve of the integral energy of adsorption which drops below the general trend is not accidental. It was repeated several times in different simulations with various simulation boxes and initial conditions. This point gives a spike in the derivative of a curve. However, we have not obtained unambiguous data for the detailed behavior in this region and the problem requires additional investigation. A comparison of experimental and simulated isotherms of adsorption is presented in Figures 2 and 3. The experimental isotherm is taken from the review by Avgul and Kiselev.22 It is a compilation of experimental data from various sources. In general, the simulated isotherm is close to the experimental one. However, there are some deviations. The deviation of experiment from the calculation (for the nonspherical model) at low coverage (Figure 2) probably results from the residual heterogeneity of the surface of the real graphitized carbon black. The existence of such a heterogeneity is clearly seen from Figure 7 of ref 17 where the same curve as in Figure 1 is compared to experimental data. A t coverages below ca. 0.7 the experimental energies are higher than the calculated ones because the real surface of graphite contains residual (21) Piper, J.; Morrison, J. A.; Peters, C.; Ozaki, Y. J. Chem. Soc., Faraday Trans. 1 1983, 7, 2863. (22) Avgul, N . N.; Kieelev, A. V. Chemistry and Physics of Carbon; Marcel Dekker: New York, 1970; Vol. 6.
Table 1. Differential Heats of Adsorption (kJ/mol) of Nitrogen on a Basal Plane of Graphite at ca. 80 K a,,, = 10.26 umol/m2
0.4 1.0 2.0 3.0 4.0 5.0 6.0 7.0
0.04 0.10 0.20 0.29 0.39 0.49 0.59 0.68
9.33 9.58 9.96 10.4 10.6 10.8 11.0 11.2
9.21 9.37 9.65 9.92 10.2 10.5 10.7 11.0
8.73 9.13 9.61 9.81 9.98 10.2 10.5 10.2
heterogeneity as against the ideally homogeneous surface of computer simulations. The experimental isotherm of adsorption measured on a more homogeneous surface23is closer to our simulated isotherm than that from ref 22. The isotherms of adsorption (for the nonspherical model) in Figure 2 make it possible to determine the isosteric heat of adsorption q,t a t an intermediate temperature ca. 84 K. The differential heat of adsorption qd which is actually the differential energy of adsorption -dU/ CW may be determined from Figure 1. The two heats of adsorption are connected by the thermodynamic relation qat = qd + RT which may be employed as a test of correctness of the GCMC simulation.20 In Table 1 the differential heats of adsorption determined by the two methods are compared to each other and to the experimental values from ref 22. The differential heats of adsorption determined in computer simulation by two intrinsically different methods (two last columns of Table 1) are reasonably close to each other. The differential heat determined as differential energy is slightly larger than that determined from qat which may be partly explained by the fact that the former was determined at the lower temperature. The experimental data in the ref 22 were compiled from different sources. Experimental error are not indicated but one may assume that they hardly are below the 5% level (due to scatter of sample properties etc.). Thus the agreement between experiment and computer simulation is quite reasonable. In the multilayer region, the simulated isotherm of adsorption is quite sensitive to the N2-N2 interaction. Our main simulated isotherm (error bars connected by dashed curve in Figure 3) was obtained with the model X1 intermolecular potential of nitrogen determined in ref 16. This is a two-center LJ 12-6 potential supplemented by the interaction between quadrupole moments. Parameters of that potential (including the quadrupole moment which differs slightly from that of a nitrogen molecule) were obtained by fitting the lattice spacing and sublimation energy of the crystal at absolute zero as well as the second virial coefficient of the gas. This model potential is successful in accounting for a wide variety of the properties of nitrogen, gaseous, liquid, and crystalline.16 It is seen from Figure 3 that the X1 model potential of nitrogen-nitrogen interaction is also successful in accounting for the adsorbate-adsorbate interaction since the simulated isotherm of adsorption is rather close to the experimental one, especially a t high coverages. At coverages above two layers (pressures above 350 Torr), the deviations of points of simulated isotherm from the experimental one are less than the errors of computer simulation which increase considerably in this region. Even in the region of the largest deviations, the agreement must be considered as good in view of the fact that even a small variation of the parameters of the interatomic potential ~~
~
(23) Cascarini de Torre, L. E.; Llanos, J. L.; Bottani, E.J. Collect. Czech. Chem. Commun. 1988,53,251.
1554 Langmuir, Vol. 10, No. 5, 1994 gives a rather large deviation of the simulated isotherm of adsorption from the experimental one. For example, the diamond points in Figure 3 were obtained with the same parameters of the nitrogen-nitrogen potential as the points marked by error bars except that the t of the LJ potential was decreased by 5 % (from 36.4 K in the latter case to 34.6 K in the former one). On the other hand, the isotherms in Figure 3 are expressed as adsorption versus absolute pressure. When a similar comparison between simulated and experimental isotherms of adsorption has been made for the case of adsorption of argon on the model amorphous ~ x i d e ,the ~ ~agreement ,~~ was not so good. It improved when isotherms were presented as adsorption vs relative pressure and the saturated vapor pressure p, were taken from computer simulation of a model Ar with the same interatomic potential as in simulation of adsorption. Thus the good agreement between the high coverage simulated and experimental isotherms of adsorption in Figure 3 suggests that the X1 model potential of N2-N2 interaction16 should be successful in a calculation of the saturated vapor pressure of liquid nitrogen at this temperature. A possible explanation for the observed small difference between the experimental and the simulated isotherms is that the many-body terms that are buried in the effective two-body X1 N2-N2 potential may be modified when the interacting pair is in a thin adsorbed film rather than in the homogeneous bulk film. Note that there are two phase boundaries in the adsorbed system: film-vapor and filmsolid. The many-body sums for pairs of adsorbed molecules will be modified by the presence of each of these. Thus, the effective two-body potential will be changed, but magnitudes and even the signs of the changes are not known in the second adsorbed layer (only the changes in the monolayer have been carefully considered to date). 3.2. Spherically Symmetric Model of Nitrogen Molecule. The rigid nonspherical model of nitrogen molecule has 2 degrees of freedom more than a simplest spherically symmetric model. These additional degrees of freedom for nitrogen on graphite provide a variety of properties absent for adsorption of monoatomic molecules like argon (whose adsorption behavior on graphite is generally very close to that of nitrogen). These properties (orientational phase transitions etc.) have been extensively studied in computer simulations (see, e.g., ref 10 and references therein). However, there exists one other, still unexploited, way to evaluate the influence of these additional degrees of freedom on adsorption phenomena-a direct comparison of results of computer simulation for spherical and nonspherical models of a nitrogen molecule. In this section we describe the first results of such a comparison. We characterize a spherically symmetric model of a nitrogen molecule by four parameters: two LJ parameters E , and a, describing interaction of nitrogen with graphite and two W parameters t and u describing interaction between nitrogen molecules. The values of t may be determined from the values of critical temperature of N2 (T,= 126.2 K26)and the value of the reduced critical temperature of the LJ liquid (TJt = 1.3327). So e = 95.0 K which coincides with the value determined from the second virial coefficient14and with (24) Bakaev, V. A.; Voit, A. V. Izu. Akad. Nauk SSSR,Ser. Khim. 1990,2007; Bull. Acad. Sci. USSR, Diu. Chem. Sci. (Engl. Transl.) 1991, 39, 1822. (25) Bakaev, V. A.; Steele, W. A. Langmuir 1992, 8, 1379. (26) Jacobsen,R. T.;Stewart, R. B.; Jahangiri, M. J.Phys. Chem. Ref. nntn -- - -1968. - - -, -1.5.-, 7.15. . - -. (27) Adams, D.J. Mol. Phys. 1976,32,647.
Bottani and Bakaev the value accepted in the lattice model of nitrogen adsorption.' The value of a may be determined in several ways. It may be obtained from the value of molecular area of nitrogen on graphite (16.2 A2)22taking it as an area per atom with the van der Waals diameter 2% in the densest packing of the atoms on a plane. Determined in this way, a = 3.85 A. Another way is to compare the saturated vapor pressure of Nzp, at, e.g., 78K26with the reduced saturated vapor pressure of LJ liquid psa3/t at the reduced temperature Determined in this way, a = 3.990A. The difference between two values of a is only 3% which, however, gives the difference in p, ca. 10% . We choose here a = 3.990 A. (The value determined from the second virial coefficient or viscosity is ca. 3.7 A.14) The values of E, and a, depend on temperature. Near 80 K the nitrogen molecules lay predominantly flat on the graphite basal planell so that a reasonable value of a, may be the same as that for a N atom and the value of tstwice as large. However, this choice yields the average energies of adsorption about 2.5% larger than those for the nonspherical model. The reason, probably, is that nitrogen molecules do not lay exactly flat on the graphite surface at about 80 K. Thus, we diminished t, so that the final values of the LJ parameters for the nitrogedgraphite interaction are e, = 61.4 K and a, = 3.36 A. The energies of adsorption for the spherical model of a nitrogen molecule are very close to those for the nonspherical model in the submonolayer region at 74 K (Figure 1). However, the isotherms of adsorption differ considerably in the same region at 90 K (Figure 2). Likely reasons may be the difference in the entropies of adsorption of the two models or the dependence of the effective parameters of the spherical model on temperature. Those parameters were chosen to match the energy of adsorption at 74 K (cf. Figure 1). The energy of adsorption for the spherical model decreases with an increase in temperature considerably less than that of the nonspherical model (this was checked by direct calculation). Thus at 90 K the energy of adsorption for the nonspherical model should be less than that for the spherical model which may be the reason for the deviation of the isotherms for both models on Figure 2. In the multilayer region (Figure 31, one again may notice the appreciable difference between the isotherms of adsorption for the two models. The isotherm for the spherical model lies below the isotherm for the nonspherical model. (In the submonolayer region, the relation of the two isotherms is opposite; cf. Figure 2.) Note that parameters of the effective LJ potential describing interaction of the spherically symmetric model of N2 with the graphite surface were chosen on the assumption that a nonsymmetric molecule lies flat on the surface. In the multilayer region, this does not take place. In effect the a parameter should depend on the distance from the surface. The entropy difference between adsorbates of symmetric and nonsymmetric model molecules also may be the reason for the above mentioned difference of isotherms. Spherically symmetric models of nonspherical molecules have been widely used in the theory of adsorption. A particular (even more simplified) case of a spherical model is a lattice gas (Ising) model which has been used in the comparative study of the adsorption isotherms of Ar and N2 on graphite (see ref 1, p 334). The results presented in Figures 1, 2, and 3 directly show to which extent the nonsphericity of a nitrogen molecule may influence its main-adsorption characteristics. However,-wecannot give
Simulation of Nitrogen on Graphite
an unambiguous explanation of this influence now. Probably, the direct calculation of the entropies and heats of adsorption in simulations with variable parameters of the spherical model may be helpful for that matter. Conclusion The GCMC computer simulation of adsorption isotherm for nitrogen is no less reliable than that for rare gases. The time of computation for the nonspherical model is about 6 times larger than that for the spherical one. Provided the adsorbate-adsorbent and adsorbate-adsorbate interactions are modeled properly, the simulated isotherm of adsorption should reproduceexperimentaldata. However, the potential of adsorbate-adsorbate interaction should take into account the many-body interactions, mainly, the mediation of the adsorbate-adsorbate interaction by the substrate. This effect somewhat distorts the simulated isotherm in the second layer. In higher layers it is weak by itself, and in the first layer where the effect is stronger than in the second layer, the adsorbate-adsorbate interaction as a whole plays a relatively minor role in deter-
Langmuir, Vol. 10, No. 5, 1994 1555 mining the shape of the isotherm of adsorption. In general, GCMC computer simulation of the isotherm of adsorption proved to be effective up to the relative pressures of about 0.9. A t higher relative pressures the errors of computer simulation increase drastically. A preliminary comparison of adsorption properties for the rigid nonspherical and spherically symmetric models of Nz has shown that the parameters of the spherical model may be chosen very easily to reproduce the energy of adsorption. However, the isotherms of adsorption for two models are considerably different. The extent of that difference may, however, depend on coverage and temperature. An understanding of these effects may be gained in comparative computer simulations of the entropies of adsorption of spherical and nonspherical model molecules.
Acknowledgment. The authors thank W. A. Steele for discussion of this paper. Support for this research was provided by Grant DMR 902 2681 of the NSF Division of Material Research. E.J.B.is a researcher of the Comisi6n de Investigaciones Cientificas de la Pcia. de Buenos Aires.
