Ziff-Gulari-Barshad model: Phase Transitions under Steady State Conditions.

Note

To follow this tutorial, either:

This example is inspired in the seminal paper: Kinetic Phase Transitions in an Irreversible Surface-Reaction Model by Robert M. Ziff, Erdagon Gulari, and Yoav Barshad in 1986 Phys. Rev. Lett. 56, (1986) 2553. The authors proposed a simple model for catalytic reactions of carbon monoxide oxidation to carbon dioxide on a surface. This model is now known as the Ziff-Gulari-Barshad (ZGB) model after their names. While the model leaves out many important steps of the real system, it exhibits interesting steady-state off-equilibrium behavior and two types of phase transitions, which actually occur in real systems. Please refer to the original paper for more details. In this example, we will analyze the effect of changing the composition of the gas phase, namely partial pressures for \(O_2\) and \(CO\), in the \(CO_2\) Turnover frequency (TOF) in the ZGB model. At variance with the example Phase Transitions in the ZGB model, here we extend the dynamics up to reaching the system’s steady state for each gas phase composition value. Most of the code is the same, except for the section regarding the setup of the steady state calculation.

First, we 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

On a typical laptop, this calculation should take no more than 10 min to complete. Here we illustrate how to run several parallel Zacros calculations. We’ll use the plams.JobRunner class, which easily allows us to run as many parallel instances as we request. In this case, we choose to use the maximum number of simultaneous processes (maxjobs) equal to the number of processors in the machine. Additionally, by setting nproc =  1 we establish that only one processor will be used for each zacros instance.

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

Running up to 8 jobs in parallel simultaneously

Now, we initialize our Ziff-Gulari-Barshad model, which by luck is available as a predefined model in pyZacros,

zgb = pz.models.ZiffGulariBarshad()

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 Ziff-Gulari-Barshad 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.0
z_sett.max_time = 10.0
z_sett.species_numbers = ('time', 0.1)
z_sett.random_seed = 953129

z_job = pz.ZacrosJob( settings=z_sett, lattice=zgb.lattice,
                      mechanism=zgb.mechanism,
                      cluster_expansion=zgb.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 96% 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. 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.nbatch = 20
ss_sett.turnover_frequency.confidence = 0.96
ss_sett.turnover_frequency.nreplicas = 4

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 )

[02.02|22:26:23] JOB plamsjob Steady State Convergence: Using nbatch=20,confidence=0.96,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.2 to 0.8 in steps of 0.01. 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.arange(0.2, 0.8, 0.01) )
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 )

[02.02|22:26:23] JOB ps_cond000 Steady State Convergence: Using nbatch=20,confidence=0.96,ignore_nbatch=1,nreplicas=4
[02.02|22:26:23] JOB ps_cond001 Steady State Convergence: Using nbatch=20,confidence=0.96,ignore_nbatch=1,nreplicas=4
[02.02|22:26:23] JOB ps_cond002 Steady State Convergence: Using nbatch=20,confidence=0.96,ignore_nbatch=1,nreplicas=4
[02.02|22:26:23] JOB ps_cond003 Steady State Convergence: Using nbatch=20,confidence=0.96,ignore_nbatch=1,nreplicas=4
[02.02|22:26:23] JOB ps_cond004 Steady State Convergence: Using nbatch=20,confidence=0.96,ignore_nbatch=1,nreplicas=4

...

[02.02|22:26:23] JOB ps_cond056 Steady State Convergence: Using nbatch=20,confidence=0.96,ignore_nbatch=1,nreplicas=4
[02.02|22:26:23] JOB ps_cond057 Steady State Convergence: Using nbatch=20,confidence=0.96,ignore_nbatch=1,nreplicas=4
[02.02|22:26:23] JOB ps_cond058 Steady State Convergence: Using nbatch=20,confidence=0.96,ignore_nbatch=1,nreplicas=4
[02.02|22:26:23] JOB ps_cond059 Steady State Convergence: Using nbatch=20,confidence=0.96,ignore_nbatch=1,nreplicas=4
[02.02|22:26:23] JOB ps_cond060 Steady State Convergence: Using nbatch=20,confidence=0.96,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. This step should take less than 10 mins!

results = ps_job.run()

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

[02.02|22:26:23] JOB plamsjob STARTED
[02.02|22:26:23] Waiting for job plamsjob to finish
[02.02|22:26:23] JOB plamsjob RUNNING
[02.02|22:26:23] JOB plamsjob/ps_cond000 STARTED
[02.02|22:26:23] JOB plamsjob/ps_cond001 STARTED
[02.02|22:26:23] JOB plamsjob/ps_cond002 STARTED

...

