Interactive Program System for Integration of Reaction Rate Equations John P. Cheslck Haverford College Haverford. PA 19041
Intemation of the reaction rate equations for multistep reactiin mechanisms has been a necessary step in research dealing with the reaction kinetics of complicated systems. For examole. oredictions of the effects of chlorofluorocarbons on the stratospheric ozone system are based on such numerical inteeration of the rate eauations for svstems of " many elementary reactions whose rate constants are more or less well known. Integration of the rate laws for multistep photochemical smog reaction systems provides much of the basis for understanding and strategies for controlling photochemical smog in the troposphere. Transport processes complicate these particular applications, but the exact integration of the couoled rate laws involvine manv reactions and many species is central to these simulations. In general, one need not include verv manv stens in a reaction mechanism before intuition pro;ides l i k e ielp in predicting chemical behavior or understandine the sensitivitv of re dictions to uncertainties in rate cons& values for individual elementary reactions. On the other hand, pedagogical treatments of chemical kinetics do not usually get much beyond discussion of the integrated forms of results for simple first- and second-order reactions, and the formulation of the coupled differential rate equations for multistep reactions. The use of the steady state approximation is introduced to eliminate, where algebra s i m ~ l vnermits. the concentrations of reactive intermediates fiomthese rate equations. In favorable cases a simple result is obtained for the final "rate law". Mathematical complexity prevents exact treatment of significantly complex systems, and intuition becomes a weaker support. I t is difficult to demonstrate the conditions for validity of various approximations without having available simulation software that is both flexible andeasy touse. "Numerical experimentation" with simulations can help to provide understanding of various approximations and their limitations, and to illustrate the role of simulations in research and practice. It is relatively easy to program the rate equations for one mechanism, to be integrated by some standard numerical analysis routine, but practical use requires a simulation system that will permit the user easily to modify the e effects of changes of reaction scheme in order to e x ~ l o r the mechanism on predicted kinetic behavior.. Each instructor will also have uniaue needs and desires for examples and time for their demonstration and/or use by students, and the special-case program will not meet these general needs. Edelson some time ago (I)provided a lucid and strong statement concerning the value of incorporating such simulations in chemical education a t the undergraduate level. We describe here our Pascal-language kinetics rate law integration package for t h e desktop microcomputer. HAVCHM allows completely general choice of reaction mechanism. UD to the dimension limits of the Dromam, and thus makespossible versatile simulations of kinetic behavior resultine from the inout choices for elementary reaction set and numerical values for rate constants and initial concentrations. No programming is required for alterations in mechanism. This implementation was designed for desk top microcomputer usage; any computer system that can operate Pascal with sufficient precision should be able to utilize this system. The precision issue will be addressed later.
-
The following examples are to provide a few illustrations of use and sueeestions oointine to the flexibility of simulation choices. f & ~ r e s1,i,and 3ihow results obtained for the reaction mechanism
F l m 1. CMIcenlratIm~. time fwswcles in mechanism 1, all rate constants
0
10
20
Figue 2. Concentrstionvs. time for specles in mechanism 1, k, = k~ = kd10.
Figure 3. Induction period for system in Figure 2,
Volume 65 Number 7 July 1988
599
Figure 1 shows the results for k, = kn = ka. No approximations hold, the steady state assumption is not even approximately correct, and the exact integrated rate law for even this relatively simple case is not found in most physical chemistry texts. The exact solution might he presented to students for comparison with numerical results from the numerical integration. Figure 2 shows the results for kl = ka = kzllO, with all rate constants scaled to give about the same half life for Aas in Figure 1.The steady state treatment for B is reasonably valid, as the maximum in that concentration is reached in a time of less than 2% of the half life for A. Figure 3 is an expansion of the calculated results picked to show that the induction period, during which the steady state concentration of B is established, is obtained as part of the calculation process. The A concentration is off scale for this graph. Behavior during this induction period is accurately calculated by the integration package even when its duration is less than of the reaction half-life. The user may compare these calculational results to those from the steady state treatment, explore variations in the rate constants to determine the requirements for the steady state treatment
to be valid, and then predict the behavior of the system if, for example, the first reaction is made bimolecular rather than unimolecular in A. A change such as this is readily made in the input file and the Droeram operated to show those results. he relation of r k e c o n s t a k s for satisfaction of the steady state approximation turns out to be different for a bimolecular first reaction. This is seen from operation of the simulation, and it is not simply predictable in advance. Other changes in the reaction mechanism are equally readily made for comparison. Parallel reactions may he added leading to different products. Reaction steps of the form may he added to the above reaction scheme, leading to description of chain reactions, most simply represented in the .~ . following mechanism:
by A as an initiator, C as a reagent, P as a product, and B a s a chain carrier species. The relation between reaction rate and rate constants for the initiation, termination, and propagation steps may be explored. Conditions for significance of other first-order termination steps may be modeled. Branchine reactions mav be added with imnlications for describing flame and explosion phenomena. Figures 4 and 5 illustrate the Oregonator (2) reaction mechanism, one of the simplest reaction sets that will provide an example of an oscillating chemical reaction. I t also produces some transients in each cycle that are truly amazing for their sharpness and range and that provide a strong test for any rate law integration package. The mechanism is A+YY-X X+YY-P X+B-X+X+Z x+x-Q Z-YY
30
100
170
Figure 4. log (canc) M. time for Onrgonator (2)oscillatory reaction mechanism.
Here P and Q are products that do not figure further in the reaction scheme and therefore mav be taken as fixed or constant in considering the hehavio; of transients. A and B are reagents. Fixing these concentrations provides the perpetual motion aspect of the system; the reaction never runs out of fuel! The variables in the integration are then the species X, YY, and Z. There is a range of initial concentra- , tions and rate constants for which the system will lock into oscillatorv behavior, as shown in Fieure 4. which shows three repetitions of the transients. Notice the log srale for concentrations: X changes . by . a factor of 10"over a verv short time. Figure 5 shows an expansion of a portion of Figure 4 to illustrate the transients for X and Z. This was selected in the use of the plotting routine from the results plotted in Figure 4. One can explore the range of rate constants and initial conditions that lead tu theoscillatory bshat,ior, and questions may be asked concerning understanding of the system. It How could one vredirr that oscillation mieht be nossihle? " ~ - - - - - -is certaintly true that a necessary test for a proposed reaction mechanism for auv real oscillatorv reaction is that the proposed scheme actually yield oscillations in a kinetic simulation. At least one published mechanism. with rate conbsci~~atory reacstant values, for the ~klousov-~habotinskii tion fails this test (3). ~
Figure 5. Detatl of Figure 4 tw lhe species X a d Z in oregonator (2)meche nism.
600
Journal of Chemical Education
~~
.~ ~
Program Description The general organization of the reaction coding and generation of the rate expression for integration was described in an earlier article (4) and will not be repeated here. Some words on the problems of integration of chemical rate equations are, however, in order. It has been previously discussed
(4,5) that the variety of time constants found in most such multisten svstems reauires variable sten inteeration techniques, df which the i ear algorithm (6,ijhas been a foundation. Rate law integration packages that are set up to allow flexihle input of multistep mechanisms must include the Gear or equivalent aleorithm to allow for the "stiffness" prohlem, tl;e mathemaiical difficulty arising from the grossly different time scales for the change of different \.ariahles in the system. Applicability of the iteady state hypothesis for some species virtually guarantees that this stiffness prohlem will arise in numerical integration of the rate equations. The Gear integration algorithm was encoded in Pascal and is huilt into the HAVCHM program system, and special provision need not he taken to cope with reaction schemes showing such sharp transient behavior for some species. Special purpose rate law integration packages that are designed for simpler cases miss interesting and useful examples. Our first version of this HAVCHM program system (4) was unique in its incorporation of the Gear integration algorithm into a kinetics simulation package that provided ease and flexihilitv in snecifvine the set of elementarv reactions and initial conditions dth-automatic generationof the rate eauation svstem. This . nroeram svstem was in FORTRAN and was designed for batch operation on a mainframe computer. In a subsequent note (8)we indicated the availability of an interactive version, also in FORTRAN, for use on a DEC 11-03 minicomputer or other terminal system which could run FORTRAN. However, FORTRAN is not currently the programming language of choice for the microcomputer, and the "interactive" version previously described was primitive in its editing and parameter storage capabilities. Field descrihed (5) a BASIC translation of a simplified version of the Gear integration algorithm, hut did not include a general program system to permit flexihle variation of reaction mechanism without reprogramming. Integration of a set of rate laws from a chosen mechanism with the HAVCHM program and the plotting of the results were deliberately separated into two companion programs. The companion plotting program CMPLO can then be run independently to provide quick illustration of the stored integration results ohtained from a variety of runs using HAVCHM with different initial conditions, different choice of rate constants, or with changes in the reaction mechanism.
" .
Program Operation The dimensioning is currently set to allow up to 30 elementary reactions with a total of up to 20 different species appearing in these reactions. A? many as three reactant molecules and also three product molecules may he specified for each reaction. Reactions of greater complexity can scarcely he termed "elementary". A reaction is defined by the input of brief character strings as labels for the reaction and for eachreactant and eachproduct species. Thereactant andlor product label is repeated for a reaction that is of second- or third-order in one molecule and/or that nroduces more than one molecule of a given species as a product. Additional inputs are the rate constant values for the reactions and the values for any non-zero initial concentrations. Anv of the s ~ e c i emav s he held constant at the input values. ~ n unigmay y he used for time, concentration& and rate constants, as long as all are consistent. The numher and identity of the reactant species define the molecularity and the rate law for each reaction with automatic generation of the total rate of formation for each snecies. Stoichiometric balances are automatically satisfied. The user is given the chance to approve and edit for changes or corrections the input reaction and initial condition parameter set. The reaction set, initial conditions, and set of operating parameters specifying output options and integration time are saved in a file. Subsequent program operation may use this file unchanged or edit it to minimize
input if the same hasic mechanism is to he run more than once. A variety of these parameter files may then he created if it is desired to have the user run through a numher of examples without taking the trouble to set up the detailed input for each case. Output Form The computational results may come in two forms; either or both may be output, with optional printing, and one form is stored as a file for plotting by the program CMPLO. Stepwise output provides the output of time and all variable concentrations whenever a specified numher of integration steps has been taken. These integration time steps are of variahle size, depending on the rapidity of change of transients and a parameter that controls the tightness of the operation of the Gear integration routine. This "equal-stepnumher" output will therefore come at varying time intervals. This provides resolution of detail in the time neighborhood of rapidly changing concentrations, for example during an induction period in which a steady-state concentration is being established. Larger time steps are taken during later stages of the reaction when changes are relatively slow. The other possible output comes at equal time intervals, ohtained by the program through an interpolation procedure. These equal-time-interval results would then he equally spaced on the time axis. CMPLO Plottlng Program This program makes plots from the HAVCM output file. From one to three species may he selected for plotting on any one graph. A linear concentration vs. time plot may be made, or alog concentration vs. time plot may he chosen. The latter form is particularly illustrative in comparing species whose concentrations differ widelv. The file is scanned for the species selected to determine the limits of time and concentration. and either (1) these full limits are used in nlottine the results for thosk species, or (2) limits for time and/& concentration or log (concentration) are selected by the user for plotting. The user is given opportunity to verify and change these choices before proceding. Printing of the plot is also optional. Thus the general scope of the variation in time of an intermediate may he illustrated or the initial transient showing the formation of this intermediate may he expanded on the time axis for study. For these comparisons HAVCHM stepwise output is to be preferred to equal-time output. lmplementatlon The system descrihed here was built using the MS-DOS operating system with Turho Pascal, version 3.0. It is recommended that the microcomputers he equipped for operation with the 8087 math coprocessing chip. The chip, readily available for $75-150, both greatly speeds numerical work and also provides the double-precision word length needed for most comfortable use of the Gear numerical analysis routine. The Pascal source code can also he compiled and run without the 8087 processor using the Turho Pascal system a ~ ~ r o n r i afor t e the standard shorter word. The Droeram . .. s h o u l d a l s ~be translarahle to a similar, structured language such as C without too much difficulty. The Turho Pawal text requires only minor modifications, chiefly in file record 110, for use as the source for other Pascal systems. Full utilization of the power of the Gear integration method requires more than a minimum-word-length microcomputer. Integrations that must include sharp transients, such as the Oregonator cyclic reaction scheme, or very short induction periods, fail when run with the shorter single precision word length of the Apple I1 system for example, although they are generally successful with the standard IBM/ AT&T word length. Such failure shows in the form of either missed sharp transients or underflow errors, depending on the setting of an operation parameter for the integration Volume 65
Number 7 July 1986
601
routine. The Apple Pascal version of this program system worked with an Apple 11+computer for many examples, hut required more selectivity in its usage t o avoid floating-point underflow termination. Successful operation with the Oregonator mechanism was not found to be possible with the Apple 11+system because of the short word length. Two similar Pascal source programs are provided for plotting. One version, CMPLO, assumes the Turbo Graphix Toolkit graphics package for drawing and laheling axes and for the actual plotting of the data array with different line types for each species plotted. This CMPLO plotting program reads the data file created by HAVCHM and must also be generated to utilize the same word length as used for the HAVCHM operation. CMPLO was written to do the maximum of data reo oar at ion and establishment of limits for plotting hefur; making graphics procedure calls. This pn)duces t h e o ~ t i m u maualitypraphicsoutput. TheTurhoGra.. . phix procedure callaare closely confined in one procedure in CMPLO. Graphics routines are particularly devicelsystem dependent. For example, it is well known that the IBM PC with high resolution BAN IBM monitor requires a "foreign" card, such as the Hercules card, for graphics use. In addition, the Turbo Graphix routines are required for convenient -eraohics . use with this monitor svstem. The CMPLO Pascal source code has been compiled i i t h the Graphix procedures and is furnished in three versions as executable code for the IBM BIW-Hercules monitor, the IBM color monitor, and the AT&T 6300 cases. The second version of the plotting program, CMPLSTD, was produced by replacing the Turbo Graphix plotting procedures with the more limited set of graphics commands available in the high-resolution graphics mode of the Turbo Pascal system. No additional graphics software package is required for compilation of this Pascal source. T h e CMPLSTD plotting program will run with the IBM color graphics monitor or the AT&T 6300 system or equivalent
602
Journal of Chemical Education
and gives a more primitive plot than is available when using the Turbo Graohix oackaee in the CMPI.0 version but is in general less fussy about tge graphics terminal used. The HAVCHM Pascal text source file, including comments, is 39.3K bytes, and that of CMPLO is 1 1 . 6 ~ b y t e s . The executable program files for HAVCHM, CMPLO, and the alternative CMPLSTD are 40.OK, 44.1K, and 22.9K bytes, respectively. The Include statements bring many routines into CMPLO from the Turbo Graphix package; storage space was sacrificed for convenience of the programmer. The Turho-87 Pascal system for use with the 8087 processor compiles programs to use and store real numbers as 8bytel64-bit words, and executable programs and files created in this format cannot he interchanged with executable programs and files using the more common 4-bytel32-bit word length. Pascal text source files, two sets of test parameter and output files, compiled HAVCHM program, three compiled versions of CMPLO plot program for different hardware configurations, and compiled CMPLSTD plot program are available on two MS-DOS formatted disks along with documentation for use. One disk contains programs and files compiled for the 8-byte real number format; the other disk contains oroerams and files com~iledfor the 4-hvte real number firmit. Program listings may be extracted from the MS-DOS Droeram text files. The HAVCHM package is ject available f i ~ ~ ~ r o SERAPHIM. Literature Clted I. Edelmn. D. J. Chom. Educ. 1915.52.642. 2. Field. R. J.; Nayes, R. M. J. Chem. Edur. 1914.60,1877. 3. Unpublished work of author; literature source not cited here toproted the guilty. 4. Stabler.R.N.:Chesiek.J.P.1nt.J.Chom.Kindicr 1978.10.461. 5. Field, d. J.J chChem. due. 1981.58,408. 6. Gear, C. W. Numerical Initial Value Problem in Omlinory Di//erenfial Equoliom; Prentics-Hall: Englewood Clim,NJ, 197L; Chapter 11. 7. Gear.C. W. Commun.ACM1971.14.176.1S5. 8. Chpaiek. J. P. J. Chem. Edur 1919,56,585.
