Langmuir-Hinshelwood model: Acceleration by Automated Rescaling of the Rate Constants.

Note

To follow this tutorial, either:

The scale disparity issue is common in Kinetic Monte Carlo simulations. It emerges as a result of the fact that some fundamental events or groups of them typically occur at vastly different time scales; in other words, their rate constants can span multiple orders of magnitude. In heterogeneous catalysis, there are typically two groups: 1) very fast events that correspond to the species’ surface diffusions and 2) slow reactions that change their chemical identity. The latter group of events is usually the one of interest because it allows the evaluation of the material’s catalytic activity. In contrast, the species’ surface diffusion does not contribute significantly to the net evolution of the slow reactions. But, as becomes the most frequent step, it also becomes the limiting factor of the simulation progress, considerably increasing the computational cost. This tutorial shows how to speed up the calculation by several orders of magnitude without sacrificing precision by automatically detecting and scaling the rate constants of fast reactions.

We will focus on the net reaction \(\text{CO}+\frac{1}{2}\text{O}_2\longrightarrow \text{CO}_2\) that takes place at a catalyst’s surface and whose reaction mechanism is described by the Langmuir-Hinshelwood model. Because this model has four very fast processes (\(CO\) and \(O_2\) adsorption, and \(O*\) and \(CO*\) diffusion) and one slow process (\(CO\) oxidation), it is an ideal prototype for demonstrating the benefits of the automated rescaling of rate constants technique. Our ultimate goal is to investigate how altering the relative percentage of the gas reactants \(CO\) and \(O_2\) (at a specific temperature and pressure) affect the rate of \(CO_2\) production under steady state conditions. This example is inspired on Zacros tutorial What’s KMC All About and Why Bother.

Let’s start! The first step is to import all packages we need:

import multiprocessing
import numpy
import scm.plams
import scm.pyzacros as pz
import scm.pyzacros.models

Then, we initialize the pyZacros environment.

scm.pyzacros.init()

PLAMS working folder: /home/user/pyzacros/examples/ZiffGulariBarshad/plams_workdir

Notice this command created the directory where all Zacros input and output files will be stored if they are needed for future reference (plams_workdir by default). Typically, the user doesn’t need to use these files.

This calculation necessitates a significant computational effort. On a typical laptop, it should take around 20 min to complete. So, in order to speed things up, we’ll use the plams.JobRunner class to run as many parallel instances as possible. In this case, we choose to use the maximum number of simultaneous processes (maxjobs) equal to the number of processors in the machine.

maxjobs = multiprocessing.cpu_count()
scm.plams.config.default_jobrunner = scm.plams.JobRunner(parallel=True, maxjobs=maxjobs)
print('Running up to {} jobs in parallel simultaneously'.format(maxjobs))

Running up to 8 jobs in parallel simultaneously

First, we initialize our Langmuir-Hinshelwood model, which by luck is available as a predefined model in pyZacros,

lh = pz.models.LangmuirHinshelwood()

Then, we must set up a ZacrosParametersScanJob calculation, which will allow us to scan the molar fraction of \(CO\) as a parameter. However, this calculation requires the definition of a ZacrosSteadyStateJob, that in turns requires a ZacrosJob. So, We will go through them one at a time:

1. Setting up the ZacrosJob

For ZacrosJob, all parameters are set using a Setting object. To begin, we define the physical parameters: temperature (in K), and pressure (in bar). The calculation parameters are then set: species numbers (in s) determines how frequently information about the number of gas and surface species will be stored, max time (in s) specifies the maximum allowed simulated time, and “random seed” specifies the random seed to make the calculation precisely reproducible. Keep in mind that max time defines the calculation’s stopping criterion, and it is the parameter that will be controlled later to achieve the steady-state configuration. Finally, we create the ZacrosJob, which uses the parameters we just defined as well as the Langmuir-Hinshelwood model’s lattice, mechanism, and cluster expansion. Notice we do not run this job, we use it as a reference for the steady-state calculation described below.

z_sett = pz.Settings()
z_sett.temperature = 500.0
z_sett.pressure = 1.000
z_sett.species_numbers = ('time', 1.0e-5)
z_sett.max_time = 100*1.0e-5
z_sett.random_seed = 1609

z_job = pz.ZacrosJob( settings=z_sett, lattice=lh.lattice,
                      mechanism=lh.mechanism,
                      cluster_expansion=lh.cluster_expansion )

2. Setting up the ZacrosSteadyStateJob

