Chemical kinetics of reactions

The study of dynamic processes in chemistry on a macroscopic scale often refers to chemical kinetics.

This discipline concerns the patterns of chemical reactions over time, their dependence on external conditions and the mechanisms of chemical transformations.

Let's implement some basic examples from this field. Since the mathematical basis for studying dynamic processes is differential equations, we will need several basic functions that are in the library. DifferentialEquations.

Pkg.add("DifferentialEquations")

As part of a small project to study simple models of chemical kinetics, we will go through four examples:

Abstract chemical reaction for three substances as a model DifferentialEquations,
Switch to the library Catalyst to simplify the syntax of the model,
Transition to realistic names of substances (hydrogen sulfide oxidation),
Adding equations to the model: accounting for the effect of temperature on the reaction rate (the Arrhenius model).

Reactions as differential equations

In the first example, we calculate a sequential chemical reaction, where the substance A turns into B, and then B turns into C. The process is described by a system of ordinary differential equations (ODES), where the rate of change in the concentration of each substance depends on the reaction constants and . The equations look like this:

(substance A is consumed),
(substance B is formed from A, but decomposes into C),
(substance C accumulates).

This model allows us to predict how the concentrations of all participants in the reaction change over time.

using DifferentialEquations

# Определяем систему ОДУ: dA/dt = -k1*A, dB/dt = k1*A - k2*B, dC/dt = k2*B
function kinetics!(du, u, p, t)
    A, B, C = u
    k1, k2 = p
    du[1] = -k1 * A          # dA/dt
    du[2] = k1 * A - k2 * B  # dB/dt
    du[3] = k2 * B           # dC/dt
end

# Начальные условия и параметры
u0 = [1.0, 0.0, 0.0]  # [A0, B0, C0]
p = [0.5, 0.2]        # константы скорости k1, k2
tspan = (0.0, 10.0)   # интервал времени

# Решаем систему
prob = ODEProblem(kinetics!, u0, tspan, p)
sol = DifferentialEquations.solve(prob, Tsit5())

# Визуализация
plot(sol, label=["A" "B" "C"], xlabel="Время", ylabel="Концентрация", lw=2)

Such graphs are observed when calculating a system with parameters = 0.5, = 0.2 and at the initial concentration ( = 1, = 0, = 0 ). The standard numerical integration method was used to solve the problem. Tsit5().

We have shown how to solve a typical example for lectures on chemical kinetics of medium complexity. Now let's explore how to simplify its description and make the entire methodology more suitable for solving realistic tasks. Let's describe this reaction using the Domain Specific Language (DSL) commands provided in the library. Catalyst.

Using the Catalyst Library

Let's set the same consistent reaction using a shorter description given by the commands from the library Catalyst.jl (by installing this library first).

Pkg.add("Catalyst")

The library contains many tools that simplify the solution of chemical kinetics problems. For example, a symbolic API for creating reaction chains (@variables, @species, Reaction, ReactionSystem), tools for analyzing reaction networks (connectivity classes, scarcity, reversibility, detection of conservation laws) and visualization of networks through GraphMakie and CairoMakie, conversion of reaction systems into ODES, stochastic dual or hopping processes, integration with packages for calculations on the GPU, etc.

But the main technique that we will use is the Catalyst DSL command (a language for describing reactions), the basic block is the construction @reaction_network begin .. end.

rn = @reaction_network begin
    k1, A --> B
    k2, B --> C
end

Each reaction is set in a separate line inside the block. begin..end. The first element of the description is the rate of the reaction, then the ratio of the chemical symbols included in the reaction. Among the technical features, we will need to set the following additional parameters:

Initial conditions,
number of steps (limited time vector),
Parameter values.

The resulting system can be designed as an ODE (object ODEProblem), calculate using the command solve and display it on a graph using the command plot.

using Catalyst, DifferentialEquations
using Catalyst: solve # Неважно, из какой библиотеки брать эту функцию

# Определяем реакционную сеть
rn = @reaction_network begin
    k1, A --> B
    k2, B --> C
end

# Задаем параметры и начальные условия
p = [:k1 => 0.5, :k2 => 0.2]           # константы скорости
u0 = [:A => 1.0, :B => 0.0, :C => 0.0] # начальные концентрации
tspan = (0.0, 10.0)                    # интервал времени

# Создаем и решаем задачу ОДУ
prob = ODEProblem(rn, u0, tspan, p)
sol = solve(prob, Tsit5())

