JBCS

INTRODUCTION

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.^1,2 Inquiry, from a teaching perspective, can be viewed as teaching where the students are actively, instead of passively, engaged in their learning process, with the teacher serving as a mediator. In this active learning environment, the activities and discussion are focused on the students, with more emphasis given on process skills, investigative activities of science questions, public communications and working as a member of a team. The practical activities thus represent a natural and important teaching strategy, particularly the use of computer simulations.³ Since the advent of computers, numerous efforts have been made to employ computers in the classroom. Teaching modules in chemical kinetics have used stochastic methods, such as the Monte Carlo method (MCM), to study chemical reactions.^4-17

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.^12,13

A simple stochastic approach based on MCM and the Ehrenfest urn model (EUM)¹² is proposed here to overcome some of the difficulties associated with teaching chemistry by simulating the chemical kinetics of elementary reactions. This hybrid model takes into account the main kinetics parameters, such as temperature and activation energy, thus making it is possible to obtain features in chemical kinetics.

EXPERIMENTAL

Mathematical model

The standard MCM,¹⁸ defined via the Metropolis update, is applied. Here the transition probability functions depend on the energy barrier (activation energy), temperature and steric factor.^11,12,19 Thus, when considering the elementary reaction between two species A B, the transition probabilities, P_f and P_r, for the forward and reverse reactions are defined as follows

where β = (k_BT)^-1, k_B is the Boltzmann constant, T is the temperature, p_f^s and p_r^s are the steric factors for the forward and the reverse reaction, respectively, and E_a and E_b are the activation energies for the forward and the reverse steps, respectively. For an irreversible reaction, Equation (2) is not taken into account, as P_r = 0, which is equivalent to assuming E_B → ∞. This part of the present approach, defined by Equations (1) and (2), is quite similar to that used in the hydrogen atom ionisation problem.^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, θ ≡ k_BT/E_a, and dimensionless energy-related variable, E_ba ≡ E_b /E_a, are defined to simplify the calculations. The computational simulations carried out in this work are thus based on the thermal energy k_BT and E_b, both of which are measured with respect to 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.^10,22,23 Consider a model comprised of two urns, one for type-A molecules and the other for type-B molecules, which contain M boxes labelled i = 1,2,...,M. The urns contain N_A and N_B molecules corresponding to the type-A and type-B molecules, respectively. If a given box A_i is occupied with a type-A molecule, the corresponding B_i box will be empty and vice versa. In addition to the definition given in the literature,^10-12 the relationship between the total number of molecules is defined here as N ≡ N_A + N_B, and the total number of available boxes 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 i ∈ [1,M] and check if there is such a molecule in the box A_i. If the choice fails, continue to step (II). If the choice succeeds, a uniform random number z ∈ [0,1] is generated and compared to P_f from Equation (1). If z ≤ P_f, the molecule indexed as i is eliminated from box A_i, and a type-B molecule is created in box B_i, which indicates that the reaction occurred, transforming A into B. If z > P_f, the type-A molecule remains in its original form.

