KINAL : A Program Package for Kinetic Analysis of Reaction Mechanisms User's manual version 2.1 Tamas Turanyi Department of Physical Chemistry E”tv”s University (ELTE) Address: H-1518 Budapest-112, P.O.Box 32, Hungary Phone : (36-1) 209-0555 ext. 1109 Fax : (36-1) 209-0602 E-mail : turanyi@garfield.chem.elte.hu WWW : www-PhCh.chem.elte.hu Abstract - A program package is provided for analysis of kinetic mechanisms on personal computers. KINAL consists of five programs called DIFF , SENS , PROC , ROPA and YRED. They require similar input data and use common subroutines. DIFF solves stiff differential equations and SENS computes the local concentration sensitivity matrix. PROC generates the rate sensitivity matrix or the quasi-stationary sensitivity matrix from concentration data or uses matrix computed by SENS and extracts the kinetic information inherent in sensitivity matrices by principal component analysis. ROPA calculates the contribution of reaction steps to the production rate of species or calculates the lifetime of species. Finally YRED provides hints for the elimination of species from the reaction mechanism. DIFFDAT, a utility program to KINAL, is an interactive program for the easy creation of data files. Introduction Simulation of complex reaction mechanisms play an important role in chemical kinetics. Investigation of mechanisms includes the solution of kinetic differential equations , the study of the effects of parameters change on the solution ( sensitivity analysis ), the exploration of important reaction pathways and the identification of redundant species and unimportant reactions. The aim of the development of package KINAL was to provide a solution of all these problems by using the same code and the same input data file. Background A homogeneous isothermal complex chemical process can be described by a system of kinetic differential equations : d c / d t = f (c,k) ; c(0) = c0 (1) Unfortunately analytical solutions are available only in the simplest cases and the numerical solution is often difficult because system (1) is usually stiff. The first successful stiff ODE solver was based on a changing order backward differentiation formula (Gear,1971). GEAR and its successors (like EPISODE ; Hindmarsh & Byrne, 1977) are still the most widely used codes. Recently the implicit Runge-Kutta methods were shown to be very efficient in solving stiff problems (Hairer&Wanner,1981). A fourth order semiimplicit Runge-Kutta method has been published for the solution of stiff ODEs occuring in chemistry and biology (Gottwald&Wanner,1981). Gottwald and Wanner (1982) claim that their code is better than GEAR and EPISODE, although there is also a contrary opinion in the literature (Seifert,1987). An important question in each modelling study is the effect of parameter change on the solution. This emerges for instance when the uncertainty of solution caused by the uncertainty of parameters is investigated. Sensitivity information is also very useful when carrying out an experimental design or a parameter estimation. An element of the local concentration sensitivity matrix S = { d ci / d kj } is the linear approximation of concentration change of species i at time t2 caused by the change of parameter j at time t1. This matrix function is the solution of the following differential equation: d Sj / d t = J Sj + Fj j=1,2,... m (2) where J is the Jacobian of equation (1) and Sj and Fj are the j-th coloumns of matrices S and F , respectively, i.e. Sj = dc / dkj and Fj= d (dc/dt) / d kj. A number of methods have been published for the effective numerical solution of (2). The best known of them is the Green's function method (Hwang et al,1978). Another family of methods utilize the fact, that the Jacobians of equations (1) and (2) are the same, and therefore the decomposition of matrix J has to be carried out only once in each step. Such methods, like (Dunker,1984), (Valko&Vajda,1984) are shown to be more accurate and more computer time saving then the Green's function method unless the number of parameters is very large. Often one is interested in the effect of a parameter change on the concentrations of a group of species. Such a measure can be introduced by utilizing objective functions. The effect of a single parameter on a group of concentrations is demonstrated by the overall sensitivities (Vajda et el,1985). A better description of parameter-concentration interdependence consists of the identification of the groups of parameters which have joint influence on a group of concentrations. This type of information is given by principal component analysis of the lognormalized concentration sensitivity matrix (Vajda et al,1985). Eigenvectors of matrix STS identifies parameter groups while the eigenvalues inform us about the effectiveness of these parameter groups for the change of the concentrations of species, occuring in the objective function of the method. A parameter is considered important, if it belongs to a large element of an eigenvector corresponding to a large eigenvalue. Redundant reactions of a mechanism can be identified by concentration sensitivity analysis (Vajda et al,1985), (Vajda&Turanyi,1986). This approach requires considerable computer time and the information gained belongs to the (t1,t2) time intervall and therefore it is impossible to relate reaction significances to a concentration set. However, this approach proved to be successful in the study of steady premixed laminar flames (Vajda et al, 1990). An alternative way for the identification of redundant reactions is the use of rate sensitivity matrix F (Turanyi et al,1989). This approach may be considered a more sophisticated handling of rate information than the often applied competitive kinetic considerations. Matrix F can be processed in a similar way to that shown for the concentration sensitivity matrix using overall sensitivities or principal component analysis. Another method for the identification of the kinetic significances of reactions at a given concentration set is to apply the principle of quasi-stationarity to the system of differential equations (2), i.e. to solve the corresponding algebraic equations instead of the original system of ODEs. This treatment is called quasi-stationary sensitivity analysis (Turanyi et al , 1988). This approach is exact for stationary systems but it is well applicable for non-stationary systems too. In some cases the objective of modelling is to describe the concentration of only a sub-group of species, which are the so called important species. A minimal reduced mechanism for the calculation of the concentration-time curves of important species can be obtained by adopting important and necessary species into and by excluding redundant species from the objective function of principal component analysis of the above described sensitivity matrices. Redundant species can be identified (Turanyi, 1990) by one of the following methods: A species can be considered redundant if the elimination of its consuming reactions does not cause significant deviation from the solution of the full model when inspecting important species. Some of the remaining species are also redundant provided that the simultaneous elimination of both the fastest forming and consuming reactions have no significant effect. The second method is based on the investigation of the Jacobian. An overall sensitivity type measure shows the effect of concentration change of each of the species on the rate of concentration change of the important species. Redundant species are obtained by an iteration procedure, described in details at the discussion of subroutine YRED. The first method directly shows the deviation between the solutions of the full and reduced mechanisms caused by the elimination of a species. The second method estimates the significances of species at a given concentration set. *************************************************************************** * * * The tools of sensitivity analysis and its applications in * * chemical kinetics including the investigation and reduction * * of reaction mechanisms were recently reviewed: * * * *-------------------------------------------------------------------------* * * * Turanyi T.: Sensitivity analysis. Tools and applications * * J.Math.Chem.,5,203-248(1990) * * * *-------------------------------------------------------------------------* * * * This review article contains 175 references. * * * *************************************************************************** The KINAL program package There are five programs in the program package KINAL. DIFF calculates concentration-time curves by solving the system of kinetic differential equations (1) or compares full and reduced mechanisms. SENS computes the local concentration sensitivity matrix S. PROC reads matrix S or reads concentrations evaluated by DIFF and generates the rate sensitivity matrix F or the quasi-stationary sensitivity matrix QSA and processes the matrices by principal component analysis. ROPA calculates the contribution of reaction steps to the rate-of-production of species from concentration data and presents this information in detail or in a concise form. ROPA can also calculate the lifetime of species. YRED also reads concentrations and identifies redundant species by the method based on the investigation of the Jacobian. The program package is written in Fortran77 language and includes five main programs and 21 subroutines. Most of the subroutines are used by at least two main programs. The program package consists of 3050 Fortran lines including 970 comment lines. Sample batch files ( DLINK.BAT, SLINK.BAT, PLINK.BAT , RLINK.BAT, YLINK.BAT and ALLPRINT.BAT) , sample data files (RED.DAT, DTEST.DAT, STEST.DAT, PTEST.DAT, RTEST.DAT and YTEST.DAT) and a sample output file (OUT.PUT) help the use of the program package. The present form of the program package can handle a mechanism of not more than 90 reactions and 50 species. In one run maximum 20 concentration sensitivity coefficients can be calculated. The program computes the solution of the system of kinetic differential equations and constructs the sensitivity matrices at maximum 60 points of time. If the investigation of a larger mechanism is required , the program should be redimensioned. In order to facilitate this, the arrays to be redimensioned are listed at the beginning of each subroutine. Detailed comments help the modification of the code. However, use of the program package doesn't require programming knowledge. Capacity of personal computers are usually enough to solve reaction kinetic problems and the application of such machines may speed up the work. KINAL has been developed on an IBM PC compatible personal computer, but the package doesn't utilize any special services of PCs (like for instance graphics) and therefore it is transferable to other computers. In case of transfer you have to change subroutine RTIME, which reads the clock of the computer and is computer and compiler dependent. 1. Input data Each KINAL program asks for the name of the input and output files. If you have an MS DOS like operation system, you can also give PRN or CON as output file names. This way you can direct the output to the printer or to the screen, respectively. The five programs use very similar input data files and only slight modifications of them are necessary when changing from one program to another. Input data includes the name of the problem , the maximum eight character identification of species , initial concentrations , parameter values , etc. Stoichiometry of a reaction system can be represented in a data file in several ways. Typing the original chemical equations is the closest to the cast of mind of a chemist user, but writing them requires considerable effort , as the formulas of species has to be typed again whenever it appears in a reaction. A machine-oriented, but more type-effective way is to relate the species to numbers. In the KINAL data files stoichiometry is provided in two blocks. The first block contains information about the left hand side of chemical equations and the second one represents the right hand side. Each row of the two blocks consists of the number of reactions in which a given species reacts or in which it is produced, the identification number of the given species and its stoichiometric number. In the first block the stoichiometric numbers are all integers , in the second one real numbers are also permitted in order to make possible the investigation of lumped mechanisms. If no stoichiometric number is given in a row of the first block , than it is assumed to be one. The blocks are closed by all-zero lines. Stoichiometry is stored in packed form , so even in case of several hundred reactions only few memory spaces are consumed. In contrast to usual kinetic ODE solvers, when applying KINAL the number of species is not limited on either side of the reaction equation. A reaction number-species number pair has to be represented only once in a data file. Hence the A + A --> type reactions have to be coded in 2 A --> form. Input data are checked by the program and if they are incredible the user is warned or the run is interrupted. Data evaluated by DIFF and SENS are written into the end of the input data file. Therefore the first part of the data file to be read by PROC (or by the other four programs) is written by the user , the second part by DIFF or SENS. Preserving the original data in the first part of data file facilitates the identification of sensitivity information. KINAL has been extended by an interactive QuickBASIC language program named DIFFDAT for the easy creation of data files. This program (in its original form) was made by Prof. G. Huybrechts and S. Merdjan ( Faculteit der toegepaste wetenschappen, Laboratorium voor fysische scheikunde, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussel, Belgium). Using this program, it is not necessary to keep in mind the input format and the sequence of data. The program asks for the data and creates a FORTRAN input file. If the reaction times are allocated equidistantly, you have to give only the initial time and the increment. If you wish to study a reaction mechanism at several temperatures, you have to give the Arrhenius parameters only at the first time and the program calculates later the rate coefficients at different temperatures as well. This program includes an easy way of typing the chemical equations in. It provides the table of species and you can select a species by arrows and the key. When typing the left hand sides, a reselection increases the stoichiometric number. On the right hand side the stoichiometric numbers are asked explicitly. A blank enter response is considered as one stoichiometric number. You can change from the left hand side to the right hand side or from the right hand side to the left hand side of the next chemical equation by pressing . It is possible to delete incorrectly typed data (only if the key has not been pressed) by using key . 2. Kinetic subroutines Chemical kinetic problems can be investigated by KINAL without changes in the computer code. Beside the main programs, the following subroutines employ the concept of mass action kinetics : FUN , JFY , JFP , JF2YY , JF2YP and EQS. EQS prints the list of reaction equations which is useful in checking the stoichiometry. Note that the printer control characters of this FORTRAN subroutine act properly only if the output of the program goes directly to the printer. FUN evaluates the system of ODEs from stoichiometry and parameters contained in ST common. JFY and JFP computes the Jacobians with respect to variables and parameters , respectively ; JF2YY and JF2YP computes the appropriate second derivative tensors. The study of some physical , chemical and biological problems requires the solution of differential equations and their analysis. KINAL is suitable for the investigation of any problems described by ODEs, if the right hand sides of the differential equations and their derivatives have been written into the subroutines listed above. KINAL can be combined with other kinetic simulation programs, too. For example concentrations of a non-isotherm and/or inhomogeneous reaction system can be calculated by a suitable computer code and than redundant species and redundant reactions of the utilized mechanism can be identified by slightly modificated versions of YRED and PROC. 3. DIFF : Solution of ODEs The solution of kinetic differential equations in KINAL is based on the ROW4A subroutine (Gottwald & Wanner, 1981,1982). The program has three operation modes. In mode zero it only solves the differential equations while in mode one DIFF also writes the computed concentrations into the input data file. In mode two the program reads the previously computed concentrations, performs the calculations of a reduced reaction system and compares the solutions of the full and reduced mechanisms. Reactions to be weeded out from the reaction mechanism should be indicated by a zero in RED.DAT data file. This option is useful for checking reduced models and for the identification of redundant species by the first method as described in the Background, section above. 4. SENS : Computation of local concentration sensitivities In a kinetic problem rate coefficients and initial concentrations are considered as parameters. SENS computes the local concentration sensitivities according to both type of parameters. It calculates sensitivities either for all of the parameters or only for a selected group of them. The latter enables the investigation of a large model in parts when using a computer of very limited memory. SENS can also continue the computation of sensitivities by reading the previously calculated sensitivity matrix. Computation of local concentration sensitivities in SENS is based on the Decomposed Direct Method (Valko&Vajda,1984). Subroutine ROW4SP of KINAL may be considered a PC adaptation of the original subroutine of Valko and Vajda. The original code was optimized for mainframe computer calculations. It preserves all the arrays, that can be used in a later stage of computations and hereby consumes a great amount of memory but saves computer time. In the original ROW4S the economy of memory is acheaved by applying a large working vector. However handling an array greater then 64K is ineffective on 16 bit PCs. Generally the limited memory space of the PCs and not the computation time is the critical factor, the ROW4SP subroutine of KINAL is organized to decrease the memory requirements on the account of computer time. The greatest part of computer time is consumed by the evaluation of matrix d2 f / d y2 . This matrix is constant provided that none of the reaction orders is greater than two and the matrix needs to be computed only once at the beginning of the run. Modify the code in that case ! 5. PROC : Processing of sensitivity matrices PROC carries out a principal component analysis of sensitivity matrices and provides a list of overall sensitivities and reaction rates. Sensitivity matrices belonging to different reaction times can be investigated by PROC either separately or simultaneously. The former is suggested in the analysis of matrices F and QSA, while the latter is used at the processing of local concentration sensitivity matrices as described by Vajda et al (1985). Principal component analysis incorporates eigenvalue - eigenvector decomposition of the cross-product of normed sensitivity matrices. In KINAL eigenvalues and eigenvectors are computed by subroutine SDIAG2, based on a QR algorithm, which is sufficiently stable and provides suitable accurate eigenvectors. PROC lists eigenvalues (in decreasing order), the corresponding eigenvectors and the reaction numbers belonging to eigenvector elements. PROC also proposes a reduced mechanism by identifying reactions which belong to large elements in eigenvectors corresponding to large eigenvalues. An eigenvalue or an eigenvector is considered "large", if it is greater than a threshold value. Recommanded threshold values are are 10(-4) and 0.1 or 0.2 for eigenvalues and eigenvectors, respectively. The best thresholds may be different for different problems and PROC provides the opportunity to try several limits in an interactive way. 6. ROPA : Rate-of-production analysis and calculation of the lifetime of species The classic method for the study of the importance of reactions is the rate-of-production analysis (see e.g. (Gelinas,1979) and (Kee et al.,1985)). Although the combination of species reduction and rate sensitivity analysis (Turanyi, 1990) seems to be a more effective way for this purpose, the rate-of-production analysis is still an important method for the exploration of important reaction pathways. The rate-of-production analysis requires the calculation of the Pij matrix elements (Gardiner,1977), which show the contribution of reaction j to the rate of production of species i. The rate-of-production information is rather difficult to look over even in an ordered form, as it is given by ROPA using mode 0. Therefore ROPA can also provide (mode=1) the rate-of-production information in a concise form. In this case ROPA do not list the actual Pij values, but provides the ordered list of the numbers of the most influental consuming and producing reactions for the rate of production of a species. The explanation of abbreviations is given in the beginning of program ROPA. A generalized lifetime of species i is the negative reciprocal of the i-th diagonal element of the Jacobian /J={d ci/ d cj}/. Program ROPA (mode 3) gives the ordered list of lifetimes of species. It has been shown by Frank-Kamenetskii (1940), that species having the shortest lifetime are candidates for the application of the quasi-steady-state approximation. Recently it was shown (Turanyi et al, 1992), that the product of the lifetime and of the production rate was a better criterion for the selection of the QSSA species. Note, that nowadays the application of QSSA became widespread for the reduction of combustion models (Peters&Rogg, 1993) 7. YRED : Identification of redundant species YRED estimates the effect of change of concentrations of each species upon the concentrations of a selected group of species. A species is considered to be redundant if both its direct and indirect effects on important species are small. (Species A affects directly the concentration of species B, if A is consumed in a reaction in which B is consumed or produced.) Direct effects are estimated by an overall sensitivity type measure, which is the sum of squares of normalized Jacobian elements. Indirect effects are taken into account by an iteration procedure. First the overall sensitivities are calculated from the normalized Jacobian elements extending summation over the important species. Then the group of important species is supplemented by species having the highest overall sensitivities among the rest of the species and therefore being considered necessary. The overall sensitivities are calculated again by taking into account the new species as well in the summation. In this way further necessary species are identified. It is recommended, that only a few new species are adopted into the list of necessary species at any one time; the number being controlled by the choice of a threshold value. To facilitate the choice of this threshold, an asterisk behind formulas in the table of overall sensitivities indicates that the marked species has been considered in the summation. The described procedure is repeated until the iteration converges and a gap appears in the series of overall sensitivity values between the necessary and redundant species. (Some important species may have little or negligible effect on other species and they are bad ranked in the overall sensitivity table.) Species taken into account in the summation at this stage are the important and the necessary species while the rest are the redundant species. 8. Test of the program package As a test of the work of the program package sample data files are enclosed for all the five programs. These files contain the highly stiff Chapman mechanism (three species and four reactions). The results of the programs are presented in the file OUT.PUT . Calculated concentrations and sensitivity informations are available for that mechanism in the literature (see e.g. Valko&Vajda,1984). Results of the investigation of large mechanisms by principal component analysis were published in the following articles: Topic: Reference: H2-Br2 reaction Vajda et al, 1985 Formaldehyde oxidation Vajda et al, 1985 EFN model of the BZ reaction Vajda&Turanyi, 1986 Formaldehyde oxidation Turanyi et al., 1989 High-temp. propane pyrolysis Turanyi et al., 1989 Atmospheric chemistry Turanyi&Berces, 1990 Low-temp. propane pyrolysis Turanyi, 1990 Briggs-Rauscher reaction Turanyi, 1991 GTF model of the BZ reaction Turanyi et al., 1993 Methanol thermolysis Borger et al., 1992 Methylene blue reaction Zhang&Field , 1991 Decomp. of silicon compounds Ahmed et al. , 1991 BZ reaction Gyorgyi&Field, 1991 Landolt reaction Ibison , 1992 Hydrogen combustion Tomlin et al., 1992a Methane combustion Tomlin et al., 1992b The first three calculations were made before the creation of KINAL. They were recalculated by KINAL giving the same results. The listed articles also contain practical hints for mechanism investigation and reduction by the principal component analysis of sensitivity matrices. References: Ahmed M., Davidson I.M.T., Morgan G.H., Simpson T.(1991), Organometallics,10,3772 Borger I., Merkel A., Lachmann J., Spangenberg H.-J., Turanyi T.(1992) Acta Chim. Hung., 129, 855-864 Dunker,A.M.(1984),J.Chem.Phys.,81,2385 Frank-Kamenetskii,D.A.(1940), Zh.Fiz.Him.,14,695 English translation is found in: Turanyi T., Toth J.(1992), Acta Chim. Hung., 129, 903-914 Gardiner,Jr.,W.C.(1977),J.Phys.Chem.,81,2367 Gear,C.W.(1971),Comm.ACM,14,185 Gelinas,R.J.(1979), Science Applications, Inc., Preprint No SAI/PL/C279. Gottwald,B.A.&Wanner,G.(1981),Computing 26,355 Gottwald,B.A.&Wanner,G.(1982),Simulation 37,169 Gyorgyi L., Field R.J.(1991), J.Phys.Chem.,95,6594 Hairer,E.&Wanner,G.(1981),SIAM J. Num.,18,1098 Hindmarsh,A.C.&Byrne,G.D.(1977),EPISODE. An effective package for the integration of systems of ordinary differential equations. Lawrence Livermore Laboratory Report UCID-30112,Rev.1. Hwang,J-T.,Dougherty,E.P.,Rabitz,S.&Rabitz,H.(1978), J.Chem.Phys.,69,5180 Ibison P.(1992), J.Phys.Chem.,96,6321 Kee,R.J.,Grear,J.F.,Smooke,M.D.&Miller,J.A.(1985), Sandia National Labs., SAND 85-8240 Peters N., Rogg B. /eds./, (1993) Reduced kinetic mechanisms for applications in combustion systems Springer Verlag, Heidelberg Seifert P.(1987),Computing,38,163 Tomlin A.S., Pilling M.J., Turanyi T.,Merkin J.H., Brindley J. (1992a), Combust. Flame, 91,107 Tomlin A.S., Pilling M.J., Turanyi T., Merkin J.H., Brindley J. (1992b) Strategies for the Reduction of Complex Chemical Mechanisms Twenty-Fourth Symposium (International) on Combustion Sydney, Australia, 5-10 July 1992 Turanyi T.,Berces T.&Toth J.(1988),J.Math.Chem.,2,401 Turanyi T.,Berces T.&Vajda S.(1989),Int.J.Chem.Kinet.,21,83 Turanyi T., Berces T. (1990) React. Kinet. Catal. Lett. , 41, 103 Turanyi T. (1990),New J.Chem,14,795 Turanyi T. (1991), React.Kinet.Catal.Lett., 45, 235 Turanyi T., Tomlin A.S., Pilling M.J.(1993) J.Phys.Chem., 97, 163-172 Turanyi T., Gyorgyi L., Field R.J. (1993) J.Phys.Chem, 97, 1931-1941 Vajda S., Rabitz H.&Yetter R.A.(1990),Combust.Flame,82,270 Vajda S.&Turanyi T.(1986),J.Phys.Chem.,90,1664 Vajda S.&Valko P.&Turanyi T.(1985),Int.J.Chem.Kinet.,17,55 Valko P.&Vajda S.(1984),Comp.&Chem.,8,255 Zhang Y-X., Field R.J.(1991), J.Phys.Chem.,95,723 *************************************************************************** * * * Please refer to the below article, if you have used KINAL : * * * * T. Turanyi * * KINAL: A program package for kinetic analysis of reaction mechanisms * * Computers & Chemistry, 14(3),253-254(1990) * * * ***************************************************************************