Draft:Computing Chemical Kinetics

• Comment: This page reads like an essay on how to model chemical kinetics at a level appropriate for an undergraduate class. This is not what Wikipedia is for, please see WP:NOTGUIDE. There is already a decent Chemical kinetics page so the topic herein already has some coverage. At most a little of this (10%) might be added to that page, but really this draft is not appropriate for Wikipedia. Ldm1954 (talk) 03:36, 7 April 2024 (UTC)

---

Here, we explain the calculation of the concentration profiles of all reacting species for a given reaction mechanism, the underlying reaction rate constants and the initial concentrations of all reacting species. Often, this process is also called modelling.

The first step in any kinetic calculation or modelling is the definition of the ordinary differential equations (ODEs) which are defined by the mechanism or the chemical model. This is best seen in an example:

Consider the reversible reaction between the species A and B to form C. The reaction mechanism is given below:

$A+B\rightleftarrows C$ (k₊/k_-)

Usually, the differential equations for this mechanism are written as:

${d[A] \over dt}={d[B] \over dt}=-{d[C] \over dt}=-k_{+}[A][B]+k_{-}[C]$

The rate constants, k₊ and k_- define the changes of the concentrations of the interacting species as a function of the present concentrations, shown in square brackets, which of course change with reaction time.

Concentration profiles as a function of time for the equilibrium reaction.

For simple mechanisms explicit integration is sometimes possible, resulting in equations describing the concentration profiles for all reacting species as a function of time. However, as for the example above, explicit integration is not possible and numerical integration is required. Numerical integration is a well-developed branch of numerical mathematics, and we refer to Runge–Kutta methods and Stiff solvers.

The figure to the right shows the three concentration profiles for the reaction where the rate constants are k₊ = 10 M^-1s^-1, k_- = 1 s^-1 and the initial concentrations are [A]₀ = 0.04 M, [B]₀ = 0.1 M, and [C]₀ = 0.01 M.

The Matlab code for the generation of the concentration profiles and their plot is given below; the Matlab ode45 solver, based on the Runge-Kutta method, is used.

% A+B=C

c0=[0.04;0.1;0.01];                     % initial conc of A,B,C
k=[10;1];                           	% rate constants k+ and k-
[t,C]=ode45('ode_ApBeqC',2,c0,[],k);    % call ode-solver
plot(t,C)                               % plotting C vs t
xlabel('time');ylabel('conc.');

function c_dot=ode_ApBeqC(t,c,flag,k)

% A+B=C

c_dot(1,1)=-k(1)*c(1)*c(2)+k(2)*c(3);       % A_dot
c_dot(2,1)=-k(1)*c(1)*c(2)+k(2)*c(3);       % B_tot
c_dot(3,1)=+k(1)*c(1)*c(2)-k(2)*c(3);       % C_tot

One potential difficulty with this program is the coding of the differential equations, as seen in the function c_dot. It is easy to make a mistake and it may be difficult to realize that.

Fortunately. there is an ever-growing list of software which performs the task. Any search with the expressions 'chemical kinetics simulation' or 'chemical kinetics modelling' will result in a substantial list. Many of these packages are similar in their capabilities, however they vary significantly in their ease of use.

Despite the development of programs to support the analysis of complex chemical kinetics, there are two main issues that have 'troubled' kineticists forever: the problem of dealing with changes in activities and/or with changes in the pH over the course of the reaction. In fact, it is the lack of appropriate software that has forced kineticist to add excesses of inert salts and/or buffers to their reaction solutions in an attempt to control the effects these variables have on their reactions. It is important to note, however, that inert salts do not prevent the absolute changes in ionic strength, only the relative change, and buffers do not suppress pH changes completely, only reduce the severity of the pH fluctuations. Further, the addition of either inert salts or buffers will always result in interference to some extent with the reactions under consideration.

Modern software packages are able to deal with both conundrums and therefore remove the need for the addition of inert salts and/or buffers, hence making kinetic processes experimentally easier to investigate and the analysis more reliable and a truer representation of the chemical processes being studied.

Integration of activity based differential equations

The thermodynamically correct version for the differential equations is based on activities of the species to define the rates of concentration change, not the concentrations as was given in the commonly used equation above. Note: this principle is well-known in equilibrium studies, but has not been widely implemented in kinetics. The correct differential equations are:

${d[A] \over dt}={d[B] \over dt}=-{d[C] \over dt}=-k_{+}\{A\}\{B\}+k_{-}\{C\}$

