I'd like to programm a toy chemical reaction program. Rough simplifications that deviate from industry and academia are totally acceptable. Not very educationaly accurate, but fun.

Are there methods to calculate rate of reaction for dozens and hundreds of reactants in a closed system? Or I can blatantly use original r = k(T)[A][B] for every pair of reactants while selectively ignoring presence of other molecules?

Or should I use something else to evaluate what will happen with hundred of reactants after two seconds?

Sample gas composition below: enter image description here

  • 4
    $\begingroup$ As a programmer, be aware that you're running straight into a combinatorial explosion of epic magnitude. Dozens of reactants means hundreds of possible pairs, many of which will not react, but many will. (How are you going to find that out, to begin with?) Some reactions may yield different products in different conditions. Then the products might want to react with something else. Then some of them might form another phase and bring an entirely new level of complexity (which is the case in your example, unless T is pretty high, and if it is, that would cause other issues). Then... $\endgroup$ Dec 8, 2017 at 9:40
  • 1
    $\begingroup$ Start with a much simple system, I say! Then see how you get on. $\endgroup$
    – Gert
    Dec 8, 2017 at 13:34
  • 1
    $\begingroup$ You may be interested in the program Cantera, which does precisely what you seek to accomplish, provided you have estimates for the reaction barriers, rate constants, and thermochemistry. Looking up any examples on microkinetic modeling will lead you to plenty of information. Another example, but in the context of heterogeneous catalysis, is CatMAP. The corresponding paper is a fairly nice reference too. $\endgroup$
    – Argon
    Dec 9, 2017 at 5:48

3 Answers 3


The questions has at least two aspects to address.

The first thing is, there is no guarantee that the reaction is first-order with respect to each of the reactants in first place. In other words, without further investigation you can not just take $r=k[\ce{A}][\ce{B}]$ because $r=k[\ce{A}][\ce{B}]^2$ or $r=[\ce{A}]^2[\ce{B}]$ or something more exotic may happen, depending on the nature of the reaction.

OK. So let's say for simplicity we take $r=k[\ce{A}][\ce{B}]$ when the system has only $\ce{A}$ and $\ce{B}$, as we don't need to be quite educationally correct. Then yes, after mixing many reactants, the equation still holds for each pair of the reactants, provided they are elementary reactions. It will be a completely different story if, for example, neither of two of the reactions are elementary and they share the same reaction intermediates.

I also have a reminder that when stepping from the expression of the rate to the expression of concentration with respect to time, you have to be careful and do the integration properly instead of just naively extrapolating the concentration expression when only two reactants exists. Let's say that we don't consider chemical equilibrium for the moment (i.e. consider all the reverse reactions to be prohibitively difficult to happen), then as far as I can see, the $r-t$ curve will not be smooth throughout the time and will have some turning point every time a reactant is fully consumed.

Side note: I suggest that you do your formatting properly in your program for the superscripts and subscripts.

  • $\begingroup$ In general, assuming that reactions are at least first order in all their reactants, the concentrations will be smooth functions of time. At least if you zoom in close enough. One can have pretty sharp-looking corners e.g. if a catalytic decomposition reaction is normally rate-limited by the saturation of the catalyst, but at some point before the reactant concentration would hit zero, the catalyst will stop being saturated and the reaction rate will become proportional to the reactant concentration. $\endgroup$ Dec 8, 2017 at 17:21
  • $\begingroup$ ... That said, I certainly agree with you that one has to do proper numerical integration, rather than just linearly extrapolating the concentrations. In fact, for these kinds of problems, usually one should use an implicit solution method that avoids numerical instability caused by extrapolation within a timestep. $\endgroup$ Dec 8, 2017 at 17:26