# Визуализация
plot(sol, label=["A" "B" "C"], xlabel="Время", ylabel="Концентрация", lw=2)

We see the familiar graph that we have obtained for the same initial conditions and reaction rates as in the first case.

Realistic variable names

Let's take a half step towards a more convenient description of the reactions. To portray a realistic reaction, you just need to change the names of the variables.

Consider the oxidation reaction of hydrogen sulfide ( ) first to sulfur (), then to sulfur dioxide () in the presence of oxygen. Simplified: ( ). This is a training example where all the stages have first-order kinetics.

using Catalyst, DifferentialEquations
using Catalyst: solve

# Определяем реакционную сеть
rn = @reaction_network begin
    k1, H2S --> S      # H₂S → S
    k2, S --> SO2      # S → SO₂
end

# Параметры и начальные условия
p = [:k1 => 0.3, :k2 => 0.1]               # константы скорости (условные, 1/с)
u0 = [:H2S => 1.0, :S => 0.0, :SO2 => 0.0] # начальные концентрации (моль/л)
tspan = (0.0, 20.0)                        # интервал времени (с)

# Решаем задачу ОДУ
prob = ODEProblem(rn, u0, tspan, p)
sol = solve(prob, Tsit5())

# Визуализация
plot(sol, label=["H₂S" "S" "SO₂"], xlabel="Время (с)", ylabel="Концентрация (моль/л)", lw=2 )

This is a simplified model, but it reflects real chemistry, for example, the processes of purifying gases from hydrogen sulfide.

We supplement the description of the reaction with the Arrhenius model

Finally, let's look at an example that shows how to complement the oxidation example. by adding additional dependencies to it, namely the dependence of the reaction rate on temperature.

To a first approximation, chemical reactions between substances occur only as a result of the collision of molecules of these substances (simple collision model). To overcome the energy barrier, molecules must have a certain minimum activation energy, hence the model , where characterizes the frequency of collisions, – temperature, – universal gas constant, – activation energy. The number of molecules involved in a chemical reaction is determined from the Boltzmann distribution as proportional .

Now our reaction has a realistic flow rate, significantly higher than the speed in the previous example. The new observation interval is reflected in the parametertspan.

using Catalyst, DifferentialEquations
using Catalyst: solve
using Plots.PlotMeasures

# Функции для констант скорости (например, уравнение Аррениуса)
k_1(T) = 1e5 * exp(-5000 / (8.314 * T))  # A * exp(-Ea/RT), условные параметры
k_2(T) = 2e4 * exp(-3000 / (8.314 * T))  # разные A и Ea

# Реакционная сеть с функциональными константами
rn = @reaction_network begin
    k_1(T), H2S --> S    # k1 зависит от температуры T
    k_2(T), S --> SO2    # k2 тоже зависит от T
end

# Параметры и начальные условия
tspan = (0.0, 0.0002)
p = [:T => 273.0]
u0 = [:H2S => 1.0, :S => 0.0, :SO2 => 0.0]

# Решаем ОДУ
prob1 = ODEProblem(rn, u0, tspan, p)
sol1 = solve(prob1, Tsit5())

# График
p1 = plot(sol1, label=["H₂S" "S" "SO₂"], xlabel="Время (с)", ylabel="Концентрация (моль/л)", lw=2,
          title="Реакция окисления H₂S при 273 K", titlefont=font(11), guidefont=font(8),
          bottommargin=30px)

# Параметры и начальные условия
p = [:T => 400.0]
u0 = [:H2S => 1.0, :S => 0.0, :SO2 => 0.0]
prob2 = ODEProblem(rn, u0, tspan, p)
sol2 = solve(prob2, Tsit5())

p2 = plot(sol2, label=[false false false], xlabel="Время (с)", ylabel="Концентрация (моль/л)", lw=2,
          title="Реакция окисления H₂S при 400 K", titlefont=font(11), guidefont=font(8))

plot(p1, p2, layout=(2,1))

Conclusion

These examples demonstrate two approaches to modeling chemical kinetics in Engee using the Julia language. The first option, manually setting the system of differential equations, gives you full control over the calculations and is useful for understanding the basics. The second method, using Catalyst.jl, allows you to describe reactions in a concise format, automatically generating equations. This is especially useful for complex systems where development speed and error minimization are important.

In addition, we have studied how to add additional equations to the model and make the entire methodology more suitable for solving real-world problems of chemical kinetics.