We also need to create a Setting object for ZacrosJob There, we ask for a steady-state configuration using a TOFs calculation with a 95% confidence level (turnover frequency.confidence), using four replicas to speed up the calculation (turnover frequency.nreplicas); for more information, see example Ziff-Gulari-Barshad model: Steady State Conditions. Then, we ask for the rate constants to be automatically scaled (scaling.enabled) using an inspection time of 0.0006 s (“scaling.max time”). In a nutshell, the scaling algorithm uses this maximum time and the original rate constants to execute a probe simulation. From there, the occurrence rates for each reaction are determined and scaled appropriately. Only reactions that are proven to be quasi-equilibrated are modified. The actual simulation is then started from the beginning using the new reaction rates following this. The ZacrosSteadyStateJob.Parameters object allows to set the grid of maximum times to explore in order to reach the steady state (ss_params). Finally, we create ZacrosSteadyStateJob, which references the ZacrosJob defined above (z_job) as well as the Settings object and parameters we just defined:

ss_sett = pz.Settings()
ss_sett.turnover_frequency.confidence = 0.95
ss_sett.turnover_frequency.nreplicas = 4
ss_sett.scaling.enabled = 'T'
ss_sett.scaling.max_time = 60*1e-5

ss_params = pz.ZacrosSteadyStateJob.Parameters()
ss_params.add( 'max_time', 'restart.max_time',
                2*z_sett.max_time*( numpy.arange(10)+1 )**2 )

ss_job = pz.ZacrosSteadyStateJob( settings=ss_sett, reference=z_job,
                                  parameters=ss_params )

[27.01|09:15:05] JOB plamsjob Steady State Convergence: Using nbatch=20,confidence=0.95,ignore_nbatch=1,nreplicas=4

3. Setting up the ZacrosParametersScanJob

Although the ZacrosParametersScanJob does not require a Setting object, it does require a ZacrosSteadyStateJob.Parameters object to specify which parameters must be modified systematically. In this instance, all we need is a dependent parameter, the \(O_2\) molar fraction x_O2, and an independent parameter, the \(CO\) molar fraction x_CO, which ranges from 0.05 to 0.95. Keep in mind that the condition x_CO+x_O2=1 must be met. These molar fractions will be used internally to replace molar fraction.CO and molar fraction.O2 in the Zacros input files. Then, using the ZacrosSteadyStateJob defined earlier (ss job) and the parameters we just defined (ps params), we create the ZacrosParametersScanJob:

ps_params = pz.ZacrosParametersScanJob.Parameters()
ps_params.add( 'x_CO', 'molar_fraction.CO', numpy.linspace(0.05, 0.95, 11) )
ps_params.add( 'x_O2', 'molar_fraction.O2', lambda params: 1.0-params['x_CO'] )

ps_job = pz.ZacrosParametersScanJob( reference=ss_job, parameters=ps_params )

[27.01|09:15:05] JOB ps_cond000 Steady State Convergence: Using nbatch=20,confidence=0.95,ignore_nbatch=1,nreplicas=4
[27.01|09:15:05] JOB ps_cond001 Steady State Convergence: Using nbatch=20,confidence=0.95,ignore_nbatch=1,nreplicas=4
[27.01|09:15:05] JOB ps_cond002 Steady State Convergence: Using nbatch=20,confidence=0.95,ignore_nbatch=1,nreplicas=4
[27.01|09:15:05] JOB ps_cond003 Steady State Convergence: Using nbatch=20,confidence=0.95,ignore_nbatch=1,nreplicas=4
[27.01|09:15:05] JOB ps_cond004 Steady State Convergence: Using nbatch=20,confidence=0.95,ignore_nbatch=1,nreplicas=4
[27.01|09:15:05] JOB ps_cond005 Steady State Convergence: Using nbatch=20,confidence=0.95,ignore_nbatch=1,nreplicas=4
[27.01|09:15:05] JOB ps_cond006 Steady State Convergence: Using nbatch=20,confidence=0.95,ignore_nbatch=1,nreplicas=4
[27.01|09:15:05] JOB ps_cond007 Steady State Convergence: Using nbatch=20,confidence=0.95,ignore_nbatch=1,nreplicas=4
[27.01|09:15:05] JOB ps_cond008 Steady State Convergence: Using nbatch=20,confidence=0.95,ignore_nbatch=1,nreplicas=4
[27.01|09:15:05] JOB ps_cond009 Steady State Convergence: Using nbatch=20,confidence=0.95,ignore_nbatch=1,nreplicas=4
[27.01|09:15:05] JOB ps_cond010 Steady State Convergence: Using nbatch=20,confidence=0.95,ignore_nbatch=1,nreplicas=4

