whitbeck1@llnl.gov DRAFT 1:58 PM 1 October 16, 1995 Numerical Modeling of Chemical Kinetic Reaction Mechanisms with Least-Squares Fitting to Experimental Data PURPOSE A program for modeling complex systems of chemical reaction is presented for use on Macintosh computers. Look elsewhere at this site for DOS and OS/2 compiled executables. The program reads in a reaction mechanism file and builds the descriptive set of differential equations which are then solved subject to integration parameters from a separate parameter file. If experimental data is available it may optionally be used to adjust one or more rate constants in a least-squares sense. HISTORY The origins of this program date back to ca 1973 to a FORTRAN program I wrote at OHIO STATE as a chemistry graduate student. I used that program for modeling flash photolysis experiments and atmospheric photochemistry. That program used Gear's ODE solver. I had no connection to the internet or arpanet yet the program found a home at several prestigious sites under various names. A brief stint in industry led to increased sophistication, non-linear least-squares using Marquardt's algorithm, sensitivity analysis to help identify key reactions and species, coupled algebraic and differential equations handled... Most of these features were only lightly used. My employer considered this version proprietary and would not distribute it outside the company. On to pseudo-academia. I left industry for an R&D position affiliated with a state university. Here the program gained the ability to handle heterogeneous chemical kinetics- namely the growth of a hydrosol droplet in air coupled with reactive scavenging of particles and gases. Still rather awkward FORTRAN. The INTERNET connection. While in this job I sent out numerous printouts of my original (grad student) FORTRAN program to eager volunteers who wanted to use this program.... I never heard back. Still the requests kept coming so I cobbled up a quick & dirty C version made available to the internet by the Ohio Supercomputing Center. Hundreds of chemists and engineers around the globe have used the program. This current release is a complete rewrite using few if any machine dependencies, accordingly look for DOS, OS/2 mutations at this site as well. I am now in pseudo-government service (i.e. a national lab) my programming and support for this project is thus limited. I have tested the program on my home computer (a Performa 6115CD, 601 risc processor) and that's about it. This application was not programmed for speed but rather convenience using an object oriented approach as much as possible. Yet with modern desktop computers it runs pretty fast- much faster than the old punched deck program ran in '73! I hope you find it useful. USAGE The program uses a master file rather than rely on UNIX conventions1 or machine-specific graphical interfaces. Input thus consists of 4 or 5 files: a master file containing a list of other files a parameter file for integration parameters a file of initial concentrations/masses a file containing the reactions and an optional name of a file of experimental data Files: master this file has one item per line in this order: name of the reaction mechanism file name of the parameter file name of the initial values file name used to save model results code word FIT or MODEL to define problem type and an optional name of a file of experimental data a sample master file: Mechanism_test_1 Parameters Initial_values_1 Output FIT observed mechanism this file has one reaction per line each reaction consists of a rate constant, a list of reactants, a reaction "ARROW" and a list of products all as whitespace (blank or tab) delimited tokens. The reaction ARROW is defined as the string "-->" (no quotes). The reactants and products may optionally be separated with a whitespace delimited " + ". e.g.: 1.03e+02 A + B + C --> D + E Chemical names should be less than 99 characters, lines less than 1999 characters long. No limit is set on the number of reactants or products in a reaction. Remember, gigo- garbage in -> garbage out: this program is no substitute for knowing chemistry! A whitespace delimited token " *" designates that reaction as one whose rate constant may optionally be adjusted to model experimental data. A limit of 200 such reactions is imposed (gigo!). a sample mechanism file: This mechanism is a test with an analytic solution 0.5 * A --> B 0.10 B --> C a silly test Lines lacking an ARROW are treated as comments. initial values file has a list of species using the names used, the initial value and a tolerance used by the solver e.g.: A 10.0 1.0e-10 The defaults used for initial value is 0.0 and 1.0e-8 for tolerance. The integration will attempt to bound the error in the dependent variables using weights of 1/(R*y+A) where y is the dependent variable, R is a relative tolerance (nominally 1e-8) and A is the absolute tolerance specified in the input; a different value (greater than or equal to zero) may be given to A for each species however A should not be zero if the integrated variable is ever zero! parameters file contains the initial time, time step, final time and the maximum number of output time steps (for safety!- gigo). If the optional least-squares fitting is being done a second line is read containing the fit parameters: a convergence term (0-1, default 0.5), an accuracy term ( default 1e-8) and a maximum number of iterations (default 100); negative entries invoke the default. This second line is ignored if the problem type is MODEL rather than FIT. a sample parameter file: 0.0 0.1 5.0 50 -99.0 -99.0 100 observed data file contains the experimental data in column format. Each column heading must be a valid species name (used in the mechanism) except for one of which must be "TIME" (no quotes). The first row is the weight to be applied to that variable in fitting (the weight for TIME is not used but a value must be given). Subsequent rows define the values observed at the specified times. a sample observations file: A B C TIME 2.0 1.0 1.0 1.0 9.048376e+00 9.468032e-01 4.821347e-03 1.000000e-01 8.187308e+00 1.794087e+00 1.860453e-02 2.000000e-01 7.408178e+00 2.551419e+00 4.040276e-02 3.000000e-01 6.703192e+00 3.227447e+00 6.936095e-02 4.000000e-01 6.065304e+00 3.829989e+00 1.047068e-01 5.000000e-01 ... 1.108118e+00 7.685633e+00 1.206249e+00 2.200000e+00 1.002670e+00 7.714074e+00 1.283257e+00 2.300000e+00 The first observation must not be the initial value (TIME=0 typically). Values must be sorted on increasing TIME. The order of appearance of variable names and values correspond but other than that any order is permitted e.g.: B TIME A C 1.0 1.0 2.0 1.0 9.468032e-01 1.000000e-01 9.048376e+00 4.821347e-03 is just as valid. OUTPUT to the screen looks like: React3 ... The mechanism is Mechanism_test_1 The parameters are from Parameters The output is to Output The initial values are from Initial_values_1 The observations are from observed ::This mechanism is a test with an analytic solution :: a silly test 2 reactions involving 3 species 1 adjustable reactions Iterating over a mechanism 5.000000e-01 adjustable A --> B 1.000000e-01 B --> C Least-squares fitting of the mechanism to 23 observations Least-squares fit used 26 iterations The new model is: Iterating over a mechanism 1.000000e+00 adjustable A --> B 1.000000e-01 B --> C finished the output file looks like: A B C t 1.000000e+01 0.000000e+00 0.000000e+00 0.000000e+00 9.048376e+00 9.468031e-01 4.821346e-03 1.000000e-01 8.187308e+00 1.794087e+00 1.860453e-02 2.000000e-01 7.408178e+00 2.551419e+00 ... As a test the model was run using synthetic "experimental" data: Another fitting example: The mechanism is tasmech The parameters are from tas2params The output is to tasmania The initial values are from tas2initials The observations are from tasdata ::This is the tas mechanism file 3 reactions involving 3 species 2 adjustable reactions 1.000000e-02 adjustable a --> b 2.000000e-02 b --> a 2.000000e-03 adjustable b --> c Least-squares fitting of the mechanism to 4 observations Least-squares fit used 26 iterations The new model is: 2.999735e-02 adjustable a --> b 2.000000e-02 b --> a 3.000302e-03 adjustable b --> c finished POC: Dr. Michael Whitbeck Chemistry & Material Science 510-424-4472 whitbeck1@llnl.gov whitbeck1@popcorn.llnl.gov CREDITS This program uses Hindmarsh and Cohen's CVODE algorithm for integration of ordinary differential equations and Johnson's implementation of the Hooke-Jeeves simplex algorithm for fitting experimental data. APPENDIX The "command line" invoked at startup can have additional entries: "React 3" master_file mxstep mxord reltol where "React 3" is the program name and prompt in the dialog box master_file is the name of the master control file mxstep is the maximum number of steps used by the ODE solver mxord is the maximum order used by the solver (5 default, reduce if necessary) reltol is the relative tolerance (scalar) applied by the solver Adjustments to the absolute tolerances (in the initial values file), relative tolerance, mxsteps, and mxord allow solution of stiff problems. 1The Macintosh interface takes the name of the master control file from the "command line" that opens as a dialog box. This dialog also provides the ability to reroute "stdout" to a file or printer. There are a few secret entries as well- see the appendix.