dynamical
Monte Theoretical
foundations
of Carlo simulations
Kristen A. Fichthorn
Department of Chemical Engineering, Pennsylvania State University, University Park, Pennsyivania 16802
W. H. Weinberg Department of Chemical Engineering, University of California, Santa Barbara, California 93106
(Received 12 December 1990;accepted 10 April 1991) Monte Carlo methods are utilized as computational tools in many areasof chemical physics. In this paper, we present the theoretical basis for a dynamical Monte Carlo method in terms of the theory of Poissonprocesses. show that if: ( 1) a dynamical hierarchy of transition We probabilities is created which also satisfy the detailed-balancecriterion; (2) time increments upon successfuleventsare calculated appropriately; and (3) the effectiveindependence of various eventscomprising the system can be achieved,then Monte Carlo methods may be utilized to simulate the Poissonprocessand both static and dynamic properties of model Hamiltonian systemsmay be obtained and interpreted consistently.
1. INTRODUCTION
Monte Carlo methods are utilized as computational tools in many areas of chemical physics. Although this * techniquehas been largely associatedwith obtaining static, or equilibrium properties of model systems, Monte Carlo methodsmay also be utilized to study dynamical phenomena. Often, the dynamics and cooperativity leading to certain structural or configurational properties of matter are not completely amenableto a macroscopiccontinuum description. On the other hand, molecular dynamics simulations describingthe trajectories of individual atoms or molecules on potential energy hypersurfacesare not computationally capableof probing large systemsof interacting particles at long times. Thus, in a dynamical capacity, Monte Carlo methodsare capableof bridging the ostensibly large gap existing between these two well-established dynamical approaches, sincethe dynamics of individal atoms and molecules are modeled in this technique, but only in a coarse-grainedway representing average features which would arise from a lower-level result. The application of the Monte Carlo method to the study of dynamical phenomenarequiresa self-consistentdynamical interpretation of the techniqueand a set of criteria under which this interpretation may be practically extended. In recent publications,3t4 certain inconsistencies have been identified which arise when the dynamical interpretation of the Monte Carlo method is loosely applied. These studies have emphasizedthat, unlike static properties, which must be identical for systems having identical model Hamiltonians, dynamical properties are sensitive to the manner in which the time seriesof eventscharacterizing the evolution of a systemis constructed. In particular, Monte Carlo studies comparing dynamical properties simulated away from thermal equilibrium have revealeddifferencesamong various sampling algorithms.3 These studies have underscoredthe importance of utilizing a Monte Carlo sampling procedure in which transition probabilities are based on a reasonabledynamical model of a particular physical phenomenon under consideration, in addition to satisfying the usual criteria for thermal equilibrium. Unless transition probabilities can be formulated in this way, a relationship
1090 J. Chem. Phys. 95 (2), 15 July 1991
between Monte Carlo time and real time cannot be clearly demonstrated.In many Monte Carlo studies of time-dependent phenomena,results are reported in terms of integral Monte Carlo steps,which obfuscatea definitive role of time. Ambiguities surrounding the relationship of Monte Carlo time to real time precluderigorous comparisonof simulated results to theory and experiment, needlesslyrestricting the technique. Within the past few years, the idea that Monte Carlo methods can be utilized to simulate the Poissonprocess has been advanced in a few publications* and some Monte Carlo algorithms which are implicitly basedon this assumption have been utilized.lS4 This is an attractive prospect, since within the theory of Poisson processes, relathe tionship between Monte Carlo time and real time can be clearly established. In this paper, we shall focus on dynamical interpretation of the Monte Carlo method. We shall show that if three criteria are met, namely, that transition probabilities reflect a dynamical hierarchy in addition to satisfying the detailed-balancecriterion, that time increments upon successful events are formulated correctly in terms of the microscopic kinetics of the system, and that the effective independenceof various events can be achieved, then the Monte Carlo method may be utilized to simulate effectively a Poisson process.Within the theory of Poisson processes, both static and dynamic properties of Hamiltonian systems may be consistently simulated with the benefit that an exact correspondence betweenMonte Carlo time and real time can be establishedin terms of the dynamics of individual species comprising the ensemble.We shall demonstratethe formalism by considering the approach to and the attainment of Langmuir adsorption-desorption equilibrium in a latticegassystem. We shall also discussstraightforward extension of the methodologyto more complicated systemsof interacting particles.
II. DYNAMICAL INTERPRETATION CARLO METHOD OF THE MONTE
Under a dynamical interpretation, the Monte Carlo method provides a numerical solution to the Master equation
0 1991 American institute of Physics
0021-9606191 /I 41 OSO-07$03.00
Downloaded 23 Nov 2005 to 128.125.4.122. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
K. A. Fichthorn and W. H. Weinberg: Dynamical Monte Carlo simulations
1091
-mu,t) = 2 W (cr -u)P(a - c W (u-+o ,t) )P(a,t), dt - a
(1) where u and uare successive statesof the system,P( u,t) is the probability that the system is in state u at time t, and W( u-+a) is the probability per unit time that the system will undergoa transition from state u to state u. The solution of the Master equationis achievedcomputationally by choosing randomly among various possibletransitions to a model system and acceptingparticular transitions with appropriate probabilities.Upon eachsuccessful transition (or, in some instances,each attempted transition), time is typically incremented in integral units of Monte Carlo steps which are related to some unit time r. At steady state, the time derivative of Eq. ( 1) is zero and the sum of all transitions into a particular stateu equalsthe sum of all transitions out of state u. In addition, the detailed-balance criterion W(u +o)P(u ,eq) in which = W(u-+u )P(u,eq),
(2)
which is basedon its rate and is independentof the events before time t. Let us consider the ramifications of these premisesfor a stationary processwith two statesrepresenting forward and reversetransitions. This processmay correspond,e.g., to the time-dependentoccupancy or vacancy of one site among many in the adsorption-desorption equilibrium of a gas-phase moleculewith a solid, singlecrystalline surface. Adopting a frequency definition of probability, the average rate t of, say, the forward transition (which becomes intermittently available via the reverse transition) can be interpreted as a time density of events. Sampling small, identical time intervals S of a larger time increment,t = na, the average is the ratio of the number rate of time intervals containing eventsns to the total number of intervals sampled n per unit time S in the limit S-0 and n-co, lim 2. (4) &.o,r-cc t In the limit S-O, each interval will contain, at most, one event. Also, consistent with our premise is that each time interval has an equal probability r8 of containing an event. Let Ne,tbe a random variablecounting the number of events which have occurred within a time t. Then, the probability that n, eventswill occur in time t is (l , PWe,t =n,) = f- (rS) -rS)- 0c and in the limit S + 0,
f (N,, =n,) =(rt) e-. n,!
r=
P( a,eq ) = Z - - H() e kBT,
(3)
must be imposed for each pair of exchanges,so that the Monte Carlo transition probabilities can be constructed to guaranteethat the system will attain a thermal equilibrium consistentwith the model Hamiltonian. In Eq. (3), Zis the partition function and N is the Hamiltonian of the system. The detailed-balance criterion doesnot, however, uniquely specify theseprobabilities. When static propertiesof model Hamiltonian systemsare sought, Eqs. (2) and (3) are the only standardswhich must be met in addition to the necessity of a sampling procedurewhich is sufficiently random to prevent statistical bias. Dynamical propertiesrequire that a more definite relationship between the Monte Carlo time We step and the transition probabilities is established. shall and show that this relationshipcan be established implemented, oncetransition probabilitiesare formulated as rateswith physical meaning,through the theory of Poissonprocesses. In a dynamical interpretation of the Monte Carlo method, it can be assumed time resolutionis accomplished that on a scaleat which no two eventsoccur simultaneously.Once this perspectivehas been adopted, the task of the Monte Carlo algorithm is to createa chronologicalsequence disof tinct eventsseparated certain interevent times. Sincethe by microscopic dynamics yielding the exact times of various eventsare not modeledin this approach,the chain of events and corresponding interevent times must be constructed from probability distributions weighting appropriately all possibleoutcomes.The distributions governing transitions and intereventtimes availableto a systemat any time t can be developedfrom considerationsfully consistentwith the mesoscopic genre. On a course-grained,mesoscopiclevel, it must be assumedthat the totality of microscopic influences underlying various transitions of a system dictate certain distinctive eventsE={e,,e2,...,en 1, which can be characterized by averagetransition rates Rsir,,r,,...,r,, 1. In the absenceof microscopicdetail, it can be held that any particular transition which becomespossibleat time c can potentially occur at any later time t + At with a uniform probability
(5)
From Fq. (6), the expected number of events occurring within a time t is (N,,,) = rt, from which the rate is recovered through division by t. In Eqs. (4)-(61, it is seenthat a mathematicaladaptionof thesevery basicassumptions leads to the characterizationof a stationary seriesof random, independenteventsoccurring with an average rater in terms of a 9*10 Poissonprocess. It can be shown that a Poissonprocessis consistent with the Master equation. Additional features are attributable to the Poissonprocessand the most significant of thesefor the purposeat hand is characterizationof the probability density of times t, between successive events (7) From the probability density, the mean time period between successive eventsis calculatedas (t, ) = l/r. A particularly useful feature of the Poissonprocessis that an ensemble independentPoissonprocesses beof will have as one, large Poissonprocesssuch that statistical properties of the ensemble be formulated in terms of the dycan namics of individual processes. ConsideringN-independent forward-reverse Poisson processes(or, keeping with the previousanalogy,the adsorption-desorptionequilibrium of an entire systemof independentmolecules),eachwith some arbitrary, but finite rate ri, let N,,, be a random variable counting the overall number of eventsin the ensemble which have occurred within a time interval t. The quantity N,,, is then the sum of random variables counting the number of f,,(t) = re-.
J. Chem. Phys., Vol. 95, No. 2,15 July or copyright, see http://jcp.aip.org/jcp/copyright.jsp Downloaded 23 Nov 2005 to 128.125.4.122. Redistribution subject to AIP license 1991
1092
K. A. Fichthorn and W. H. Weinberg: Dynamical Monte Carlo simulations
events which have occurred in each of the individual processes, i.e., Neat = 5 New
i=l
(8)
The overall probability of n, eventsin time t is given by the convolution of the individual probability mass functions characterizing each of the individual processes RN,,, = n,1 = P(N,,,,)*P(N,,,,)*...*P(N,,,) , (9) and the probability massfunction characterizing the overall distribution of eventsis obtainable by the method of characteristic functions *
r
Here, 19 the fractional surface coverageof A. The absence is of adsorbate-adsorbateinteractions in our example allows Eq. ( 12) to be solved readily for the approach to and the attainment of adsorption-desorption equilibrium once r, and rD are known in terms of the intensive and extensive properties of the system. With the initial condition @t=O) =o,
D A and, in the limit as t+ ~0,
e(t) =
r yr
(l-e-~+rq
(13)
(14)
rA + rD
P(N,,,
=n,)
=@& n,!
--It ,
where A= 2 ri.
i= I
(11)
A final point which should be emphasizedis that the basicpremisesleading to Eqs. (6)) (7), and ( 10) are equally applicableto systemswhich are nonstationary and evolving toward equilibrium. In the nonstationary Poisson process, the overall rate simply becomesa function of time. Thus, if the Monte Carlo algorithm can be madeto simulate the Poisson process,then the relationship betweenMonte Carlo time and real time can be given a firm basis in both static and dynamic situations. In the following section, we shall demonstrate, through an example of the approach to the attainment of Langmuirian adsorption equilibrium, the criteria which must be applied so that the Poisson processmay be effectively simulated with Monte Carlo methods.
III. ADSORPTION-DESORPTION EQUILIBRIUM
Adsorption equilibrium of a gas-phase speciesA with a solid, single-crystalline surface occurs when the chemical potentials of gas-phase chemisorbedA are equal. From a and kinetic point of view, adsorption equilibrium (steady state) is establishedwhen the net rate of chemisorption of gasphaseA is equal to the net rate of desorption of chemisorbed A to the gasphase.Within the context of a lattice-gas model in which adsorption is unactivated, each chemisorbed A molecule requires one adsorption site. We assumethat chemisorbed A molecules do not interact appreciably with one another and envision a collection of gas-phasemolecules whose temperature, pressure, and intrinsic partition functions dictate a seriesof independentarrivals of moleculesto a surface containing a uniform and periodic array of adsites. The arrivals occur at random, uncorrelated times and can be characterized by an average rate r,. A similar scenario is applicableto moleculeschemisorbedon the surface-the totality of microscopic influences (e.g., surface phonons and electron-hole pair creation) acting on an ensembleof chemisorbed moleculesinduces desorption events which occur with an averagerate r, . The appropriate kinetic expression for this balanceis
-$-=r,(l-0) -r,t9. (12)
Equation ( 14) reflects the detailed balance of the simple adsorption-desorption model. If the rates are cast in terms of appropriate partition functions, 5 then both the desired 3- kinetics and thermal equilibrium of the system are ensured. To maintain generality, we shall retain generic rate expressions for the elementary steps. It should be stressed,however, that both kinetic and equilibrium behavior can and should be incorporated in the rates of elementary stepsindigenous to a particular system, so that both aspectsof the system can be modeled consistently. Figure 1 depicts the general features of one algorithm for simulating as a Poissonprocessthe adsorption equilibrium of a gas-phasespeciesA with a two-dimensional lattice containing N sites.A trial in this algorithm beginswhen one of the N sites is selectedrandomly. If the site is vacant, adsorption occurs with probability W, ; and desorption occurs with probability W, if the site is occupied. Time is advanced by an increment ri upon successfulrealization of an event at trial i and, for illustrative purposes,we shall also count the overall number of trials T, which accumulate over repetition of the algorithm. We shall show that through a proper definition of W, and W, , the utilization of an appropriate ri , and the random selectionprocess,the Monte Carlo algorithm of Fig. 1 simulatesthe Poissonprocessand provides the correct
FIG. 1. A flow diagram for simulating as a Poisson process the approach to and the attainment of Langmuirian adsorption equilibrium. 7 the (inteis gral) number of trials, f represents real time, ris a uniform random number between 0 and 1, W, is the transition probability for event i [i = A (adsorption) or D (desorption) 1, and rr is the real time increment at trial T.
J. Chem. Phys., Vol. 95, No. 2.15 July 1991
Downloaded 23 Nov 2005 to 128.125.4.122. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
K. A. Fichthorn and W. H. Weinberg: Dynamical Monte Carlo simulations
1093
solution to Eq. (12) for the N-site ensemble. The first criterion which must be met to simulate effectively the Poissonprocessis that the transition probabilities W,, and W, must be chosenso that the Monte Carlo simulation obeysdetailed balance,as specified,e.g.,by Eq. ( 14). To demonstrate the manner in which this is achieved,let us first consider the discrete stochastic processof the Monte Carlo algorithm of Fig. 1. By performing this algorithm, we simulate a sequenceof independent Bernoulli trials in which the probability per trial of a successfuladsorption event is WA ( 1 - 19~ the probability per trial of a successfuldesorption ), event is W, 8, and the total probability of success trial is per WA ( 1 - 6Ji + W, 19~ 1. Here, 19~ the fractional cover) < is ageof A at trial i, which for Mi occupied lattice sites is given by ei = M,/N. When the system has reached steady state, (8, ) + t9, the simulated equilibrium fractional surfacecoverage of A. In general, 0, will be different from 0, the continuum equilibrium fractional surface coverage,since the continuum coverageshould arise from time-weighted (and not trial-weighted) Bi. However, in systems which are sufficiently large, 8, should approach 8, (uide infra) . The statistics of the equilibrium system may be obtained readily. Let NA,r be a random variable counting the number of successful adsorption eventsin T trials. Then the averageprobability of nA successfuladsorption events in T trials is given by P(N,,,,=n,)=(T/n,)[W,(l--8,)] (15) for O(n, <T. From Eq. ( 15), the expected number of adsorption events in T trials is given by (N~.~) = w,u -e,u-. (16) Similar results may be obtained for desorption events and the overall statistics of both adsorption and desorption. When steady state has occurred in the simulation, the averagerate of adsorption [obtained from Eq. ( 16) ] is equal to the averagerate of desorption, i.e., W,(I -es) = w,e,. (17) Detailed balance is satisfied at equilibrium (steady state) if WA and W, are defined in a way which allows the physical model represented by Eq. ( 14) to be recovered from Eq. ( 17). The transition probabilities satisfying this relationship are not unique and, as noted previously, in the traditional application equilibrium of the Monte Carlo method, these probabilities are often formulated without regard for the behavior of the system away from equilibrium. To simulate dynamical phenomena, an additional criterion is necessary so that transition probabilities reflect unique transition rates (and, hence,simulate dynamics). Theseprobabilities should be formulated so that a dynamical hierarchy of transition rates is establishedin terms of appropriate models for the rates of microscopic events comprising the overall process. Generally stated, a dynamical hierarchy of transition probabilities is created when these probabilities are defined, for a transition i, as wi = ri/&., 9 (18) where ri is the rate at which event i occurs and gm,, &sup {ti}. A dynamical hierarchy is not achieved,e.g., in the stan-
dard Metropolis algorithmI applied to systems approaching equilibrium, becauseall transitions of the system to lower or equivalent energy states are considered to have a probability of unity. The Kawasaki transition probabilities -19 on the other hand, create a dynamical hierarchy among transition rates. However, it has been pointed out3s4 that this hierarchy is not appropriate for most physical processes. the algorithm of Fig. 1, e.g., the transition probaIn bilities could be constructed through normalization of the rates of adsorption and desorption by the larger of the two. If, say, r, > rD, then transition probabilities could be defined
WA=1
and
W,=z,
rA
(19)
X[I-
wA(i-es)lT-+,
(i.e., iL,, = r, ). With these definitions, the relative frequenciesof adsorption and desorption events in the Monte Carlo simulation will satisfy the detailed-balancecriterion [ Eq. ( 14) 1 with 0, + 8, (uide infra). Furthermore, with this choice of transition probabilities, the success-to-trial ratio will be optimized for an algorithm such as the one depicted in Fig. 1. It should be noted, however, that while this particular algorithm is reasonably effective if the time scales of various processes the system are similar, its efficiency dein clines as the stiffnessof the system increases.In stiff systems, where the majority of events are generally confined to a minority of sites, many trials will have to be attempted before a successfulevent is selected.Other, more efficient algorithms are availableL~20~2 should be utilized to simulate these and systems.Regardlessof the implementation, the relative frequencies with which various events are performed must comply with the detailed-balancecondition for both the dynamics and the equilibrium of the physical system if any meaningful comparison of the simulation of a physical system is intended. A second criterion which must be satisfied to simulate the Poisson process is proper correspondence of Monte Carlo time to real time. To this end, it should be noted that the Poisson process is, in actuality, a continuous-time version of the discrete Bernoulli processwhich is simulated by the Monte Carlo algorithm when time is measuredin terms of trials, By replacing the discrete interevent times with appropriate continuous values, the Monte Carlo algorithm produces a chain of events which is a Poisson process.The continuous interevent times are constructed through the developments leading to Eqs. (7)-( 11). Upon each trial i at which an adsorption or desorption event is realized, time should be advanced with an increment ri selectedfrom an exponential distribution with parameter ri = (NMi)rA f Afir,. (20)
Here, Mi is the number of sites occupied at trial i (i.e., ei = Mi /N) . The selection of a time increment in this way yields consistencywith Eq. (7)) which provides the distribution of interevent times for the Poisson process.Over many successfultrials at steady state, the average time between successive events is
(to)=q (N-&f,;I
A
+&f.r D . I
(21)
Here, h is the fraction of successfultrials at which Mi sites
J. Chem. Phys., Vol. 95, No. 2,15 July or copyright, see http://jcp.aip.org/jcp/copyright.jsp Downloaded 23 Nov 2005 to 128.125.4.122. Redistribution subject to AIP license 1991
1094
K. A. Fichthorn and W. H. Weinberg: Dynamical Monte Carlo simulations
are occupied. Equation (2 1) representsa time weighting of various configurations of the simulated N-site ensemble which is consistent with that dictated by the detailed-balance criterion for the equilibrium ensemble. finite-size latA tice can approximate the continuum ensembleto the extent that its sizeallows resolution of the ensemble.In general,the simulated ensembleat a particular point in time [characterized, in our simple example, by t9,(t) ] could fluctuate about the true continuum ensembleat that time [8(t), in our example] without ever achieving exactly the continuum value. Thus, time-weighted averages (which correspond to thermodynamic averagesat equilibrium if the detailed-balance criterion for thermal equilbrium is fulfilled) must be computed to estimate the true continuum ensembleat any point in time to within the desired degreeof accuracy. When time is incremented according to the procedure delineatedabove and transition probabilities are chosen to satisfy the detailed-balancecriterion, the correct macroscopic rates of adsorption and desorption can be measuredfrom the simulation at steady state. Let us consider, e.g., the rate of adsorption. Over Ssuccessfultrials resulting in either adsorption or desorption, (N,,s) adsorption events will have occurred, on the average,where
W4s) = [ w (yy ;;Iy,B, A -s = [ T Wfl? r)J;lxei]
]s s*
(22)
The correspondingamount of time which has passedduring the S successes is At
Defming the rate of adsorption as the number of adsorption eventsoccurring per site per unit time and utilizing the definition of IV,, and W, in Eq. ( 19) [or any definition satisfying Eq. ( 18) 1, the simulated rate of adsorption (R, ) becomes (K4) =pu i where
Af;:/[r" (N4i = ~jif;./[r"(N-il!f,)
i IA(N-Mj) +rDi 1
A S.
a lattice sufficiently large to representthe full systemensemble and to simulate independent events. Nevertheless, if transition probabilities are chosento satisfy the detailed-balancecriterion and if time is incremented in a procedureanalogous to that outlined in Eqs. (20) and (2 1) in thesesituations, the time seriescan be interpreted in terms of a Poisson process of one of the events (analogous to the continuum rate-limiting step approximation), and accurate ensemble averages be obtained through time averaging.For examcan ple, in the single-site adsorptiondesorption processintroduced previously, successiveadsorptions and desorptions are correlated. Nevertheless,a Poisson processcan be constructed consisting of one of the events,e.g., adsorption occurring with a rate r, ( 1 - 8). The criterion of a random selectionof lattice sitesfor potential eventspreventscorrelations from developing among specific sites and is usually achievedby utilization of an adequaterandom number generator. In certain applications, the intersite correlations induced by an inadequaterandom number generator may lead to erroneous results. The selection of an appropriate random number generator is, therefore, an issuerequiring careful consideration. Thus, through the example of the lattice gas, we have outlined the basicelementscomprising a formalism through which Monte Carlo simulations may be utilized to simulate dynamical phenomenawithin the context of the lattice-gas model. We have utilized this methodology to simulate the adsorptiondesorption algorithm of Fig. 1. Simulations were run on 128X 128 square lattices with r, = 1.0 (site s) - and r, = 2.0 (site s) - I. Transition probabilities were defined by normalization of each rate by the rate of desorption (i.e., W, = l/2, IV, = 1.0). Figure 2 depicts the fractional surface coverageof adsorbateas a function of time for an initially empty surfacefor both the transient analytical (exact) solution of Eq. ( 13) and the Monte Carlo simulation. The Monte Carlo curve is the result of one run
-eiv,
=r,(l---8,),
(24)
Mi)
+
r&f,] =+rDMj]
Ari
(25)
/ 1.2 I , I
(to>
provides the time weighting of eachconfiguration characterized, in this simple system, by 19~. similar result can be A obtained for the desorption events. A final and perhapsmore subtle criterion which must be fulfilled for the Poissonprocessto be simulated effectively is that independence eventscomprising the time sequence of of the process is achieved. Strictly speaking, the formalism which we have presented is valid only when independent eventsare simulated. In general,the independence succesof sive trials is ensuredboth by utilization of a system which is sufficiently large that single site and intersite correlations are lost, and by random selection of sites on which potential eventsmay occur. The former condition is not always possible to achieve,particularly if the system is stiff. In such systems, it may not always be feasiblecomputationally to define
0.0
0.5
1 .o t s
1:5
2.0
2.5
FIG. 2. The transient solution of Eq. (12) for the initial condition of an empty surface [ 0( t = 0) = 0] provided by both the analytical form of Eq. (13) and the MonteCarlo algorithm of Fig. 1 with r, = 1.0 (sites)- and r,, = 2.0 (site s) - ( WA = l/2 and W, = 1). The inset depicts the rate of adsorption measured from the Monte Carlo simulation at steady-state and the continuum steady-state rate.
J. Chem. Phys., Vol. 95, No. 2,15 July 1991
Downloaded 23 Nov 2005 to 128.125.4.122. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
K. A. Fichthorn and W. H. Weinberg: Dynamical Monte Carlo simulations
1095
only. There is excellent agreementbetweenthe two solutions in the approach to equilibrium and the steady-state fractional surface coverage.The inset of Fig. 2 depicts both the analytical (exact) and Monte Carlo steady-staterate of adsorption (desorption) . The rate of adsorption on a per site basiswas measuredas the reciprocal of the variable amount of time At required for 50 adsorption events (i.e., r, = 50/ NAt, and N is the size of the lattice). It can be seenthat this rate fluctuates about the predicted continuum rate of r, = 0.666...(site s) - Of course, the amplitude of the fluctu. ation dependson the number of eventsaveraged.The mean rate measured from the simulation over the depicted time interval of 2.0 s was 0.67 f 0.06 adsorptions/site/s.
be partitioned among the various possibletransition events asN={n 1, 2,...,nk1, where ni is the number of speciescapan ble of undergoing a transition with a rate ri and i%l Thus, a particular configuration of the systemat a particular time can be characterized by the distribution of N over R. This distribution is constructed by a Monte Carlo algorithm which selectsrandomly among various possibleeventsavailable at each time step and which effects the events with appropriate transition probabilities W = {w1,w2,...,wk1. The transition probabilities should be constructed in terms of R so that detailed balance is achieved at thermal equilibrium and a dynamical hierarchy, as expressedin Eq. ( 18)) of transition rates is preservedaway from equilibrium. If a sufficiently large system is utilized to assurethat the independence of various events is achieved, then the Monte Carlo algorithm effectively simulates the Poissonprocess,and the passage real time can be maintained in terms of R and N. of To accomplish this, at each trial i at which an event is real...