Where we use the curly brackets, {A}, to represent the activity of the species A where {A} = $\gamma _{A}[A]$ and $\gamma _{A}$ is the activity coefficient of A.

Traditionally there were many reasons to use concentration based differential equations, foremost because the activities of species are difficult to determine and if the ionic strength of the reaction solution changes over time so do the activity coefficients and hence activities of the reacting species. This variation in activity coefficient value with changing ionic strength made calculations complex, beyond the capability of existing software. The standard and widely applied method to deal with this was to add a large excess of an inert salt to the reaction solution as this maintains an approximate constant ionic strength. As the activity coefficients, $\gamma _{A},\gamma _{B}$ and $\gamma _{C}$, are largely dependent on ionic strength. With an excess of inert salts, the activity coefficients are approximately constant and are ‘taken’ into the rate constants. Returning to the example:

${d[A] \over dt}={d[B] \over dt}=-{d[C] \over dt}=-k_{+}\gamma _{A}[A]\gamma _{B}[B]+k_{-}\gamma _{C}[C]$

$=-k_{+}^{'}[A][B]+k_{-}^{'}[C]$

$k_{+}^{'}=-k_{+}\gamma _{A}\gamma _{B}\quad and\quad k_{-}^{'}\gamma _{C}$

Consequently, the experimental rate constants, k₊’ and k_-‘, are reported at a particular ionic strength.

The thermodynamic rate constants, which by definition are at zero ionic strength, can be determined by determining multiple k₊’ and k_-‘ values at different ionic strengths and extrapolating these back to zero. The dependence of the observed kinetics on ionic strength is also known as the kinetic salt effect . This process is very expensive in time as well as reagents as multiple experiments are required for an acceptable extrapolation to zero ionic strength.

Accounting for changes in ionic strength

An alternative approach to the addition of excess inert salts is to estimate the activity coefficients by calculating the ionic strength at any time during the reaction and use the estimated activity coefficients to compute the differentials. The Debye-Hückel approximation and its derivatives are relatively easily computed and can be used in the numerical integrations. At any time during the numerical integration the concentrations of all reacting species are known, in this example they are [A]_t, [B]_t, and [C]_t. Using their known charges, the ionic strength, I, is defined and allows the approximate calculation of the activity coefficients and thus the activities of all species at any time interval in the reaction. It is worth stressing that these activities are an approximation and generally, the higher the ionic strength the more imprecise the approximation. However, compared with the traditional approach which assumes that all activity coefficients are constant, this approach is superior as it takes into account the changes in ionic strength.

Computing the changing ionic strength and subsequent approximations of the activity coefficients does mean that the numerical integration is more complex. Fortunately, there is freely available software which performs such calculations, ReactLab KINSIM, which was used to simulate all the reactions below.

Concentration profiles as a function of time with different sets of charges for the species.

To demonstrate the effect that the inclusion of activity coefficients and variation in ionic strength throughout a kinetic reaction has, we have considered the above reaction but varied the charges on the species in the reaction: Set 1: A^2-, B⁰, C^2-(dotted lines); Set 2: A^2-, B²⁺, C⁰ (dashed lines); and Set 3: A⁰, B⁰, C⁰ (full lines).

The Figure to the right summarizes the results. The plot includes the concentration profiles for the species A (light blue), B (orange) and C (dark blue) as well as the ionic strengths for the Sets 2, and 3 (in black, using the secondary y-axis on the right). The Debye-Hückel approximation was used to approximate the activity coefficients for all sets of data.

The following observations can be made:

The solid lines, Set 3, are the same as in the first Figure above as species with a zero charge by definition have an activity coefficient of one so using either concentration or activity based differential equations result in the same concentrations.

For Set 1, the ionic strength is constant as the 2- species A is transformed into the 2- species C, meaning there is no overall change in charge as the reaction progresses. So for the given initial concentrations the ionic strength is a constant value of 0.15 (dotted black on the secondary y-axis). Even so, the reactions for this case are significantly slower compared with uncharged species as the activity coefficients for the charged species are much less than one, demonstrating the real effect that charged species can have on the rate of a reaction.

Set 2 has astonishing effects as now the reaction goes ‘the other way’. The equilibrium at the end of the reaction is on the side of the ‘starting’ materials A and B and the concentration of C gets lower with the reaction. The reason is that the activities of the ions are much lower than their concentrations while the activity of the uncharged C equals the concentration.