[02.02|22:29:57]           CO2        0.44695        0.02674        0.05983     False
[02.02|22:29:57] JOB plamsjob/ps_cond020/ss_iter004_rep000 FINISHED
[02.02|22:29:57] JOB plamsjob/ps_cond020/ss_iter004_rep000 FAILED
[02.02|22:29:57] Waiting for job ss_iter004_rep001 to finish
[02.02|22:29:57] JOB plamsjob/ps_cond020/ss_iter004_rep001 FINISHED
[02.02|22:29:58] JOB plamsjob/ps_cond020/ss_iter004_rep002 FINISHED
[02.02|22:29:58]    Average
[02.02|22:29:58] JOB plamsjob/ps_cond020/ss_iter004_rep001 SUCCESSFUL
[02.02|22:29:58] Waiting for job ss_iter004_rep003 to finish
[02.02|22:29:58]       species            TOF          error          ratio     conv?
[02.02|22:29:58]            CO       -0.42282        0.01585        0.03748      True
[02.02|22:29:58]            O2       -0.21130        0.00793        0.03751      True
[02.02|22:29:58]           CO2        0.42282        0.01585        0.03748      True
[02.02|22:29:58] JOB plamsjob/ps_cond021 Steady State Convergence: CONVERGENCE REACHED. DONE!
[02.02|22:29:58] JOB plamsjob/ps_cond020/ss_iter004_rep002 SUCCESSFUL
[02.02|22:29:58] JOB plamsjob/ps_cond020/ss_iter004_rep003 FINISHED
[02.02|22:29:58] JOB plamsjob/ps_cond021 FINISHED
[02.02|22:29:59] JOB plamsjob/ps_cond020/ss_iter004_rep003 SUCCESSFUL
[02.02|22:29:59] WARNING: Trying to obtain results of crashed or failed job ss_iter004_rep000
[02.02|22:29:59] WARNING: Trying to obtain results of crashed or failed job ss_iter004_rep000
[02.02|22:29:59] Obtaining results of ss_iter004_rep000 successful. However, no guarantee that they make sense
[02.02|22:29:59] Obtaining results of ss_iter004_rep000 successful. However, no guarantee that they make sense
[02.02|22:29:59] JOB plamsjob/ps_cond020/ss_iter004_rep000 Steady State Convergence: RESTART ABORTED
[02.02|22:29:59] JOB plamsjob/ps_cond020/ss_iter004_rep000 Steady State Convergence: JOB REMOVED
[02.02|22:29:59] JOB plamsjob/ps_cond020/ss_iter004_rep001 Steady State Convergence: JOB REMOVED
[02.02|22:29:59] JOB plamsjob/ps_cond020/ss_iter004_rep002 Steady State Convergence: JOB REMOVED
[02.02|22:29:59] JOB plamsjob/ps_cond020/ss_iter004_rep003 Steady State Convergence: JOB REMOVED
[02.02|22:30:00] JOB plamsjob/ps_cond020 FINISHED
[02.02|22:30:00] JOB plamsjob/ps_cond021 SUCCESSFUL
[02.02|22:30:00] JOB plamsjob/ps_cond020 SUCCESSFUL
[02.02|22:30:03] JOB plamsjob FINISHED
[02.02|22:30:14] 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 10 steps in the simulation to calculate the average coverages. Notably, the loop iterates over the internal indices (results.indices()) of each condition. i is simply a sequential number for each condition in this case, but idx is an id for the corresponding child job. This id may be a complex object for n-dimension (n>1) scanning parameters, so if you want to access the properties of one of the child jobs, we recommend using a loop like the one we use here. In the lines that follow, we use this idx to get the maximum time that the simulation required to achieve the steady state for that specific composition (“max time”).

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

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

for i,idx in enumerate(results.indices()):
    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'] )
    max_time.append( results.children_results( child_id=idx ).history( pos=-1 )['max_time'] )

print( "-----------------------------------------------------------" )
print( "%4s"%"cond", "%8s"%"x_CO", "%10s"%"ac_O", "%10s"%"ac_CO", "%12s"%"TOF_CO2", "%10s"%"max_time" )
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], "%10.3f"%max_time[i] )

