Autores: A partir do fascículo 39/9 a revista Química Nova adotou a licença CC-BY. Mais informações a respeito dessa licença podem ser obtidas aqui.

An approach to the kinetics and thermodynamics of elementary chemical reactions using a stochastic model |

Francis Pereira Nascimento; Baraquízio Braga do Nascimento Junior; Luiz Augusto Martins Cardoso; Rodrigo Veiga Tenório de Albuquerque; Nemesio Matos Oliveira-Neto Departamento de Ciências e Tecnologias (DCT), Universidade Estadual do Sudoeste da Bahia, 45208-409 Jequié - BA, Brasil Recebido em: 16/11/2017 Endereço para correspondência
*e-mail: nmon@uesb.edu.br RESUMO
A mathematical model is proposed to investigate the kinetics and equilibrium of homogeneous elementary first- and second-order chemical reactions. The dynamics is defined by the Monte Carlo method (MCM), with the Metropolis update, and the Ehrenfest urn model (EUM). MCM is an important step that accesses the kinetic and thermodynamic properties of the system, while the EUM defines the orders of the reactions studied in this work. The main parameters, such as temperature, the activation energy and the steric factor, were taken into account in the calculation of the transition probabilities between the reactants and products. It is thus possible to reproduce the kinetic profiles of the reactions and to evaluate the influence of temperature and the steric factor. Furthermore, it is possible to simulate the behaviour of the system by modifying the activation energy barrier, thus simulating a catalytic process. The effect of the addition of molecules was also investigated, with the system returning to a thermodynamic equilibrium condition after partial consumption of the added reactant, as predicted by Le Chatelier's principle. All the simulated data are in agreement with the theoretical results present in physical chemistry textbooks. Palavras-chave: chemical reaction; kinetics; thermodynamics; Monte Carlo method; Ehrenfest urn model.
Since John Dewy recommended the inclusion of inquiry in the K-12 science curriculum in 1910, science education has undergone several modifications and adaptations around the world. Stochastic methods offer a number of advantages. First, due to the lack of infrastructure at some school laboratories and universities, certain important topics in physical chemistry are either neglected or only treated theoretically. Another advantage of such numerical approaches is that it is not always possible to achieve the necessary high pressure and temperature regimes of some chemical reactions in the laboratory. Also, from a theoretical point of view, it is not always possible to analytically solve the system of ordinary differential equations that describe a chemical reaction due to its complexity. A simple stochastic approach based on MCM and the Ehrenfest urn model (EUM)
The standard MCM, where ^{-1}, k is the Boltzmann constant, _{B}T is the temperature, p and _{f}^{s}p are the steric factors for the forward and the reverse reaction, respectively, and _{r}^{s}E and _{a}E are the activation energies for the forward and the reverse steps, respectively. For an irreversible reaction, Equation (2) is not taken into account, as _{b}P = 0, which is equivalent to assuming _{r}E → ∞. This part of the present approach, defined by Equations (1) and (2), is quite similar to that used in the hydrogen atom ionisation problem._{B}^{20,21} Since the transition probabilities are well known, the important features that are intimately related to the rate law of elementary reactions can be determined, namely the Arrhenius' Law of kinetics.The dimensionless temperature-related variable, θ ≡ E, and dimensionless energy-related variable, _{a}E ≡ _{ba}E/_{b} E, are defined to simplify the calculations. The computational simulations carried out in this work are thus based on the thermal energy _{a}k and _{B}TE both of which are measured with respect to _{b},E._{a}Equations (1) and (2) define one part of the chosen applied dynamics of the elementary chemical reaction. The other part of the dynamics is defined via EUM as follows. N molecules corresponding to the type-A and type-B molecules, respectively. If a given box _{B}A is occupied with a type-A molecule, the corresponding _{i}B box will be empty and vice versa. In addition to the definition given in the literature,_{i}^{10-12} the relationship between the total number of molecules is defined here as N ≡ N + _{A}N, and the total number of available boxes _{B}M is limited by N, such that N ≤ M. The effects of the number of molecules on the model can be investigated when N = N + _{A}N < _{B}M, for instance, by adding molecules during the simulation, with the case where N = M corresponding to the maximum capacity (full-filled case) of the containers (urns).The following steps were carried out to complete the dynamics of a first-order reversible reaction: (I) Try to select a type-A molecule by choosing a uniform random number integer z ∈ [0,1] is generated and compared to P from Equation (1). If _{f}z ≤ P, the molecule indexed as _{f}i is eliminated from box A, and a type-B molecule is created in box _{i}B, which indicates that the reaction occurred, transforming A into B. If _{i}z > P, the type-A molecule remains in its original form._{f}(II) Try to select a type-B molecule, following the same approach employed for the type-A molecule in step (I). If the choice fails, the step ends. If the choice succeeds, a uniform random number z ≤ P, the reaction occurs and B is transformed into A. If _{r}z > P, the type-B molecule remains in its original form._{r}By running steps (I) and (II) For an irreversible second-order reaction, A + B → C, only Equation (1) is followed, thus representing the transition probability of the type-A and type-B molecules to overcome the energy barrier, (III) Try to select two molecules, one of each type (A and B). If at least one of them fails to be chosen during the selection, the step ends. If the choices succeed, a uniform random number z ≤ P the reaction occurs and A reacts with B to form C. If _{f}z > P, the A and B molecules remain unreacted. As previously mentioned above, by performing step (III) _{f}M times, one MCS is completed.It should be noted that the simulation started with given molar fractions of the type-A and type-B molecules, χ _{B}_{,} respectively. For an arbitrary choice of θ, the dynamics of the reaction and the influence of model parameters can be studied. Unless otherwise stated, all simulations will be performed with N = 10^{6} molecules. The influence of N will be treated at the end. The time is measured in Monte Carlo Steps (MCSs).In the next section, the irreversible reaction, A → B, will be studied, followed by the reversible A Breaction and finally the second-order A + B → C transformation.
The irreversible first-order reaction A → B, will be studied first. Figure 1 shows χ as a function of time, in MCSs. The simulation started with all type-A molecules, where χ_{B}^{0}_{A}= 1.0, χ^{0}_{B} = 0,0, θ = 1.0 and steric factor p = _{f}^{s}p = 1.0. As can be seen, χ^{s} decreases until these species is completely consumed, whereas χ_{A} approaches unity. The inset figure shows the behaviour of χ_{B} on a log scale, showing the exponential time dependence. This result fits well to_{A}
Figure 1. χ_{A} and χ_{B} as a function of time, in MCSs
where The effect of temperature on the rate constant is shown in Figure 2, where χ k/_{B}TE, this behaviour can be interpreted in two different ways._{a}^{11,12} One is that the increase in θ means that the k term has to increase when using a constant activation energy. We thus explore the effect of temperature on a given system by defining a specific activation energy. However, the increase in θ can also occur by fixing _{B}Tk and decreasing _{B}TE. Therefore, we also investigated different systems (different _{a}E) with the same temperature. Considering the set of parameters (θ values) used in Figure 2, the same relationship between the rate constant (_{a}k) and temperature, in the form of the Arrhenius' Law, is observed
Figure 2. χ_{A} evolution as a function of time for different temperatures
where = 1.002 ± 0.003. This is an expected result, because the Metropolis update used to simulate the data has the same mathematical structure as the Arrhenius' Equation, which can be seen in the comparison between Equations (1) and (4). Another important feature here is the bhalf-life (τ_{1/2}). Using this model, the simulation yielded^{11,12}where The influence of the steric factor in the system was also addressed. This parameter was introduced to take into account the dependence of the rate constant on the mutual orientation of the reactant molecules. The effect of the steric factor on the reaction is shown in Figure 3, where χ p decreases ^{s}k, which leads to a slower consumption rate, thus highlighting the fact that the molecules are not able to react in all directions. This implies that there is a favourable collision orientation, from a geometric point of view, that becomes increasingly limited. Here the case where p = 1.0 means that all collisions are geometrically favourable.^{s}
Figure 3. The effect of the steric factor on the irreversible first-order reaction
An important issue that deserves attention is the fact that the steric factor does not change the main aspects of the kinetics. A simulation using The rate constant ^{/θ}, with A = 0.49881 and b' = 0.9992 ± 0.0004. It is thus worth noting that, within the standard error limit, follows the Arrhenius' Law, with the only change being the pre-exponential factor A, as demonstrated in Equation (4). For other p values, it is seen that ^{s}A ~ p, thus indicating that this approach is useful in dealing with the kinetics of irreversible chemical reactions. As mentioned above, it is expected that ^{s}A = p and ^{s}b¢ = 1, which represent the parameters and in Equation (4), since the transition probabilities have been defined from a mathematical structure similar to the Arrhenius' Equation.
Figure 4. (a) The rate constant k, on a log scale, as a function of the reciprocal of θ and (b) The half-life τ_{1/2} as a function of 1/k
The half-life τ _{1/2} on 1/k is approximately linear, as indicated by the solid line on Figure 4(b). The fitted data yield τ_{1/2} ~ c'/k, with c' = 0.699 ± 0.001, which is in good agreement with Equation (5) which was obtained without steric factor.This approach is applied in the kinetics of catalytic reactions by considering an energy change in the following manner: E/E. The kinetic profiles of the catalytic and non-catalytic reactions, with θ = 0.9 and no steric factor, are shown in Figure 5. As can be seen for a fixed value of θ, as the energy decreases (Δ_{a}E < 0), a given type-A molecule has greater chance (probability) of overcoming the potential energy and being transformed into a type-B molecule. This analysis thus shows that by decreasing the energy barrier, the rate constant of the catalytic reaction increases in comparison to the non-catalytic one. An analogue analysis can be made for Δ_{cat}E > 0._{cat}
Figure 5. The kinetic profiles of the catalytic and non-catalytic reactions
Here the reversible reactions of the type A B are discussed, where both Equations (1) and (2) are considered. For simplicity, the steric factors are set to p = 1. Unlike the irreversible reactions, it is necessary to define the dimensionless energy parameter, _{r}^{s}E ≡ _{ba}E/_{b}E, which represents the reverse reaction energy barrier (rationed) to _{a}E. Some results are shown in Figure 6, where _{a}E = 2.0 are considered, with χ_{ba}^{0}_{A} = 0.8 and χ^{0}_{B} = 0.0. The behaviours of χ and χ_{A} are shown in Figure 6(a) for two distinct temperatures, θ = 1.0 (empty squares and circles) and θ = 1.5 (filled squares and circles)._{B}
Figure 6. (a) The time evolution of χ_{A} and χ_{B} for two distinct temperatures and (b) The influence of the addition of the type-A reactant on the kinetics of the reaction
When the system reaches thermodynamic equilibrium, χ are altered by increasing θ. The increase in the temperature from θ = 1.0 to θ = 1.0 promotes a shift of the thermodynamic equilibrium position to lower χ_{B} and higher χ_{B} values. The simulated system behaves similarly to an exothermic reaction. As it is defined here, the proposed model considers the relationship between the activation energies of the forward and reverse reactions (_{A}E = 2, which implies _{ba}E = 2_{b}E), similar to an exothermic reaction._{a}^{19,25} This behaviour can be easily observed in Figure 6(a). Calculation of the reaction equilibrium constants, in terms of χ, for θ = 1.0 and θ = 1.0 yields K = 2.720 ± 0.006 and _{e}K = 1.949 ± 0.006, respectively. Thus, the results obtained from the simulation are consistent with several textbooks that discuss Le Chatelier's Principle._{e}^{24,25} The influence of the addition of the type-A reactant on the kinetics of the reaction is shown in Figure 6(b), where empty symbols represent the system without adding reactant and filled symbols represent the system with the addition of Δχ = 0.2 at _{A}t = 15. These results were obtained for θ = 1.0, and showed that the system quickly promoted the consumption of a portion of type-A molecules and the formation of type-B molecules soon after the addition of the A reactant (Δχ = 0.2). However, the system was shifted to a new thermodynamic equilibrium position with the new values of χ_{A} and χ_{A}. Similar to the results in Figure 6(a), these results are in agreement with Le Chatelier's principle, with the added reactant disturbing the thermodynamic equilibrium of the system._{B}Even though the equilibrium values of χ before and after adding the type-A molecules to the system are different, their ratios remain constant and equal to the thermodynamic equilibrium constant, _{B}K ≡χ_{e}/χ_{B}. The thermodynamic equilibrium constant was evaluated by measuring _{A}K at _{e}t ∈ [50,100] (not shown on the plot). The simulations without extra reactant added yielded K = 2.720 ± 0.006, whereas _{e}K = 2.718 ± 0.005 for the simulation with extra type-A reactant added. These results are in agreement with the expected theoretical value of the equilibrium constant, _{e}K = 2.718..., which is defined as the ratio between the kinetic constants of the forward and the reverse reactions, _{e}K = _{e}k/ _{f} k. Similarly, _{r}K can be defined here as the ratio between the forward and reverse transition probabilities, _{e}K = _{e}P/ _{f} P._{r}Finally, the results of the simulations for the catalytic and non-catalytic reversible reactions are compared in Figure 7. ^{0}_{A} = 0.8, χ^{0}_{B} = 0.0 and θ = 1.0 were considered in this case. The simulation of the catalytic reaction was accomplished by shifting the energy barrier by ΔE = -0.5. It was observed that the rate constant of the catalytic reaction was higher than the non-catalysed one. It was also observed that, regardless of the decrease in the energy barrier, the equilibrium position of the system converged to the same value for both cases. These results show an important aspect regarding the use of catalysts in chemical reactions, whereby these species accelerate the reaction rate (by increasing the rate constant) but do not affect the equilibrium position of the system._{cat}
Figure 7. Simulations for the catalytic and non-catalytic reversible reactions
The last examples we discuss are second-order reactions. Second-order reactions are categorised into two types, those that consist of two first-order reactants of the form A + B → C and those that consist of a single second-order reactant of the form 2A → C. and χ_{A}, our model simulates the kinetic behaviour of two second-order reactions, A + B → C, if χ_{B}^{0}_{A} ≠ χ^{0}_{B}, and 2A → C, if χ^{0}_{A} = χ^{0}_{B}.The time evolution of χ for the second-order reaction, 2A → C, is shown in Figure 8. The type-B molecule has the same behaviour as that observed for the type-A molecule since χ_{C}^{0}_{A} = χ^{0}_{B}= 0.5. Thus, it is expected that the reciprocal of the molar fraction exhibits a linear behaviour with time.^{24,25} The inset shows the dependence of as a function of time. The data provide a remarkable fit to , where = 2.006 ± 0.005 and k = 0.6059 ± 0.0001, which is in good agreement with the theoretical calculation (k = 0.606530...).^{24,25} This kinetic characteristic is valid for other values of q and initial molar fractions for χ and χ_{A}, with the constraint χ_{B}^{0}_{A} = χ^{0}_{B}.
Figure 8. The time evolution of χ_{A} and χ_{C} for the second-order reaction,2A → C
The behaviour of χ (B symbols) and χ_{B} (C symbols) for a second-order chemical reaction, A + B → C , with the initial condition χ_{C}_{A}^{0} ≠ χ_{B}^{0}, is shown in Figure 9(a). The predicted theoretical behaviour is that exhibits a linear dependence on time as follows: . We thus carried out simulations considering χ_{A}^{0} = 0.5, χ_{B}^{0} = 0.2, χ_{C}^{0} = 0.0 and θ = 1.0 to investigate this situation. The simulated results are plotted in Figure 9(b), where is plotted on a log scale as a function of time. The simulated data provided an exceptional fit to , which was approximately equal to the theoretical result of .^{19,24,25} These initial parameters, as well as the other tested parameters, yielded results that reproduced the theoretical results.
Figure 9. (a) The time evolution of χ_{A}, χ_{B} and χ_{C} for the second-order reaction,A + B → C, (b) -χ on a log scale as a function of time
Here we discuss how the number of the molecules ( Figure 10 shows χ ^{2} (filled circles) and 10^{5} (empty circles) molecules. The simulated results were obtained using χ^{0} = χ_{A}^{0} = 0.5, χ_{B}^{0} = 0.0, θ = 1.0 and _{C}N = 1. The fluctuations on the temporal evolution of χ_{s} diminishes as _{A}N increases, and the simulated data merge toward the expected theoretical result (straight line). As seen in Table 1, the simulated rate constant (k) possesses an increasingly better agreement with the theoretical one (_{s}k) as _{t}N increases.
Figure 10. χ_{A} versus time for N = 10^{2} (filled circles) and 10^{5} (empty circles) molecules
The effect of multiple simulations is presented in Table 2. The data were obtained with N = 10 ^{0} = 0.5, χ_{B}^{0} = 0.0, θ = 1.0. As can be seen, _{C}k approaches _{s}k as _{t}N increases, evidencing a similar effect when compared to the increase in _{s}N.
A simple mathematical model has been used to investigate the equilibrium and kinetic properties of elementary homogeneous chemical reactions. The reaction dynamics between the reactants were proposed via the Monte Carlo method and the Ehrenfest urn model. The Ehrenfest urn model was used to define the order of the reaction, while the Monte Carlo method, with the Metropolis update, was used to access the main physical chemistry features as a function of activation energy, steric factor and temperature. All the basic characteristics of the reactions, such as the reaction order, the Arrhenius' Law and Le Chatelier's principle, were recovered with this simple model. Although the methodology presented here was only applied to first- and second-order reactions, we argue that it can be also used to study more complex (higher order) reactions. Here we demonstrated the effectiveness of our model as a useful tool for chemistry-related class activities at all education levels, as well in scientific research that covers physical chemistry topics related to thermodynamics and reaction kinetics, which highlights the wide applicability of such a model. It is worth pointing out the advantages of this present model. One of them is the low cost to carry out such a study when compared with the experimental laboratory procedures. Furthermore, the lack of laboratory facilities limits the activities and research in physical chemistry at some schools and universities. Additionally, from a mathematical point of view, the model seems to be useful as a complementary tool to the traditional methodologies that employ differential ordinary equations and integrated rate laws amplifying mathematical concepts necessary to a better understanding of the phenomena in chemical kinetics. It is important to highlight that the Monte Carlo step, defined and used in this work, does not have at first a direct relationship with the real time. However, it is possible to propose a relationship between them, in order to use the simulations to treat real problems in chemical kinetics. One possibility would be to establish a connection between the MCS and some characteristic time of the chemical kinetics in a real problem, as proposed by López-Castilho and Souza Filho in which the step was related to the half-life time. Finally, we highlight that it is possible to develop a computational software interface with the algorithm reported here to facilitate its manipulation by users, be they teachers or students. This software development, as well as other results, are currently being undertaken.
This work was partially supported by UESB, CAPES and CNPq (Brazilian Agency). Also, N. M. Oliveira-Neto thanks the ICTP (International Center for Theoretical Physics) in Trieste,Italy, where this work was partially developed.
1. Barrow, L. H.; 2. Zômpero, A. F.; Laburú, C. E.; 3. Finkelstein, N. D.; Adams, W. K.; Keller, C. J.; Kohl, P. B.; Perkins, K. K.; Podolefsky, N. S.; Reid, S.; 4. Johnson, K. J.; 5. Dixon, D. A.; Shafer, R.; 6. Para, A. F.; Lazzarini, E.; 7. Moebs, W. D.; Haglund, E. A.; 8. Tellinghuisen, J.; 9. Bentenitis, N.; 10. López-Castilho, A.; de Souza Filho, J. C.; 11. Farias, R. R.; Cardoso, L. A. M.; Oliveira-Neto, N. M.; 12. Farias, R. R.; Cardoso, L. A. M.; Oliveira-Neto, N. M.; Junior, B. B. N.; 13. Novo, J. B. M.; Júnior, L. C. D.; 14. Hollauer, E.; 15. Winnischofer, H.; Araújo, M. P.; Júnior, L. C. D.; Novo, J. B. M.; 16. Silva-Dias, L.; López-Castilho, A.; 17. Júnior, L. C. D.; Novo, J. B. M.; 18. Binder, K.; Heermann, D. W.; 19. Houston, P. L.; 20. Oliveira-Neto, N. M.; Curado, E. M. F.; Nobre, F. D.; Rego-Monteiro, M. A.; 21. Oliveira-Neto, N. M.; Curado, E. M. F.; Nobre, F. D.; Rego-Monteiro, M. A.; 22. Salinas, S. R. A.; 23. Karlin, S.; McGregor, J.; 24. Atkins, P. W.; 25. Levine, I. N.; |

On-line version ISSN 1678-7064 Printed version ISSN 0100-4042

Química Nova

Publicações da Sociedade Brasileira de Química

Caixa Postal: 26037
05513-970 São Paulo - SP

Tel/Fax: +55.11.3032.2299/+55.11.3814.3602

Free access