The ionic strength in both Sets does not change much as the concentrations of the counter ions to the charged species do not change with time.

Taking pH changes into account

Many reactions in aqueous solution include Lewis acids and/or Lewis bases and consequently they include protonation equilibria which are affected by pH and often influence the concentrations of the reactants. To control such interferences buffers have traditionally been added; they approximately maintain pH and thus the relative positions of the protonation equilibria. Buffers have similarities with inert salts in that they need to be added to maintain some potentially interfering properties of the solutions, in this case the pH. Again, these additions can be avoided if the changes are simply taken into the computations. This is again best illustrated with an example, this time a real example, the complexation reaction of cyclam with Cu²⁺ .

The ligand cyclam.

Cyclam features four secondary amines in a 14-membered macrocycle, see the figure to the right. All amine groups are involved in protonation equilibria, defined by four protonation constants (K₁ to K₄ in the scheme below). Under the conditions of the experiment only the species LH and LH₂ react significantly with the metal ion. LH₃ is not sufficiently reactive while LH₄ is not reactive at all. Further, the highly reactive unprotonated state of the ligand, L, exists only in very low concentrations. At the pH of the solution, the formation reactions with Cu²⁺ are essentially irreversible. The complete reaction scheme thus is:

$\qquad L+H^{+}\longleftrightarrow LH^{+}\qquad \qquad \qquad logK_{1}=11.6$

$LH^{+}+H^{+}\longleftrightarrow LH_{2}^{2+}\qquad \qquad \qquad logK_{2}=10.6$

$LH_{2}^{2+}+H^{+}\longleftrightarrow LH_{3}^{3+}\qquad \qquad \qquad logK_{3}=1.56$

$LH_{3}^{3+}+H^{+}\longleftrightarrow LH_{4}^{4+}\qquad \qquad \qquad logK_{4}=2.34$

$Cu^{2+}+LH^{+}\longrightarrow CuL^{2+}+H^{+}\qquad \qquad k_{LH}=1.1\times 10^{6}M^{-1}s^{-1}$

$Cu^{2+}+LH_{2}^{2+}\longrightarrow CuL^{2+}+2H^{+}\qquad k_{LH_{2}}=0.14\ M^{-1}s^{-1}$

The differential equation describing the kinetics is:

${d[Cu^{2+}] \over dt}={d[L]_{tot} \over dt}=-{d[CuL^{2+}] \over dt}=-k_{LH}\{Cu^{2+}\}\{LH^{+}\}-k_{LH2}\{Cu^{2+}\}\{LH_{2}^{2+}\}$

Cyclam concentration profiles complexing with Cu(II)

With the additional constraints at any time (x = 1 - 4):

$K_{x}={\frac {\{LH_{x}\}}{\{LH_{x-1}\}\{H^{+}\}}}$

[L]_tot is the sum over the concentration of all differently protonated forms of L:

$[L]_{tot}=[L]+[LH^{+}]+[LH_{2}^{2+}]+[LH_{3}^{3+}]+[LH_{4}^{4+}]$

In a real experiment, where the initial concentrations used were [CuCl₂] = 0.00394 M, [cyclam] = 0.00433 M and [HCl] = 0.00960 M, and the rate and equilibrium constants are as given in the above reaction scheme, the following concentration profiles are calculated. Here, for simplicity, all activity coefficients were set to a value of one.

Time plot of the changing activity coefficients over the course of the reaction.

We observe the general decrease of the concentration of the free metal and the combined ligand concentrations; however, note that the [LH₄⁴⁺] increases with time as the pH decreases (black line on the secondary y-axis).

It is also instructive to plot the time dependence of the activities (using the Davies equation) as well as the ionic strength of the solution. The time behaviour is similar to the concentration profiles but the activities are considerably lower as the activity coefficients are significantly lower that one even at the moderate ionic strength of about 0.027 M. The activity coefficient for highest charged LH₄⁴⁺ is approximately 0.08. The ionic strength does not change much as only cations react with each other and thus the overall charges do not change much. The concentrations of the anions, the counterions of Cu²⁺ and H⁺, do not change at all. More information on the Cu²⁺/cyclam system is available in the freely downloadable manual for ReactLab Kinetics PRO.

Both example demonstrates how incorporating the pH and ionic strength changes into the calculations when analyzing kinetic data enables the scientist to more robustly analyze the system being investigated while at the same time simplifying the experimental procedure.