The parameters scan calculation setup is ready. Therefore, we can start it by invoking the function run(), which will provide access to the results via the results variable after it has been completed. The sentence involving the method ok(), verifies that the calculation was successfully executed, and waits for the completion of every executed thread. Go and grab a cup of coffee, this step will take around 15 mins!

results = ps_job.run()

if not ps_job.ok():
    print('Something went wrong!')

[27.01|09:15:05] JOB plamsjob STARTED
[27.01|09:15:05] Waiting for job plamsjob to finish
[27.01|09:15:05] JOB plamsjob RUNNING
[27.01|09:15:05] JOB plamsjob/ps_cond000 STARTED
[27.01|09:15:05] JOB plamsjob/ps_cond001 STARTED
[27.01|09:15:05] JOB plamsjob/ps_cond002 STARTED

...

[27.01|09:31:03]       species            TOF          error          ratio     conv?
[27.01|09:31:03]            CO      -52.52688        4.84537        0.09225     False
[27.01|09:31:03]            O2      -26.90450        2.78063        0.10335     False
[27.01|09:31:03]           CO2       52.91640        4.65244        0.08792     False
[27.01|09:31:03]    Replica #2
[27.01|09:31:03]       species            TOF          error          ratio     conv?
[27.01|09:31:04]            CO     -225.40015       12.98476        0.05761     False
[27.01|09:31:04]            O2     -115.66183        4.44209        0.03841      True
[27.01|09:31:04]           CO2      228.95111        8.73665        0.03816      True
[27.01|09:31:04]    Average
[27.01|09:31:04]       species            TOF          error          ratio     conv?
[27.01|09:31:04]            CO      -54.87733        1.63967        0.02988      True
[27.01|09:31:04]            O2      -26.85096        1.08631        0.04046      True
[27.01|09:31:04]           CO2       54.73375        1.56100        0.02852      True
[27.01|09:31:04] JOB plamsjob/ps_cond000 Steady State Convergence: CONVERGENCE REACHED. DONE!
[27.01|09:31:04] JOB plamsjob/ps_cond000 FINISHED
[27.01|09:31:05]    Replica #3
[27.01|09:31:05]       species            TOF          error          ratio     conv?
[27.01|09:31:05]            CO     -216.58663       13.48940        0.06228     False
[27.01|09:31:05]            O2     -106.99160        5.19478        0.04855      True
[27.01|09:31:05]           CO2      214.75940       10.88629        0.05069     False
[27.01|09:31:06]    Average
[27.01|09:31:06]       species            TOF          error          ratio     conv?
[27.01|09:31:06]            CO     -221.49798        5.07706        0.02292      True
[27.01|09:31:06]            O2     -110.17810        2.42785        0.02204      True
[27.01|09:31:06]           CO2      220.65690        4.39982        0.01994      True
[27.01|09:31:06] JOB plamsjob/ps_cond010 Steady State Convergence: CONVERGENCE REACHED. DONE!
[27.01|09:31:06] JOB plamsjob/ps_cond010 FINISHED
[27.01|09:31:06] JOB plamsjob/ps_cond000 SUCCESSFUL
[27.01|09:31:07] JOB plamsjob/ps_cond010 SUCCESSFUL
[27.01|09:31:11] JOB plamsjob FINISHED
[27.01|09:31:15] JOB plamsjob SUCCESSFUL

If the execution got up to this point, everything worked as expected. Hooray!

Finally, in the following lines, we just nicely print the results in a table. See the API documentation to learn more about how the results object is structured, and the available methods. In this case, we use the turnover_frequency() and average_coverage() methods to get the TOF for the gas species and average coverage for the surface species, respectively. Regarding the latter one, we use the last ten steps in the simulation to calculate the average coverages.

x_CO = []
ac_O = []
ac_CO = []
TOF_CO2 = []

results_dict = results.turnover_frequency()
results_dict = results.average_coverage( last=10, update=results_dict )

for i in range(len(results_dict)):
    x_CO.append( results_dict[i]['x_CO'] )
    ac_O.append( results_dict[i]['average_coverage']['O*'] )
    ac_CO.append( results_dict[i]['average_coverage']['CO*'] )
    TOF_CO2.append( results_dict[i]['turnover_frequency']['CO2'] )

