Ziff-Gulari-Barshad model: Phase Transitions under Steady State Conditions.¶
Note
To follow this tutorial, either:
Download
PhaseTransitions-SteadyState.py(run as$AMSBIN/amspython PhaseTransitions-SteadyState.py).Download
PhaseTransitions-SteadyState.ipynb(see also: how to install Jupyterlab)
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()
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