(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 ∈ [0,1] is generated and compared to P_r from Equation (2). If z ≤ P_r, the reaction occurs and B is transformed into A. If z > P_r, the type-B molecule remains in its original form.

By running steps (I) and (II) M times, one Monte Carlo step (MCS) is defined. If the reaction is irreversible, step (II) is not performed.

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, E ≡ E_a, and then transform A and B into product C. The following step is performed in this case.

(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 ∈ [0,1] is generated and compared to P_f from Equation (1). If z ≤ P_f the reaction occurs and A reacts with B to form C. If z > P_f, the A and B molecules remain unreacted. As previously mentioned above, by performing step (III) 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, χ_A and χ_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⁶ 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.

RESULTS AND DISCUSSION

Irreversible reaction

The irreversible first-order reaction A → B, will be studied first. Figure 1 shows χ_A and χ_B as a function of time, in MCSs. The simulation started with all type-A molecules, where χ⁰_A= 1.0, χ⁰_B = 0,0, θ = 1.0 and steric factor p_f^s = p^s = 1.0. As can be seen, χ_A decreases until these species is completely consumed, whereas χ_B approaches unity. The inset figure shows the behaviour of χ_A on a log scale, showing the exponential time dependence. This result fits well to

Figure 1. χ_A and χ_B as a function of time, in MCSs

where k is the rate constant and t is the time in MCSs.

The effect of temperature on the rate constant is shown in Figure 2, where χ_A is plotted as a function of time for different temperatures. Since θ = k_BT/E_a, this behaviour can be interpreted in two different ways.^11,12 One is that the increase in θ means that the k_BT 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 k_BT and decreasing E_a. Therefore, we also investigated different systems (different E_a) with the same temperature. Considering the set of parameters (θ values) used in Figure 2, the same relationship between the rate constant (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 A = 0.99844 and b = 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 half-life (τ_1/2). Using this model, the simulation yielded^11,12

where c = 0.69122 ± 0.0009, which is in good agreement with the theoretical first-order behaviour that yielded an expected value of ln(2) = 0.69314... . These results thus align with the theory of chemical kinetics that is widely discussed in many elementary textbooks.^19,24,25

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 χ_A is plotted as a function of p^s for θ = 1.0. Decreasing p^s decreases 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^s = 1.0 means that all collisions are geometrically favourable.

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 p^s = 0.5 was thus performed to check this phenomenon.

The rate constant k, on a log scale, is a function of the reciprocal of θ as shown in Figure 4(a). The solid line is the fitted function k = Ae^-b'/θ, 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^s values, it is seen that A ~ p^s, thus indicating that this approach is useful in dealing with the kinetics of irreversible chemical reactions. As mentioned above, it is expected that A = p^s and 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 is shown as a function of 1/k in Figure 4(b). Even though k is different for other p^s values, the dependence of τ_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 + ΔE. If ΔE > 0, there is a "negative" catalysis or inhibition. Otherwise, if DE < 0, there is a "positive catalysis".²⁴ The kinetics are thus defined by the parameter ΔE_cat ≡ ΔE/E_a. 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 (ΔE_cat < 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 ΔE_cat > 0.

Figure 5. The kinetic profiles of the catalytic and non-catalytic reactions

Reversible reaction

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_f^s = p_r^s = 1. Unlike the irreversible reactions, it is necessary to define the dimensionless energy parameter, E_ba ≡ E_b/E_a, which represents the reverse reaction energy barrier (rationed) to E_a. Some results are shown in Figure 6, where E_ba = 2.0 are considered, with χ⁰_A = 0.8 and χ⁰_B = 0.0. The behaviours of χ_A and χ_B are shown in Figure 6(a) for two distinct temperatures, θ = 1.0 (empty squares and circles) and θ = 1.5 (filled squares and circles).

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, χ_A and χ_B 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 χ_A 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 (E_ba = 2, which implies E_b = 2E_a), similar to an exothermic reaction.^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_e = 2.720 ± 0.006 and K_e = 1.949 ± 0.006, respectively. Thus, the results obtained from the simulation are consistent with several textbooks that discuss Le Chatelier's Principle.^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 Δχ_A = 0.2 at 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 (Δχ_A = 0.2). However, the system was shifted to a new thermodynamic equilibrium position with the new values of χ_A and χ_B. 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.

Even though the equilibrium values of χ_A andχ_B before and after adding the type-A molecules to the system are different, their ratios remain constant and equal to the thermodynamic equilibrium constant, K_e ≡χ_B/χ_A. The thermodynamic equilibrium constant was evaluated by measuring K_e at t ∈ [50,100] (not shown on the plot). The simulations without extra reactant added yielded K_e = 2.720 ± 0.006, whereas K_e = 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, K_e = 2.718..., which is defined as the ratio between the kinetic constants of the forward and the reverse reactions, K_e = k_f / k_r. Similarly, K_e can be defined here as the ratio between the forward and reverse transition probabilities, 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. E_ba = 2.0, χ⁰_A = 0.8, χ⁰_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_cat = -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.

Figure 7. Simulations for the catalytic and non-catalytic reversible reactions

Second-order 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.¹⁹ These irreversible reactions only require Equation (1) in the mathematical model, with p^s = 1, following step (III) in the experimental section. Depending on the initial conditions of χ_A and χ_B, our model simulates the kinetic behaviour of two second-order reactions, A + B → C, if χ⁰_A ≠ χ⁰_B, and 2A → C, if χ⁰_A = χ⁰_B.

The time evolution of χ_A and χ_C 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 χ⁰_A = χ⁰_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 χ_A and χ_B, with the constraint χ⁰_A = χ⁰_B.

Figure 8. The time evolution of χ_A and χ_C for the second-order reaction,2A → C

The behaviour of χ_A (A symbols), χ_B (B symbols) and χ_C (C symbols) for a second-order chemical reaction, A + B → C , with the initial condition χ_A⁰ ≠ χ_B⁰, 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.5, χ_B⁰ = 0.2, χ_C⁰ = 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

Effect of the number of molecules (N) and the sampling (N_s) on rate constant for a second- order reaction (2A → C)

Here we discuss how the number of the molecules (N) and the sampling (N_s) influence the results for a given rate constant and a second-order reaction of the type 2A → C.

Figure 10 shows χ_A versus time for N = 10² (filled circles) and 10⁵ (empty circles) molecules. The simulated results were obtained using χ⁰_A = χ⁰_B = 0.5, χ⁰_C = 0.0, θ = 1.0 and N_s = 1. The fluctuations on the temporal evolution of χ_A diminishes as N increases, and the simulated data merge toward the expected theoretical result (straight line). As seen in Table 1, the simulated rate constant (k_s) possesses an increasingly better agreement with the theoretical one (k_t) as N increases.

Figure 10. χ_A versus time for N = 10² (filled circles) and 10⁵ (empty circles) molecules

The effect of multiple simulations is presented in Table 2. The data were obtained with N = 10³ molecules and χ⁰_A = χ⁰_B = 0.5, χ⁰_C = 0.0, θ = 1.0. As can be seen, k_s approaches k_t as N_s increases, evidencing a similar effect when compared to the increase in N.

CONCLUSIONS AND PROSPECTS

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.¹⁰ Another possibility would be to propose a proportionality between the Monte Carlo step and the time between successive molecular collisions, given by the collision theory, as discussed by Oliveira-Neto et al., in the study of hydrogen atom ionization problem.²⁰

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.

ACKNOWLEDGEMENTS

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.

REFERENCES

1. Barrow, L. H.; J. Sci. Teacher Educ. 2006, 17, 265.

2. Zômpero, A. F.; Laburú, C. E.; Revista Ensaio 2011, 13, 67.

3. Finkelstein, N. D.; Adams, W. K.; Keller, C. J.; Kohl, P. B.; Perkins, K. K.; Podolefsky, N. S.; Reid, S.; Physical Review Physics Education Research 2005, 1, 1.

4. Johnson, K. J.; J. Chem. Educ. 1970, 47, 819.

5. Dixon, D. A.; Shafer, R.; J. Chem. Educ. 1973, 50, 648.

6. Para, A. F.; Lazzarini, E.; J. Chem. Educ. 1974, 336.

7. Moebs, W. D.; Haglund, E. A.; J. Chem. Educ. 1976, 53, 506.

8. Tellinghuisen, J.; J. Chem. Educ. 2005, 82, 1709.

9. Bentenitis, N.; J. Chem. Educ. 2008, 85, 1146.

10. López-Castilho, A.; de Souza Filho, J. C.; Quim. Nova 2007, 30, 1759.

11. Farias, R. R.; Cardoso, L. A. M.; Oliveira-Neto, N. M.; Revista Científica do Departamento de Química e Exatas 2011, 2, 7.

12. Farias, R. R.; Cardoso, L. A. M.; Oliveira-Neto, N. M.; Junior, B. B. N.; Quim. Nova 2013, 5, 729.

13. Novo, J. B. M.; Júnior, L. C. D.; Quim. Nova 2011, 34, 707.

14. Hollauer, E.; Quim. Nova 1988, 11, 186.

15. Winnischofer, H.; Araújo, M. P.; Júnior, L. C. D.; Novo, J. B. M.; Quim. Nova 2010, 33, 225.

16. Silva-Dias, L.; López-Castilho, A.; Quim. Nova 2015, 38, 1232.

17. Júnior, L. C. D.; Novo, J. B. M.; Quim. Nova 2014, 37, 361.

18. Binder, K.; Heermann, D. W.; Monte Carlo Simulations in Statistical Physics, 2^nd ed., Springer-Verlag: New York, 1992.

19. Houston, P. L.; Chemical Kinetics and Reaction Dynamics, 2^nd ed., Dover Publications, Inc.: New York, 2006.

20. Oliveira-Neto, N. M.; Curado, E. M. F.; Nobre, F. D.; Rego-Monteiro, M. A.; Physica A 2007, 374, 251.

21. Oliveira-Neto, N. M.; Curado, E. M. F.; Nobre, F. D.; Rego-Monteiro, M. A.; J. Phys. B: At., Mol. Opt. Phys. 2007, 40, 1975.

22. Salinas, S. R. A.; Introduçao à Física Estatística, 1^st ed.; EDUSP.: Sao Paulo, 1999.

23. Karlin, S.; McGregor, J.; Journal of Applied Probability 1965, 2, 352.

24. Atkins, P. W.; Físico-Química, 6^th ed.; LTC.: Rio de Janeiro, 1999.

25. Levine, I. N.; Físico-Química, 6^th ed.; LTC-Gen.: Rio de Janeiro, 2012.