print( '------------------------------------------------' )
print( '%4s'%'cond', '%8s'%'x_CO', '%10s'%'ac_O', '%10s'%'ac_CO', '%12s'%'TOF_CO2' )
print( '------------------------------------------------' )
for i in range(len(x_CO)):
    print( '%4d'%i, '%8.2f'%x_CO[i], '%10.6f'%ac_O[i], '%10.6f'%ac_CO[i], '%12.6f'%TOF_CO2[i] )

------------------------------------------------
cond     x_CO       ac_O      ac_CO      TOF_CO2
------------------------------------------------
   0     0.05   0.665594   0.028219    54.733748
   1     0.14   0.615375   0.082844   135.299544
   2     0.23   0.582250   0.126594   198.215992
   3     0.32   0.532375   0.178719   257.936395
   4     0.41   0.497031   0.221625   300.382949
   5     0.50   0.442750   0.272969   329.214754
   6     0.59   0.402031   0.318562   352.259276
   7     0.68   0.351687   0.372437   351.805662
   8     0.77   0.297906   0.422687   332.247334
   9     0.86   0.232156   0.488156   310.444388
  10     0.95   0.144281   0.554344   220.656898

The results table above demonstrates that when \(CO\) coverage rises, the net \(CO\) oxidation reaction tends to progress more quickly until it reaches a maximum of about 0.7. At this point, it begins to decline. Notice that the maximal \(CO_2\) generation occurs simultaneously as the \(CO*\) coverage reaches parity with the \(O*\) coverage. This makes intuitive sense since, at that point, we maximize the likelihood of discovering an \(O*\) and a \(CO*\) that are close enough to react in accordance with the proper reaction’s stoichiometry. The number of \(O*\) or \(CO*\) observers rises as we move away from that critical point, decreasing the likelihood that the \(CO\) oxidation process could take place.

Now, we can close the pyZacros environment:

scm.pyzacros.finish()

[27.01|09:31:31] PLAMS run finished. Goodbye

Comparing with Traditional Kinetic Models

Note

To follow the second part of this tutorial, either:

The goal of this second part of the tutorial is to demonstrate how to resume the calculation from the first part and visualize the results. In addition, we will compare these results to an analytical model to assess their quality. Before we begin, make sure we have the working directory generated by CoverageAndReactionRate.ipynb or CoverageAndReactionRate.py. We’ll assume it’s called plams_workdir, which is the default value. However, keep in mind that each time you run any of the above scripts, PLAMS will create a new working directory by appending a sequential number to its name, for example, plams_workdir.001. If this is the case, simply replace plams_workdir with the appropriate value.

So, let’s get started!

First, we load the required packages and retrieve the ZacrosParametersScanJob (job) and corresponding results object (results) from the working directory, as shown below:

import scm.pyzacros as pz
import scm.pyzacros.models

scm.pyzacros.init()

job = scm.pyzacros.load( 'plams_workdir/plamsjob/plamsjob.dill' )
results = job.results

scm.pyzacros.finish()

PLAMS working folder: /home/user/pyzacros/examples/ZiffGulariBarshad/plams_workdir.002
[27.01|09:45:00] PLAMS run finished. Goodbye

To be certain, we generate and print the same summary table from the end of the first part of the tutorial. They must be exactly the same:

x_CO = []
ac_O = []
ac_CO = []
TOF_CO2 = []

results_dict = results.turnover_frequency()
results_dict = results.average_coverage( last=10, update=results_dict )

for i in range(len(results_dict)):
    x_CO.append( results_dict[i]['x_CO'] )
    ac_O.append( results_dict[i]['average_coverage']['O*'] )
    ac_CO.append( results_dict[i]['average_coverage']['CO*'] )
    TOF_CO2.append( results_dict[i]['turnover_frequency']['CO2'] )

print( '------------------------------------------------' )
print( '%4s'%'cond', '%8s'%'x_CO', '%10s'%'ac_O', '%10s'%'ac_CO', '%12s'%'TOF_CO2' )
print( '------------------------------------------------' )
for i in range(len(x_CO)):
    print( '%4d'%i, '%8.2f'%x_CO[i], '%10.6f'%ac_O[i], '%10.6f'%ac_CO[i], '%12.6f'%TOF_CO2[i] )

