B. Rabinovitch
University of Oklahoma Medical School Oklahoma City, 73104
The Monte Carlo Method Plotting the course o f complex reactions
I n the course of reading in reaction kinetics, a very useful tool has come to our attention which can help enormously in solving certain problems in mechanisms of organic and biological reactions ( 1 ) . Moreover, it seems to us that the tool can be readily adapted to many other areas including mechanisms of inorganic reactions, catalysis, surface chemistry, and even such physical processes as rates of evaporation and condensation, rates of adsorption, and growth phenomena. The method is far from being new but up till now, we believe, it has been of restricted value, both as a research tool and for educational purposes, because it has required the usc of a digital computer. The present paper describes a technique that adapts this method for use a t thc desk with nothing more than pencil, paper, and a table of random numbers. I n addition, its use for educational purposes will be very clear, presenting as it does a model for the reacting system very easily grasped by the student. The technique will provide data describing the course of any complex reaction and the rate of appearance of any or all of the reactants or products. I t will not solve differential equations describing the reaction nor supply closed form rate equations. The technique involves the consideration of the statistics of random events and the probability of their occurring in concert, simultaneous!y, consecutively, or repeatedly. This technique is often referred to as the Monte Carlo (2) approach and may he applied quite rigorously to the problem of the probability that a given molecule will spontaneously decompose in a monomolecular reaction process, or that t.wo molecules of the same or different types will react together in a bimolccular reaction process. The Monte Carlo Method
I n this technique, digits play the role of molecules, molecules of the same kind by the same digit, of different kinds by different digits. The digits are placed in a matrix of some sort and their positions coded in some way. The matrix is then pin-pointed in a completely random fashion and repeatedly. When a given pinpointed position is already occupied by a digit representing an urireacted molecule, it is capable of reacting, and this is a L'fruitful" selection. If the process is a monomolecular one, the given molecule represented by the selected digit will react, and the digit will be changed to somc other digit representing some other molecule. If the process is a bimolecular one, the selection of a first digit must be accompanied by the simultaneous sclcctio~~ of a sccond digit of thc appropriate kind reprcsel~thga molecule which can react with the first, i n order to hc fruitful. Both digits are then changed to 262 / Journal of Chemical Education
some other digits representing reaction products. I n the monomolecular case, if the random selection of a position in the matrix yields a digit representing a reacted molecule, or a molecule of the wrong kind for the process we are considering, or a reaction product, or even the absence of a molecule, then such a selection will represent an unfruitful event, or simply the passage of an element of time without the appropriate reaction. I n the bimolecular case, unless both random selections of a position in the matrix yield digits representing unreacted molecules of the appropriate kind, then no elementary reaction process will take place. The method, when used with a digital computer, presents relatively few problems, since each mcmory bank used as a matrix has a capacity for storage of digits of many thousands. Moreover, it may generate its own random numbers one a t a time and as many times as is consistent with good statistical sampling and the desired accuracy. With the manual technique, used a t a desk, the matrix used will, in general, have a much lower capacity for storage and therefore the statistics will be subject to a much larger error. I n spite of the fact that the accuracy in computing the course of a reaction is lower by the manual technique to be described here, it can be improved by successive repetitions of the "experiment" and subsequent averaging of the results. Method of Computation
I n this manual technique, the matrix is simply a square grid representing the reaction vessel, with a predetermined number of positions each of which is numbered in sequence (see Fig. 1). For illustrative purposes, each position is conceived of as representing a molecule so that we may place on this grid, in a random fashion, a chosen number 0 I 2 3 4 5 6 7 8 9 of, say, digits 1representing o a given reactant of a pre10 selected concentration. The 20 remainder of the positions 30 may be considered to be 40 filled with solvent moleso cules, accounting for a mole 60 fraction of 1.0. The dism position of the digits 1 need 80 not be random, since the 90 elementary acts of selection of sites leading to reaction Figure 1. Grid matrix of 100 are themselves random, positions representing o reaction vessel containing 20 mole yo of nevertheless such random digits 1 and on equimolar soncen- disposition simulates the tratian of digits 2. The readion is bimoleculor ond shown to b e 25y0 randomdeploymentof molecules in solution. complete.
Thus, if all positions on the grid are occupied with reactant, then we can consider that we are dealing with a pure substance. If a second reactant is involved, then either a second grid of proper size is used with the appropriate number of digits placed arbitrarily on it, or the appropriate number of digits 2 may be placed on the first grid t,o give the preselected conccntration of the second reactant. Again, if bot,h reactants are pure substances, i.e., not in solution, then all positions on both grids will be occupied by digits. I n the examples that we will represent below, we have used a grid of 100 positions, i.e., no solvent. If a system in solution is being considered, a larger grid should be used, say 1000 or 10,000 positions. The positions corresponding to the set of random numbers of 3 or 4 digits, dependiug on the size of the grid used, are marked with, say, a 1 unt,il the number of marked positions corresponds to t,he concentration required. Thus, 100 marked positions in a grid of 10,000 corresponds to a solution of 1 mole pcrccnt. To revert to the conditions we have used, namely pure reactants, a conveniently small "set" (5 or 10) of 2 digit random numbers (3) is marked off on the grid, and at the end of each set a tally is made of the total numher of fruitful selections (unit reactions) against the total number of "tries" (N). The set of tries, then, represents a unit of time, and the quantity N the total elapsed time, since, in a monomolecular reaction, such as radioactive decay, it corresponds to the number of "attempts" of the atom to spontaneously disintegrate, while in a bimolecular reaction, it corresponds to the number of collisions, both of which processes are going on at a constant rate, a t least, statistically, under fixed conditions. If the reaction being considered is bimolecular, a conveniently small set of pairs of two digit numbers is read from the table, the first of the pair applied to one grid, the second to the other grid. If the reaction scheme involves the further reaction of a product of a previous step (i.e., in a consecutive reaction), then the marking off on a given grid of a fruitful selection in the previous step is tantamount to the introduction of a new digit 3 which can become involved in the consecutive step. If two react,ionsare going on simultaneously at different rates, the liberty is taken of giving each branch a number of tries in order, using small sets of numbers, each branch receiving the proper number of tries in relation to the other branch, as described below. The smaller the size of the set of numbers used, the closer to simultaneity one gets, but the greater the inconvenience. I n this way, a computation of the course of a reaction may be made in terms of elapsed time represented by N = k.t/y where k is the reaction rate constant; t, time; and 7 , a constant, being the Monte Carlo equivalent to a rate constant, in terms of units of concentration of reactants of [I] = a[A], where [A] is the concentration of reactant in moles/l; [I], the ratio of number of digits 1 to total number of positions in the matrix representing the reaction vessel; and a, a proportionality constant obtained as follows. According to the definition above, [ I ] corresponds to a mole fraction, since each position on the grid corresponds to a molecule. Thus, it is easily shown that for a molar concentration [A],, the mole fraction of A =
+
Plo = [AIoM./([AloM, - [AIoMA 1000d), where MA and M , are the molecular weights of A and solvent respectively, and d is the density of the solution. Therefore, we can say or = M,/([AIoM, - [A],MA 1000 cl). Such a procedure may be extended to reactions of any order, to consecutive reactions, and to simultaneous competitive reactions such that any reaction scheme of no matter what complexity may be written down in these terms and the rates of the various steps or branches computed. One only needs to relate the ratio of the number of tries in each branch or step directly to the ratio of reaction rate constants for these steps. Unfortunately, these two ratios are equal only if all the branches and steps in the scheme are of the same order. If not, then i t is quite easy to determine the ratio of tries, as follows. Let us take an arbitrary reaction scheme such as the following
+
kl
A-B+C k.
A+B-D+C
*I
B-E+F B,
B+E-G+F
step one step two step three step four
and let the digits 1, 2, and 3 represent the species A, B, and E, respectively. Then we may write that the digits 1 will disappear from the matrix according to the following expression
where Nl is the number of tries applied to step one; dN,, a small increment in this number of tries; and P, a proportionality constant which is a property of random statistics and independent of the order of the reaction. On integration this equation gives
where [lIo is the value of [I] before any tries have been made, i.e., when N1 = 0. But since 111 = a [A], we can write that
which immediately identifies ON1
-- kd
Similarly, for the rate of disappearance of digits 2 from the matrix by step three, we can identify ONBN3 = k3twhere N3 is the number of tries applied to step three. This tells us, as already stated, that for these two firstorder steps NI
kt
N,=k; For step two, however, we will obtain
which describes the disappearance of digits 1 and 2 from the matrix by the process of step two, a bimolecVolume 46, Number 5, May 1969
/
263
ular process. This expression gives
which identifies &N% = kd
Similarly, for the disappearance of digits 2 and 3 from the matrix by the second-order process of step four, we will identify mpN4 = k4t, which again demonstrates for steps of the same kinetic order that
%=!? N,
k,
However, i t is now clear that the ratio of tries for steps of different kinetic order will no longer be simply the ratio of rate constants for these steps. Thus,
Results of Computations
Ordinarily, to obtain an accuracy of better than 1% in data describing the course of reaction, the number of digits representing molecules that we need will be 10,000 or more. The number of tries required to give say 60% reaction will be in the order of several tens of thousands, and to accomplish 90% reaction will perhaps be as much as 100,000. This represents no prohlem if a digital computer is used, provided that it has a large enough memory for storage of random numbers. However, for many purposes an accuracy, in describing the course of the reaction, of perhaps 5% is perfectly acceptable to investigations of reaction mechanisms, at least for separating out possible from impossible schemes. Moreover, accuracy in this computation may be materially improved by repetition of the process and subsequent averaging of the results. This is demonstrated in some of the results given be-
low, but it is worthwhile to point out that accuracy will only improve by a factor which is the square root of the factor by which the number of digits, representing molecules, is increased. First-Order Reaction
We have taken here the case of the decomposition of a pure substance by a unimolecular process. One hundred digits were used but the reaction was repeated separately three times. Figure 2 shows the data plotted as for a first-order reaction, the first curve representing the first reaction, the second curve the mean of all three reactions. The data conform very well to firstorder kinetics, the improvement in the data in the second curve being noticeable but not marked. Second-Order Reaction
Again, two pure substances in equimolar amounts were considered to he reacting in a himolecular reaction. One hundred digits for each substance were taken, and the reaction repeated a second timc. The rate equation for such a simple reaction is
and the data as shown in Figure 3 conform very well
Figure 3 . a = 0.5.
Second-order reaction. Pure rvb+mces of equimolor omovnh, A, results of one experiment using 100 digib of each reoctont, abscissa lower scale; B, results averaged for two experiments, eoch of I00 digits, absciwo is upper rcale.
to this expression. The improvement in the data as represented in the second curve on this figure (mean of two reactions) is again noticeable but not marked. First-Order Competitive (Simultaneous) Reoction
This reaction can be represented by the scheme 11
kc
C-A-B
Figure 2. Firsborder reoction. Pure substonce, a = 1.0 A, rowlts of one experiment using 100 digits, obscisro is lower scale; B, rerulh overogsd over three experiments, each of 100 digits, abscissa is upper scale.
A pure substance A represented by 100 digits 1 was allowed to decompose simultaneously by two path ways, but monomolecular, into products B and C represented by digits 2 and 3, respectively. The ratio kl/k2 of rate constants was taken as 2, and the reaction repeated two times. The number of counts for branch 1 was taken in sets of 10, for branch 2 in sets of 5. The rate equations representing this reaction can readily be shown to be
the reaction scheme. However, if one cannot obtain a useful rate equation for the reaction scheme, one is often forced t,o make simplifying assumptions which are either questionable or difficult to validate. This shows the clear advantage of the present technique being described, where one can readily plot the course of a reaction without proceeding through a rate equation. By way of example, the presently considered scheme of a first-order consecutive reaction gives rise to the following rate equations
and
and
Figure 4 shows the data obtained as a mean of the two reactions, plotted appropriately. The time scale for product B is Nz while for the disappearance of A is NI N2, and this is simply to separate the two curves. If the time scale for B had been N , N2,the two curves would coincide very well as required by the above rat,e equations. The conformity to these equations is very good.
+
+
Figure 4. First-order competitive reaction. Pure wbrtonce, a = 1.0, k d k * = 2.0. A, rewlts for appearance of digitr 2, ordinate 1 - ( 3 . 1 2 ] / 2 ~ [ l ] a ) ; B, rervitr for diroppeorance of digitr I . ordinote 1.
Although quite a simple reaction scheme, the derived rate equations cannot be easily aud rigorously used to derive k2. If we proceed a little further in complexity to a simple competitive second-order scheme--this will be examined in thc next section-we already are denied the nse of a closed form rate equation. A reaction scheme such as the following
gives rise, with great difficulty, to such a complex rate cquat,ion that I hesitate to reproduce it here, containing as it docs Bcssel fuuctions and Hankel functions of both the zeroth- and first-order ( I ) . I n the computation of the present reaction, we have considered a pure substance A, represented by 100 digits 1, dccomposiog uuimolecularly and successively to B and C, represented by digits 2 and 3. The ratio of rate constants lcl/k2 was taken as 4, the number of t,ries in step 1 being taken in sets of S, in step 2 in sets of 2. Figure 5 shows the results of this computation wit,h the time scale taken as N,, t,he number of tries in the first step. The curves have been calculated from the above expressions and the time scalc adjusted to give the best fit, for the case of the appearance of B, and the same scale used for the disappearance of A and the appearance of C. The fit for B is very good throughout, but thcre is clearly a deviation between
Firsf-Order Consecutive Reaction
We may represent this reaction in the followiug way I,
In
A-B-C
I t is appropriate a t this point to consider what happens to the rate equations describing a reaction as the reaction scheme becomes morc and more complex. I t very rapidly becomes extremely difficult to obtain a closed-form rate equation, if at all, and if oue is successful in doiug this, the derivation of rate constants or ratio of rate corrstauts from such equations is more often than not prohibitively difficult by a rigorous mathematical technique. One is forced to verify the validity of the scheme under consideration by showing that the experimental data for the reaction conform to the theoretical course of the reaction as derived by the rate equation. The more reactants and products one call follow in this way, the stronger is the support for
Figure 5. First-order consecutive reoction. Pure rub+mce, a = 1.0, k n / k 9 = 4.0. Curves colcvlated fmm theory. A, rervltr for dirappearonce of digitr 1 , ordinote [I]/[la];B, rervitr for appearance of intermediate digits 2, ordinote 12]/[l]o; C, rewltr for appearance of product digitr 3, ordinate [3]/[lIa.
Volume 46, Number 5, May 7 969
/
265
calculated and computed in the last parts of the A and C curves. These may well disappear if the computation is repeated several times, or these deviations might have arisen because of the finite intervals for N that are involved. The agreement, however, is considered a substantial confirmation of the technique. Second-Order Competitive, Consecutive Reaction
This reaction, although remaining relatively simple in its reaction scheme, already gives differential equations that have not been solved, except under special circumstances (4). The reaction may he represented as follows
pure. However, this should make no difference to the results since both branches are second order. The results of these computations are shown in Figure 6 where the points represent the computed course of the reaction for the first reaction and the mean of the two reactions. The curves are the courses of the appearance and disappearance of products and reactants calculated from the equations above and the time scale adjusted to give the best fit to the computed data for the disappearance of A. No further adjustment was made for curves B and C. The curves are seen to fit the computed data very well, with a noticeable improvement for the mean data. Pyrolysis of a Hydrocarbon
For the most general case, where the ratio k l / k 2 is an arbitrary value and the concentrations [Ale and [BIo, the starting materials, also have no simple relationship to each other, the problem of obtaining a closed-form rate equation remains unsolved. However, we have chosen the special case where [Ale = 2[BIo and kl/lcz = 2. Under these circumstances it is quite easy to show that
We have taken A and B as substances represented by 100 digits 1 and 50 digits 2, respectively. The reaction was carried out twice with the number of tries in branch 1 taken in sets of 8, in branch 2 in sets of 4. I n the first computation, the 50 digits 2 were placed randomly on a grid of 100, whereas in the second computation they filled a grid of 50. This means that in the first computation the reactants are not pure hut have mole fractions of 0.5 and 0.25, respectively, whereas in the second computation the reactants are
By way of example of an experimental procedure, we have taken the case of the pyrolysis of neopentane. The proposed reaction scheme is (5)
Usually, in such a scheme, the stationary state .approximation is applied to the concentrations of all intermediate species, i.e., the free radical concentrations are considered to be constant and very low. As a first step we have taken pure neopentane, represented by 100 digits 1 and taken different relative rate constants for the separate branches and steps. Moreover, insofar as we know experimentally that very little ethane is produced, a t least in the early part of the reaction, we have ignored branch 5. Under these conditions it is instructive to see how the relative amounts of the various products and the intermediates vary with differing relative rate constants, and to determine under what conditions the steady state approximation is valid. We have taken no pains to pursue this problem to its ultimate and proper conclusion, since this is not our purpose. In fact, where relative rates for the branches of the above scheme are small, the kinetic chain length of the chain reaction will he small, and properly we should include six more branches representing the various recombination reactions of the three intermediate free radicals. We have avoided this since we use the system only as an example, but the added work is by no means prohibitive. Representing the reaction scheme as
"i, 10
u
Do
w
1b2 224 256 8 N1 Figure 6. Second order competitive consecutive reoction, [I]" = 2[210,k l l k s = 2.0. Open circles ore resulh far 100 digih 1; solid circles are results averaged for 2 reactions. Curve calculated from lheory. A, curve for the disappearonce of digit, I , ordinate [ I ] ; B , w n e for lhe dimppeorance of digits 2, ordinate [2]; C, curve for the appearonce of digits 3, ordinate[3].
266
/
32
M
118 1k0
Journal o f Chemical Education
we have taken A as a pure substance, represented by 100 digits 1 and set up 4 grids representing the amounts of A, B, C, and E throughout the reaction. Reaction in branch one gives rise to digits in grids 2 and 3;
Pyrolyrir of o hydrocarbon. Pure wbrtonce, a = 1.0,k1 = Acc~mulotionof stable products. Curve A, rervlh for the disoppeoronce of digits I, ordinate [ I ] ; B, results for the appeoronce of digits 4, ordinate 141; C, results for the oppcorance of digits 6, ordinate 161; D, results for the appearance of digit. 7,ordin.k 171; E, rewlk for the change in pressure in a completely gaseous system, ordinote APIP;.
Figure 7.
kr
=
As for Figvre 7 with
Figure 9.
ki:k2:ks:kn = 1 : 4 : 4 : 4 .
ka = k,.
reactions in branches two and three give rise to digits appearing in grid 4; and reaction in step four gives rise to further digits appearing in grid 2. I n this way, we were able to compute the change in amounts of the various species for differellt relative rates, and the results are shown in Figures 7-12. The smoothness in the data is quite good for case 1 (Figs. 7 and 8) hut gets p~ogwssively worse from case 2 to case 3. The reason for this lies in the disparity and absolute size of the number of tries per set for each branch or step, and this is dependent upon the relative rate constants. However, the accuracy of the data can be improved by repetition of thc experiment.
It can he seen that the course of the appearance and disappearance, and thc relative amounts, of reactant and products can be changed radically by controlling relative rates. We have included in Figures 7, 9, and 11 the change in pressure A P / P 1 where P, is the initial pressure of the system, assuming a completely gaseous reaction. It is apparent that by so changing the relative rate contants for the different branches and steps and comparison with experimental data, enough information is available from such computations to draw conclusions as to the mechanism of the reaction. Conclusions We have show11 that the computation of the course of a reaction by a manual Monte Carlo technique will conform to known rate equations for given reaction paths. The accuracy of such computations can be governed by the amount of energy one is prepared to put into them and can be made as good as experimental
-tI
m
S
eE
' 10 * 0
U
:L 10
Figvre 8. Some reaction or in Figure 7 ; occvmvlotion of intermedioter. A, rerultr for appearance of digits 2, ordinote 121; B, r e d t r for appearance of digits 3, ordinote 131; C, rewlh for appearonce of digits 5, ordinate 151.
0
Figure 10.
NI
As for Figure 8 with
k1:k9.kr.kr = 1 .4:4:4.
Volume 46, Number 5, Moy 1969
/
267
Figure 11.
As for Figure 7 with
h:h:ka:kr
= 1 :10:10:4.
data usually obtained. Application of such a technique to a complex reaction schemc can, by comparison with experimental data, afford much useful iuformat,ion, and, in the final analysis, may demonstrate the correctness of a chosen mechanism and evaluate rate constants. Other properties of the reacting system can be computed as wcll, as shown by the change it1 pressure of a gaseous system. The dependencc of, say, limiting conceutrations of a product, or maximyn concent~ration of an intermediate, upon concentration of the reactants could readily be dctermined. Application of thc technique to problems such as catalysis and rates of adsorption can readily be made by consideration of the solid surface as being made up of reactive sites of varying concentration, varying reactivity. The problem of the particle size distribution of alumina in a combustion react,ion of aluminum should also be amenable to attack by this technique, by consideration of the various possible reaction pathways for the oxidation of aluminum coupled together with the phenomena of nucleation, condensation, and evaporat,ion. I n fact,
wherever a chance phenomenon is involved, such as in the collision of molecules or atoms with each other or with a solid surface, the arrangement of molecules or atoms in a physical system (entropy), and where this can be related to the probability of selection of locations in a matrix, then useful results can be obtained by this technique. Literature Cifed
SCHAAD, L. J., J . Am. Chem. Sac., 85, 3.588 (1963). Housa~o~o~sn, A. S., FORSYTHII, G . E., A N D GRBMOND, H. H., "Monte Carlo Met,hod," National Bureau of Standards, Applied Mathematies Series 12, Washingt,on,1). C., 1951. (3) ARKIN, H., AND COI.TON, R. R., "Tables for St,atist,icians," Colleee Outline Series 7.5, Barn- and Noble, New York, 1950,~.142. (4) FROST, A. A,, A N D PEARSON, 13. G., "Kinetics and Mechanism," John Wiley & Sons, New York, 1961. (5) MINKOFF, G. J., A N D TIP PI:^, C.F. H., "Chemistry of Combustion Reactions," Bobterworth, London, 1962.
(1) (2)
+
268
/
lournol of Chemicol Education
