Langmuir-Hinshelwood Model: Acceleration by Automated Rescaling of the Rate Constants
*************************************************************************************

Note:

  To follow this tutorial, either:

  * Download "CoveragesAndReactionRate.py" (run as "$AMSBIN/amspython
    CoveragesAndReactionRate.py").

  * Download "CoveragesAndReactionRate.ipynb" (see also: how to
    install Jupyterlab)

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:

  * Download "CoveragesAndReactionRate_ViewResults.py" (run as
    "$AMSBIN/amspython CoveragesAndReactionRate.py").

  * Download "CoveragesAndReactionRate_ViewResults.ipynb" (see also:
    how to install Jupyterlab)

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()

[image]

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:

   \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}

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:

   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}

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()

[image]

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.