------------------------------------------------
cond     x_CO       ac_O      ac_CO      TOF_CO2
------------------------------------------------
   0     0.05   0.665594   0.028219    54.733748
   1     0.14   0.615375   0.082844   135.299544
   2     0.23   0.582250   0.126594   198.215992
   3     0.32   0.532375   0.178719   257.936395
   4     0.41   0.497031   0.221625   300.382949
   5     0.50   0.442750   0.272969   329.214754
   6     0.59   0.402031   0.318562   352.259276
   7     0.68   0.351687   0.372437   351.805662
   8     0.77   0.297906   0.422687   332.247334
   9     0.86   0.232156   0.488156   310.444388
  10     0.95   0.144281   0.554344   220.656898

Additionally, you can see the aforementioned results visually if you have installed the package matplotlib. Please review the code below.

import matplotlib.pyplot as plt

fig = plt.figure()

ax = plt.axes()
ax.set_xlim([0.0, 1.0])
ax.set_xlabel('Molar Fraction CO', fontsize=14)
ax.set_ylabel('Coverage Fraction (%)', color='blue', fontsize=14)
ax.plot(x_CO, ac_O, marker='$\u25CF$', color='blue', linestyle='-.', markersize=4, zorder=2)
ax.plot(x_CO, ac_CO, marker='$\u25EF$', color='blue', markersize=4, zorder=4)
plt.text(0.3, 0.60, 'O*', fontsize=18, color='blue')
plt.text(0.7, 0.45, 'CO*', fontsize=18, color='blue')

ax2 = ax.twinx()
ax2.set_ylabel('TOF (mol/s/site)',color='red', fontsize=14)
ax2.plot(x_CO, TOF_CO2, marker='$\u25EF$', color='red', markersize=4, zorder=6)
plt.text(0.3, 200.0, 'CO$_2$', fontsize=18, color='red')

plt.show()
/scm-uploads/doc/pyzacros/_images/CoveragesAndReactionRate_ViewResults_5_0.png

Starting from the left side of the figure, it shows that by increasing the x_CO; the net \(CO\) oxidation reaction tends to progress more quickly (like a first-order reaction). Then the reaction rate peaks at around x_CO=0.7, becoming almost independent of x_CO (like a zero-order reaction). And finally, it then begins to decline (like a negative-order reaction). Furthermore, we can see that as we increase the coverage of \(CO*\), we decrease the coverage of \(O*\). As a result, when we mostly cover the surface with one of the reactants (either \(CO*\) or \(O*\)), the rate of \(CO_2\) production becomes slow because there isn’t enough of the complementary reactant on the surface. So, it is not a surprise, that the maximum \(CO_2\) generation coincides with the point at which \(CO*\) coverage equals \(O*\) coverage.

In the Zacros tutorial What’s KMC All About and Why Bother?, there is a thorough discussion of Langmuir-Hinshelwood-type models, their approximations with respect to a real system, and how to obtain analytical expressions for reactant coverages and the net reaction rates concerning gas phase composition. In particular, these expressions are reached by making the fundamental assumption that both diffusion and adsorption/desorption events are quasi-equilibrated (being the later one the rate-limiting factor), occurring on a much faster timescale than the oxidation event. We use the same argument to speed up our calculations. Thus, we can use the analytical expressions we obtained to assess the overall quality of our results. Our results should be indistinguishable from the Langmuir-Hinshelwood deterministic equations, which are shown below:

\[\begin{split}\begin{gather} \theta_\text{O} = \frac{ \sqrt{B_{\text{O}_2}x_{\text{O}_2}} }{ 1 + B_\text{CO}x_\text{CO} + \sqrt{B_{\text{O}_2}x_{\text{O}_2}} } \\ \theta_\text{CO} = \frac{ B_\text{CO}x_\text{CO} }{ 1 + B_\text{CO}x_\text{CO} + \sqrt{B_{\text{O}_2}x_{\text{O}_2}} } \\[7mm] \text{TOF}_{\text{CO}_2} = 6 \, A_\text{oxi}\theta_\text{CO}\theta_\text{O} \end{gather}\end{split}\]

Here \(B_\text{CO}\)/\(B_{\text{O}_2}\) represent the ratio of the adsorption-desorption pre-exponential terms of \(CO\)/\(O_2\) (pe_ratio in Zacros), \(x_\text{CO}\)/\(x_{\text{O}_2}\) the molar fractions of \(CO\)/\(O_2\); \(\theta_\text{CO}\)/\(\theta_\text{O}\) the coverage of \(CO*\)/\(O*\), \(A_\text{oxi}\) the pre-exponential factor of the CO oxidation step (pre_expon in Zacros), and TOF\(_{\text{CO}_2}\) the turnover frequency or production rate of \(CO_2\). The number 6 is because, in our lattice, each site has 6 neighbors, so the oxidation event is “replicated” across each neighboring site.

