Stoichiometry

Note: these classes are designed for internal use in class ReactionStoichManager. More...

Classes
class	C1
	Handles one species in a reaction. More...
class	C2
	Handles two species in a single reaction. More...
class	C3
	Handles three species in a reaction. More...
class	C_AnyN
	Handles any number of species in a reaction, including fractional stoichiometric coefficients, and arbitrary reaction orders. More...

Detailed Description

Note: these classes are designed for internal use in class ReactionStoichManager.

The classes defined here implement simple operations that are used by class ReactionStoichManager to compute things like rates of progress, species production rates, etc. In general, a reaction mechanism may involve many species and many reactions, but any given reaction typically only involves a few species as reactants, and a few as products. Therefore, the matrix of stoichiometric coefficients is very sparse. Not only is it sparse, but the non-zero matrix elements often have the value 1, and in many cases no more than three coefficients are non-zero for the reactants and/or the products.

For the present purposes, we will consider each direction of a reversible reaction to be a separate reaction. We often need to compute quantities that can formally be written as a matrix product of a stoichiometric coefficient matrix and a vector of reaction rates. For example, the species creation rates are given by

\[ \dot C_k = \sum_k \nu^{(p)}_{k,i} R_i \]

where \( \nu^{(p)_{k,i}} \) is the product-side stoichiometric coefficient of species k in reaction i. This could be done by straightforward matrix multiplication, but would be inefficient, since most of the matrix elements of \( \nu^{(p)}_{k,i} \) are zero. We could do better by using sparse-matrix algorithms to compute this product.

If the reactions are general ones, with non-integral stoichiometric coefficients, this is about as good as we can do. But we are particularly concerned here with the performance for very large reaction mechanisms, which are usually composed of elementary reactions, which have integral stoichiometric coefficients. Furthermore, very few elementary reactions involve more than 3 product or reactant molecules.

But we can do even better if we take account of the special structure of this matrix for elementary reactions involving three or fewer product molecules (or reactant molecules).

To take advantage of this structure, reactions are divided into four groups. These classes are designed to take advantage of this sparse structure when computing quantities that can be written as matrix multiplies.

They are designed to explicitly unroll loops over species or reactions for operations on reactions that require knowing the reaction stoichiometry.

This module consists of class StoichManager, and classes C1, C2, and C3. Classes C1, C2, and C3 handle operations involving one, two, or three species, respectively, in a reaction. Instances are instantiated with a reaction number, and n species numbers (n = 1 for C1, etc.). All three classes have the same interface.

These classes are designed for use by StoichManager, and the operations implemented are those needed to efficiently compute quantities such as rates of progress, species production rates, reaction thermochemistry, etc. The compiler will inline these methods into the body of the corresponding StoichManager method, and so there is no performance penalty (unless inlining is turned off).

To describe the methods, consider class C3 and suppose an instance is created with reaction number irxn and species numbers k0, k1, and k2.

• multiply(in, out) : out[irxn] is multiplied by in[k0] * in[k1] * in[k2]
• power(in, out) : out[irxn] is multiplied by (in[k0]^order0) * (in[k1]^order1) * (in[k2]^order2)
• incrementReaction(in, out) : out[irxn] is incremented by in[k0] + in[k1] + in[k2]
• decrementReaction(in, out) : out[irxn] is decremented by in[k0] + in[k1] + in[k2]
• incrementSpecies(in, out) : out[k0], out[k1], and out[k2] are all incremented by in[irxn]
• decrementSpecies(in, out) : out[k0], out[k1], and out[k2] are all decremented by in[irxn]

The function multiply() is usually used when evaluating the forward and reverse rates of progress of reactions. The rate constants are usually loaded into out[]. Then multiply() is called to add in the dependence of the species concentrations to yield a forward and reverse rop.

The function incrementSpecies() and its cousin decrementSpecies() is used to translate from rates of progress to species production rates. The vector in[] is preloaded with the rates of progress of all reactions. Then incrementSpecies() is called to increment the species production vector, out[], with the rates of progress.

The functions incrementReaction() and decrementReaction() are used to find the standard state equilibrium constant for a reaction. Here, output[] is a vector of length number of reactions, usually the standard Gibbs free energies of reaction, while input, usually the standard state Gibbs free energies of species, is a vector of length number of species.

Note the stoichiometric coefficient for a species in a reaction is handled by always assuming it is equal to one and then treating reactants and products for a reaction separately. Bimolecular reactions involving the identical species are treated as involving separate species.