2315
Computer Simulation to Complicated Chemical Systems
Application of Computer Simulation to Complicated Chemical Systems Richard M. Noyes Deparfment of Chemistry, University of Oregon, Eugene, Oregon 97403 (Received April 26, 1977) Publication costs assisted by the National Science Foundation
If all rate processes and parameters are known, existing computational techniques are adequak to model chemical change in very complicated uniform systems and also in isothermal systems of high symmetry such that position can be defined by a single parameter. However, rate constants and even the identities of the essential elementary chemical processes are often unknown for chemical systems of interest. Attempts to model such systems should confine permissible steps to elementary processes or to certain specifically defined combinations of such processes. Permissible steps in a model can also be usefully restricted by recognizing that the direction of an elementary process is defined by strictly thermodynamic considerations even though thermodynamics cannot predict the direction of more complicated processes which may be coupled with other processes having free energy changes of opposite sign.
If the detailed state of a chemical system is defined at a specified time, and if the differential equations describing rates of change are known as unique functions of state, the subsequent behavior is calculable in principle. In practice, such calculations often encounter nontrivial computational difficulties. Transient intermediates at low concentrations are often both created and destroyed so rapidly the equations may become very “stiff’ unless increments of time are short compared to those appropriate to describe changes in major reactants. Excessive computational periods can be avoided only if the program can continuously adjust time increments to take advantage of opportunities to reduce computational effort. Fortunately such programs are now well established, and a hundred or more simultaneous equations can be handled even if they are quite stiff. If all chemical species, elementary processes, and rate constants can be assigned, the subsequent variation of a uniform system can be modeled with confidence. If concentrations can vary in space as well as in time, the situation is much more complicated. For the special cases of linear, cylindrical, and spherical symmetries, spatial variations can be described by a single position parameter. Especially if diffusion coefficients are independent of concentration, the simultaneous second-order equations describing such systems can be reduced to a number of first-order equations no more than twice what would be necessary if the system were uniform. If a system also contains thermal gradients, the modeling problems become much more severe. Thermal conductivity is described by second-order equations much like those for diffusion, but rate constants for chemical reactions are often strong functions of temperature. Thermal gradients often arise because the reaction of interest is strongly exothermic. It is obviously premature to claim that we could model each chemical system for which all rate constants, activation energies, heats of reaction, thermal conductivities, and diffusion coefficients are known, but such problems are being attacked with ever more sophistication. The problem facing the practicing chemist is still more severe because the necessary reaction parameters are often unknown for many of the important processes in a real system. It may not even be known which of the many conceivable component processes are significant. If the reactions themselves are associated with virtual discontinuities in rates such as occur in some chemical oscillators,
usual procedures for best fit of parameters are likely to be of little help. A particularly dramatic example was recently observed by Edelsonl during an effort to model the Sharma-Noyes2 mechanism of oscillations in iodine and iodide concentrations for the iodate catalyzed decomposition of hydrogen peroxide, For one set of rate constants assigned to about 20 potentially reversible reactions, the concentration of iodate never changed by more than a factor of 2. When two pairs of forward and reverse rate constants were changed in the same direction by a factor of 2 each, the concentration of iodate suddenly plunged by a factor of loll when a critical condition was attained during the reaction! When comparatively minor changes in rate constants can throw a system into an entirely different region of configuration space, techniques for fitting to experiment will be difficult a t best. Computer modeling certainly offers the best available method for testing the extent to which a complicated chemical mechanism is indeed understood. However, a moderately good fit to experimental data offers little confidence in the uniqueness of the model. We have found our confidence in modeling specific systems was enhanced if we restricted steps to elementary processes resembling established reaction types and if thermodynamic requirements were identified and respected. The application of these restrictions will be discussed briefly here. Restriction to Elementary Processes. It is generally believed that even the most complicated chemical change is the sum of the consequences of elementary processes each one of which involves no more than two or three individual molecules (with possible intervention of additional molecules of solvent). Moreover, most such processes can be assigned to a comparatively limited number of classes such as bond dissociations, recombinations, atom abstractions, proton or electron transfers, additions to or isomerizations of multiple bonds, etc. An ever growing body of data limits the ranges of plausible rate constants for many such reaction types. Reaction steps that involve the making and breaking of several chemical bonds at the same time cannot be elementary and should be replaced by sequences. The restriction to elementary processes can be relaxed for computations when equilibria between certain species are established very rapidly compared to the time scale of interest. Thus, in acidic aqueous solution it is usually permissible to assume that protons are transferred between The Journal of Physical Chemistry, Vol. 81, No. 25, 1977
2316
Richard M. Noyes
oxygens so rapidly that species differing only in degree of protonation can be assumed in equilibrium at all times. Each such assumption of equilibration should be carefully examined. It is also permissible to combine elementary processes that involve what Clarke3 has called a “flow-through reactant”. Such an intermediate never attains a significant concentration compared to any reactant from which it is formed, and the permissible processes for the intermediate are either reversion to starting materials or continuation to a single set of products. As an example, consider the sequence A t B2I I t C - P t Q
-
For this pair of elementary processes, I could be ignored and the single processes A + B + C P Q could be treated as elementary provided that I never attains a significant concentration compared to A, B, or C. If I can also undergo another process such as I + D M N, it is not a flow-through reactant and all steps must then be included in the modeling, Thermodynamic Restrictions. Our experiences with modeling various chemical oscillator^^^^^^^^ have greatly enhanced our appreciation of the applicability of thermodynamics to dynamic processes in systems far from equilibrium. For any closed chemical system at constant temperature and pressure, the total Gibbs free energy indubitably decreases monotonically with time. However, in a complicated system (such as a living organism and its immediate surroundings), certain chemical processes may go in a direction of increasing free energy provided they are coupled to other processes that involve an even greater decrease. It is not generaIly recognized that thermodynamic criteria again become valid when they are applied to elementary processes. By its very nature, a truly elementary process is not directly coupled to any other chemical process in the system. If chemical potentials of all reactant and product species are known, the direction of change is uniquely defined for an elementary process regardless of whatever else is happening in the system. These considerations lead us to recommend two procedures when complicated systems are modeled. First, assign standard free energies of formation to all reactant, product, and intermediate species included in the model. If such information is unavailable for some species, it should be estimated as well as possible. Flow-through reactants3 may be ignored in such a treatment, but all other intermediates must be included. These free energy data can often be used to derive electrochemical potentials6 for various conceivable oxidation-reduction processes, and many ostensibly plausible reactions can thereby be eliminated. We have found that this approach greatly narrows the range of permissible processes. Second, include the reverse of each elementary process. This simple procedure can protect a modeler from a lot of nonsense. T h e ratio of forward and reverse rate constants is always determined strictly by thermodynamics. If the model is not intended to cover the ultimate approach to equilibrium, this recommendation may be ignored for processes that are clearly virtually irreversible in the
+
The Journal of Physical Chemistry, Voi. 8 1, No. 25, 1977
-
+
situation of interest, but the burden of proof should be to justify exclusion of a reverse process rather than inclusion of on%. We must concede to not always having recognized the desirability of such a procedure. Only a few years ago, we published a complicated model5 containing three steps whose stoichiometric sum generated no net chemical change; we assigned the value of zero to one and only one of the six rate constants! Let that example be a warning of the pitfalls of modeling chemical reactions by computer. The above discussion does not pretend to define a ‘‘best” way to use a computer to model the chemistry of a complicated system for which many of the most desirable reaction parameters are experimentally inaccessible. It does suggest a few pitfalls to be avoided. Our understanding of chemical change must ultimately be able to lead back to a set of elementary processes describable by a set of dynamic parameters which are consistent for all the many systems in which those same elementary processes might occur. That long range goal is still very far from attainment.
Acknowledgment. The above ideas were developed during work on research problems supported in part by the National Science Foundation and by the Energy Research and Development Administration. Understanding has been sharpened by my interactions with Drs. Richard J. Field, David Edelson, Bruce L. Clarke, Kumud R. Sharma, Kenneth Showalter, Kedma Bar-Eli, Herman L. Liebhafsky, and John J. Tyson, but I take sole responsibility for the opinions expressed here.
References and Notes (1) D. Edelson and R. M. Noyes, to be published. (2) K.R. Sharma and R. M. Noyes, J. Am. Chem. Soc., 98,4345 (1976). (3) 6 . L. Clarke, J . Chem. Phys., 64, 4165 (1976). (4) R, J. Field, E. Koros, and R. M. Noyes, J. Am. Chem. Soc., 94, 8649 (1972). (5) D. Edelson, R. J. Field, and R. M. Noyes, Int. J . Chem. Kinet., 7, 417 (1975). (6) W. M: Latlkr, “Oxidation Potentials”, 2nd ed,Prentice-Hall, Englewood Cllffs, N.J., 1952.
Discussion W. A. REINHARDT (NASA, Ames Research Center). How important are convection and transport properties in the modeling of the complex mechanisms you have experimentally observed?
R. M. NOYES.The systems I discussed were deliberately stirred to maintain uniformity and to eliminate the need to complicate the modeling with transport effects. In an unstirred system, transport properties would need to be considered carefully in any effort to model experimental observations. P. HAAKE(Wesleyan University). I would like to see a brief presentation of the mathematics used for the computer simulation. A brief algorithm would be useful. In addition, if it does not take too much space, I think a typical application to a set of kinetic equations would be useful t o the readers. R. M. NOYES.I did not want to try to discuss computational procedures. Dr. Edelson developed the programs a t the Bell Laboratories, and he can provide appropriate references.
