/* The permission is granted to use the software as is and redistribute it in its original and complete form. If you find this program useful and publish results obtained by it, please cite the program as: Michael Whitbeck, Desert Research Institute, 1991, "REACT - Program to solve kinetic equations for chemical systems". Version 1.00, this code is available from the authors or from Dr. Whitbeck. */ This is the REACT package release 2. A while back someone asked in sci.chem about chemical kinetics simulation software. Well I hacked out a program that just does the simulation (and parameter optimization). It is written in C for the Sun 3|4 computers (should port to other platforms- I ported an intermediate version to GEMDOS (1040ST) with no great difficulty). The code uses lsoda as the core integrater = livermore solver for ordinary differential equations, with automatic method switching for stiff and nonstiff problems by Linda R. Petzold and Alan C. Hindmarsh, computing and mathematics research division, l-316 lawrence livermore national laboratory (lsoda.f is available at netlib@research.att.com, the C port was provided by Hon Wah Tam Wolfram Research, Inc. tam@wri.com). INPUT: two files root.m where root can be your choice of names This file contains a reaction mechanism (free format) like: 2.0 a -> b 0.15 b -> c 0.80 b + c -> d + d the first field is a rate constant in 'F' or 'E' format (e.g. 1.00e-09) then a list of reactants and products separated by '->' (required). The specie names can be anything, eg CH3O2, NO2, PABA etc. There are some restrictions on name lengths (10 chars) and some tokens are 'special'. root.p is a file containing the integration parameters again in a free format. It looks like: 0.0 6 .1 1.0E-9 a 1.0 1.0e-10 1.0e-6 b 0.0 1.0e-10 1.0e-8 You simply provide (first line) starting time, stopping time, output interval, and a solution tolerance value (a default value for solution accuracy- to be explained in the package). This line is followed by one or more lines to initialize specie values (a,b,c,d in the example). No entry is required if the initial value is zero and the default tolerance is acceptable. Otherwise give the species name, starting value, absolute tolerance, relative tolerance (explained in the package). The program REACT sucks in these two files and spits out a file named root.r (again root is from your input file name; i.e. REACT A produces A.r from A.[m,p]). This file contains the answers! To facilitate printing and plotting the package contains 3 post processors: SELECT reads root.e, a file containing the species to be printed. if root.e contains the lines a b d c b then a table is printed with the headings TIME a b d then a table is printed with the headings TIME c b a sample from the above example might be: 3 reacions involving 4 species a b c d time a b c 0.0000e+00 1.0000e+00 0.0000e+00 0.0000e+00 1.0000e-01 8.1873e-01 1.7987e-01 1.3925e-03 2.0000e-01 6.7032e-01 3.2439e-01 5.1475e-03 3.0000e-01 5.4881e-01 4.3989e-01 1.0669e-02 4.0000e-01 4.4933e-01 5.3151e-01 1.7433e-02 .... 6.0000e+00 6.1442e-06 2.4471e-01 1.7167e-01 EXTRACT same as select but no table headings, useful with unix tools like dm, awk, etc. EXTRACTER dumps all the data (prints to stdout) in ASCII starting with time. A fifth file root.s contains the list of species that goes with root.r so that root.m does not have to be reprocessed to make these tables. =============================== to make it first compile the lsoda integrater then edit this makefile to reflect where you are keeping the integrater code then make react. If that succeeds make select extract extracter. NOTE: as time goes buy ports may be made to other platforms; check the makefile for defines for your system! Send correspondance, bug reports to ___________________________________________________________ |Mike Whitbeck | whitbeck@unssun.unr.edu | |Desert Research Inst. | whitbeck@wheeler.wrc.unr.edu | |POB 60220 | whitbeck@sanjuan.UUCP | |RENO, NV 89506 | 702-673-7348 | |__________________________|______________________________| ======================================================================== ============================UPDATE====================================== ======================================================================== ========================================================================= ============additional info for least-squares fitting==================== ========================================================================= REACT can optimize the rate parameters used in the model if it is given a file root.f containing the measured data for fitting. Add the keyword 'fit' to invoke this option. a sample root.f file for react version 2 (pre-release) (do not enter the <- marker or anything after it) 100 <- max number of iterations allowed, CAREFUL! 3 3 <- no. parameters to fit no. observed species to use 1 2 3 <- list of subscripts for k_i's a b c <- list of specie names that are observed 0.01 0.01 0.001 <- step sizes to be used in hunting for the answer 1e-2 1e-2 1e-2 1e-1 <- max. error (fractional) allowed/sought 0.0000e+00 1.0490e+01 4.1000e+00 0.0000e+00 < d 2.0000e+01 7.5018e+00 6.7479e+00 3.4029e-01 < a 4.0000e+01 6.3005e+00 7.5160e+00 7.7354e-01 < t 6.0000e+01 5.7445e+00 7.6161e+00 1.2294e+00 < a 8.0000e+01 5.4238e+00 7.4833e+00 1.6829e+00 < the input data for the observables are in the order time and then concentrations for the measured species in the order given in the list above. ==========================sample output======================= equinox:141 src> time react fit tas2 1) 9.00e-02 a -> b 2) 9.00e-02 b -> a 3) 9.00e-03 b -> c SIMPLEX Optimization Program exited after 113 iterations. The mean is: 3.000023e-02 2.000013e-02 2.999969e-03 5.258171e-09 The estimated fractional error is: 1.747761e-05 2.325398e-05 1.065852e-05 7.440557e-02 The standard deviation is 6.985651e-05 The estimated error of the function is 2.016584e-05 4.030u 0.220s 0:05.61 75.7% 0+235k 3+7io 0pf+0w ==========================sample output======================= ----explanation----- IN THIS EXAMPLE I specified that k's k_1, k_2, and k_3 were to be optimized using measured data for [a], [b], and [c]. The initial k's used to start the model were: .09, .09, and .009 The model found the best values to be : 3.000023e-02 2.000013e-02 2.999969e-03 (from the line 'the mean is:') With a new model first get it running with some rough guesses for the parameter(s) to be fitted. Next try a loose fit (few iterations low maxerr criteria) just a few data points closely spaced (e.g. at t= 0, 10, 20, 30, 40). This will give you a better starting point. Improve the starting guesses for the fitted k's using this output then fit with as much data as you have. If convergence is unsatisfactory re-think the model, re-check the data entry, etc. Ther may/propbably will be changes and possibly bugs will be uncovered- report bugs to whitbeck@unssun.unr.edu