For truly elementary reactions, the law of mass action can generally be assumed to hold fairly accurately, as predicted by the collisional theory of reaction kinetics. (This is something of a tautology, since a deviation from the law of mass action generally indicates that the reaction isn't really elementary.) In particular, this implies that:

  • elementary unimolecular reactions occur at a rate proportional to the concentration of the reactant (or more generally to its activity, if we're treating e.g. the effects of a solvent on the reaction rate implicitly);
  • elementary bimolecular reactions with distinct reactants occur at a rate proportional to the product of the concentrations (or activities) of the reactants;
  • elementary reactions between two molecules of the same species occur at a rate proportional to the square of the concentration (or activity) of the reactant; and
  • elementary reactions involving more than two reactants don't occur, except at negligible rates, since they would have to involve a perfectly simultaneous collision of three molecules. (Any such reactions that are observed to occur must actually proceed via one or more ephemeral transition states produced in bimolecular collisions, and are thus not actually elementary.)

Thus, in principle, any* such system of simultaneous elementary reactions can be modeled using a system of ordinary differential equations of the form:

$$ \frac{d}{dt} [A_x] = \left( \sum_{y=1}^n k_{xy}[A_y] \right) + \left( \sum_{y=1}^n \sum_{z=1}^y k_{xyz}[A_y][A_z] \right), $$

where $[A_1]$ to $[A_n]$ are the concentrations of the $n$ distinct chemical species present in the system (some of which might only exist as temporary transition states), $k_{xy}$ is the sum of the rate constants for reactions of the form $A_y \to A_x + (\text{other products})$, and $k_{xyz} = k_{xzy}$ is the sum of the rate constants for reactions of the form $A_y + A_z \to A_x + (\text{other products})$.**

So yes, if you know the rate constants for all the elementary mono- and bimolecular reactions your reactants can undergo, then you can indeed plug them into a differential equation like the one above and solve it.

However, as a fair warning from someone who's worked on numerically solving these kinds of differential equation systems, they are notorious for often exhibiting high stiffness and consequent numerical instability. This is particularly often the case if some of the reaction constants are much larger and/or some of the reactant concentrations are much higher than others. Efficiently solving such stiff ODE systems is a hard problem, and there's a whole bunch of advanced (and often jargon-filled) mathematical literature on it.

In the end, the practical ways of solving such systems usually come down to timescale separation or other stiffness-reducing approximations, where we essentially assume that certain "fast" reactions occur so rapidly that the concentrations of their products and reactants are always at a local quasi-equilibrium, where the overall net reaction rate of the fast reactions and their reverse reactions is zero. This allows us to then only focus on the remaining "slow" reactions, which are presumably the ones we're actually interested in.

Alas, as a side effect, these approximations also often tend to destroy the mathematical simplicity of the basic elementary reaction model by reintroducing non-elementary reaction laws and/or nontrivial dependencies between the reactant concentrations. Also, unfortunately, at least I know of no general method for automatically constructing such approximations, just as there's no fully general method for efficiently solving the original stiff ODE systems in the first place.

*) At least assuming that the reactions occur in a homogeneous well-mixed environment. Things like reactions occurring on the surface of a substrate can be included by treating the substrate-bound and free forms of a molecule as distinct species and, if necessary, also treating free surface bonding sites as a separate species with its own concentration. Similar tricks can also be used to model e.g. reactions occurring in two immiscible solvents, provided that mixing within each solvent phase can be assumed to be rapid. However, if the reactant concentrations may vary spatially to a significant degree, then one would have to use a reaction–diffusion model instead.

**) If the products of a reaction include $m$ molecules of species $A_x$, that reaction's rate constant should be included $m$ times in $k_{xy}$ or $k_{xyz}$. Also note that, to avoid double counting $A_y + A_z \to A_x + \dots$ and $A_z + A_y \to A_x + \dots$ as distinct reactions, the last sum in the differential equation is only taken over $z \le y$.


Allright, guys. All I wanted happens to be stochastic simulation.

Here is the list of sample algorithms:

  • Gillespie SSA - original.
  • Tau-leaping - evaluate time step and then out of that reactions and substances amount.
  • R-leaping - vise versa. Good for big systems with only a few relatively highly reactive species.
  • S-leaping - hydrid of R-leaping and Tau-leaping.

Tau-leaping can be explicit and implicit. The second one is better for stiff systems. Adaptive Tau-leaping switches between one and the other.

Many algorithms for Tau-leaping time step evaluation exist:

  • The first original algorithm to evaluate time step in Tau-leaping is flawed in performance - Cao, Y.; Gillespie, D. T.; Petzold, L. R. (2005). "Avoiding negative populations in explicit Poisson tau-leaping"

  • Upgraded one is better all around - Cao, Y.; Gillespie, D. T.; Petzold, L. R. (2006). "Efficient step size selection for the tau-leaping simulation method"

  • Don't bother with Binomial distribution. No gain over Cao-Gillespie-Petzold-2006 upgrade but greater algorithm complexity.

  • There is also Hybrid Chernoff Tau-leaping paper

  • And Anderson, David F - 2005 and 2018 papers

Or should I use something else to evaluate what will happen with hundred of reactants after two seconds?

The great mistery is why none of you have mentioned a single algorithm above.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.