-----------------------------------------------------------
cond     x_CO       ac_O      ac_CO      TOF_CO2   max_time
-----------------------------------------------------------
   0     0.20   0.999250   0.000000     0.002567     20.000
   1     0.21   0.999900   0.000000     0.000000     20.000
   2     0.22   0.999590   0.000000     0.000700     20.000
   3     0.23   0.999130   0.000000     0.002100     20.000
   4     0.24   0.998630   0.000000     0.004267     20.000
   5     0.25   0.998720   0.000000     0.003267     20.000
   6     0.26   0.998760   0.000000     0.003267     20.000
   7     0.27   0.998920   0.000000     0.003233     20.000
   8     0.28   0.998090   0.000000     0.007100     20.000
   9     0.29   1.000000   0.000000     0.003241     80.000
  10     0.30   1.000000   0.000000     0.004973     80.000
  11     0.31   1.000000   0.000000     0.005157     80.000
  12     0.32   0.992570   0.000020     0.026933     20.000
  13     0.33   1.000000   0.000000     0.007310     80.000
  14     0.34   1.000000   0.000000     0.011559     80.000
  15     0.35   1.000000   0.000000     0.016751     80.000
  16     0.36   1.000000   0.000000     0.022879     80.000
  17     0.37   1.000000   0.000000     0.041531     80.000
  18     0.38   0.999260   0.000000     0.069112     80.000
  19     0.39   0.985810   0.000140     0.146791     80.000
  20     0.40   0.975090   0.000140     0.214842    320.000
  21     0.41   0.905970   0.000780     0.422818    320.000
  22     0.42   0.856420   0.001410     0.613742     80.000
  23     0.43   0.806960   0.002770     0.795856     80.000
  24     0.44   0.778850   0.003330     0.975367     20.000
  25     0.45   0.739230   0.004710     1.117933     80.000
  26     0.46   0.728010   0.005240     1.278163     80.000
  27     0.47   0.688730   0.008080     1.455912     80.000
  28     0.48   0.624570   0.013870     1.625328     80.000
  29     0.49   0.607870   0.014540     1.831873     80.000
  30     0.50   0.568050   0.019630     2.021464     80.000
  31     0.51   0.532990   0.034100     2.234967     20.000
  32     0.52   0.472360   0.050070     2.461300     20.000
  33     0.53   0.270100   0.359930     2.502505     80.000
  34     0.54   0.000000   1.000000     0.191716     80.000
  35     0.55   0.000000   1.000000     0.066715     80.000
  36     0.56   0.000000   1.000000     0.000000     20.000
  37     0.57   0.000000   1.000000     0.000000     20.000
  38     0.58   0.000000   1.000000     0.000000     20.000
  39     0.59   0.000000   1.000000     0.000000     20.000
  40     0.60   0.000000   1.000000     0.000000     20.000
  41     0.61   0.000000   1.000000     0.000000     20.000
  42     0.62   0.000000   1.000000     0.000000     20.000
  43     0.63   0.000000   1.000000     0.000000     20.000
  44     0.64   0.000000   1.000000     0.000000     20.000
  45     0.65   0.000000   1.000000     0.000000     20.000
  46     0.66   0.000000   1.000000     0.000000     20.000
  47     0.67   0.000000   1.000000    -0.000000     20.000
  48     0.68   0.000000   1.000000     0.000000     20.000
  49     0.69   0.000000   1.000000     0.000000     20.000
  50     0.70   0.000000   1.000000     0.000000     20.000
  51     0.71   0.000000   1.000000    -0.000000     20.000
  52     0.72   0.000000   1.000000    -0.000000     20.000
  53     0.73   0.000000   1.000000     0.000000     20.000
  54     0.74   0.000000   1.000000     0.000000     20.000
  55     0.75   0.000000   1.000000     0.000000     20.000
  56     0.76   0.000000   1.000000     0.000000     20.000
  57     0.77   0.000000   1.000000    -0.000000     20.000
  58     0.78   0.000000   1.000000     0.000000     20.000
  59     0.79   0.000000   1.000000     0.000000     20.000
  60     0.80   0.000000   1.000000     0.000000     20.000

The above results are the final aim of the calculation. However, we can take advantage of python libraries to visualize them. Here, we use matplotlib. Please check the matplotlib documentation for more details at matplotlib. The following lines of code allow visualizing the effect of changing the \(CO\) molar fraction on the average coverage of \(O*\) and \(CO*\) and the production rate of \(CO_2\):

import matplotlib.pyplot as plt

fig = plt.figure()

ax = plt.axes()
ax.set_xlabel('Molar fraction CO', fontsize=14)
ax.set_ylabel("Coverage Fraction (%)", color="blue", fontsize=14)
ax.plot(x_CO, ac_O, color="blue", linestyle="-.", lw=2, zorder=1)
ax.plot(x_CO, ac_CO, color="blue", linestyle="-", lw=2, zorder=2)
plt.text(0.3, 0.9, 'O', fontsize=18, color="blue")
plt.text(0.7, 0.9, '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, color="red", lw=2, zorder=5)
plt.text(0.37, 1.5, 'CO$_2$', fontsize=18, color="red")

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

As shown in the Figure, we also found three regions involving two phase transitions, which are similar to the ones obtained in the example Phase Transitions in the ZGB model, where we did not converge the calculation up to the steady-state condition. However, in this case, the first transition, corresponding to the limit of \(O^*\) poisoning, occurs at x_CO=0.40 rather than x_CO=0.32, and the second transition, corresponding to the beginning of CO poisoning, occurs roughly at the same value x_CO=0.55 but as a more abrupt change.

Now, we can close the pyZacros environment:

scm.pyzacros.finish()

[02.02|22:31:38] PLAMS run finished. Goodbye