Notice that to get the above expressions based on the ones shown in the Zacros tutorial, you need the following equalities:

\[\begin{split}\begin{gather} K_s P_s = \left( \frac{A^\text{ads}_s P^{-1}}{A^\text{des}_s} \right) x_s P = B_s x_s \qquad \therefore\qquad s=\text{CO},\text{O}_2 \\ k_\text{oxi} = A_\text{oxi} \end{gather}\end{split}\]

The final equality follows from the fact that in pyZacros, the activation energy for all elementary reactions is equal to zero.

The code below simply computes the coverages and TOF of \(CO_2\) using the analytical expression described above:

import numpy

lh = pz.models.LangmuirHinshelwood()

B_CO = lh.mechanism.find_one( 'CO_adsorption' ).pe_ratio
B_O2 = lh.mechanism.find_one( 'O2_adsorption' ).pe_ratio
A_oxi = lh.mechanism.find_one( 'CO_oxidation' ).pre_expon

x_CO_model = numpy.linspace(0.0,1.0,201)

ac_O_model = []
ac_CO_model = []
TOF_CO2_model = []

for i in range(len(x_CO_model)):
    x_O2 = 1 - x_CO_model[i]
    ac_O_model.append( numpy.sqrt(B_O2*x_O2)/( 1 + B_CO*x_CO_model[i] + numpy.sqrt(B_O2*x_O2) ) )
    ac_CO_model.append( B_CO*x_CO_model[i]/( 1 + B_CO*x_CO_model[i] + numpy.sqrt(B_O2*x_O2) ) )
    TOF_CO2_model.append( 6*A_oxi*ac_CO_model[i]*ac_O_model[i] )

Additionally, if you have installed the package matplotlib, you can see the aforementioned results visually. Please look over the code below, and notice we plot the analytical and simulation results together. The points in the figure represent simulation results, while the lines represent analytical model results. They are nearly identical.

import matplotlib.pyplot as plt

x_CO_max = (B_O2/B_CO**2/2.0)*(numpy.sqrt(1.0+4.0*B_CO**2/B_O2)-1.0)

fig = plt.figure()

ax = plt.axes()
ax.set_xlim([0.0, 1.0])
ax.set_xlabel('Molar Fraction CO', fontsize=14)
ax.set_ylabel('Coverage Fraction (%)', color='blue', fontsize=14)
ax.vlines( x_CO_max, 0, max(ac_O_model), colors='0.8', linestyles='--',)
ax.plot(x_CO_model, ac_O_model, color='blue', linestyle='-.', lw=2, zorder=1)
ax.plot(x_CO, ac_O, marker='$\u25CF$', color='blue', lw=0, markersize=4, zorder=2)
ax.plot(x_CO_model, ac_CO_model, color='blue', linestyle='-', lw=2, zorder=3)
ax.plot(x_CO, ac_CO, marker='$\u25EF$', color='blue', markersize=4, lw=0, zorder=4)
plt.text(0.3, 0.60, 'O*', fontsize=18, color='blue')
plt.text(0.7, 0.45, 'CO*', fontsize=18, color='blue')

ax2 = ax.twinx()
ax2.set_ylabel('TOF (mol/s/site)',color='red', fontsize=14)
ax2.plot(x_CO_model, TOF_CO2_model, color='red', linestyle='-', lw=2, zorder=5)
ax2.plot(x_CO, TOF_CO2, marker='$\u25EF$', color='red', markersize=4, lw=0, zorder=6)
plt.text(0.3, 200.0, 'CO$_2$', fontsize=18, color='red')

plt.show()
/scm-uploads/doc/pyzacros/_images/CoveragesAndReactionRate_ViewResults_10_0.png

As a final note, we included in the code above the value of the \(CO\) molar fraction (\(x_\text{CO}^*\)) on which we get the maximum \(CO_2\) production rate. The figure shows this value as a vertical gray dashed line. It is simple to deduce it from the preceding analytical expressions:

\[x_{CO}^* = \frac{B_{\text{O}_2}}{2 B_\text{CO}^2}\left( \sqrt{1+\frac{4 B_\text{CO}^2}{B_{\text{O}_2}}} - 1 \right)\approx 0.656\]

Notice that the position of this maximum \(x_\text{CO}*\) depends exclusively on the ratio of the pe_ratio parameters for \(CO\) and \(O_2\) in the Langmuir-Hinshelwood model.