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.